[cig-commits] r11502 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . UTILS UTILS/estimate_best_values_runs
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Mar 22 08:45:43 PDT 2008
Author: dkomati1
Date: 2008-03-22 08:45:42 -0700 (Sat, 22 Mar 2008)
New Revision: 11502
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/sort_array_coordinates.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt
Log:
- added Cuthill-McKee sorting (debugged by David Michea last week)
- suppressed saving iMPIcut_xi and iMPIcut_eta to disk (which had become useless)
- added directory "estimate_best_values_runs" in UTILS with a program to determine the best sets of parameters to maximize memory use on a given machine
- fixed a bug in the Makefile
- updated the todo list
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2008-03-22 15:45:42 UTC (rev 11502)
@@ -83,7 +83,6 @@
$O/define_derivation_matrices.o \
$O/define_superbrick.o \
$O/euler_angles.o \
- $O/exit_mpi.o \
$O/get_MPI_1D_buffers.o \
$O/get_MPI_cutplanes_eta.o \
$O/get_MPI_cutplanes_xi.o \
@@ -96,6 +95,7 @@
$O/get_jacobian_boundaries.o \
$O/get_jacobian_discontinuities.o \
$O/get_model.o \
+ $O/get_perm_cuthill_mckee.o \
$O/get_shape2D.o \
$O/get_shape3D.o \
$O/get_value_parameters.o \
@@ -193,13 +193,13 @@
####
# rules for the main programs
-XMESHFEM_OBJECTS = $O/meshfem3D.o $(LIBSPECFEM)
+XMESHFEM_OBJECTS = $O/meshfem3D.o $O/exit_mpi.o $(LIBSPECFEM)
xmeshfem3D: $(XMESHFEM_OBJECTS)
## use MPI here
${MPIFCCOMPILE_CHECK} -o xmeshfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
# solver also depends on values from mesher
-XSPECFEM_OBJECTS = $(SOLVER_ARRAY_OBJECTS) $(LIBSPECFEM)
+XSPECFEM_OBJECTS = $(SOLVER_ARRAY_OBJECTS) $O/exit_mpi.o $(LIBSPECFEM)
xspecfem3D: $(XSPECFEM_OBJECTS)
## use MPI here
${MPIFCCOMPILE_NO_CHECK} -o xspecfem3D $(XSPECFEM_OBJECTS) $(MPILIBS)
@@ -538,6 +538,9 @@
$O/create_serial_name_database.o: constants.h $S/create_serial_name_database.f90
${FCCOMPILE_CHECK} -c -o $O/create_serial_name_database.o ${FCFLAGS_f90} $S/create_serial_name_database.f90
+$O/get_perm_cuthill_mckee.o: constants.h $S/get_perm_cuthill_mckee.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_perm_cuthill_mckee.o ${FCFLAGS_f90} $S/get_perm_cuthill_mckee.f90
+
## use MPI here
$O/read_arrays_buffers_solver.o: constants.h $S/read_arrays_buffers_solver.f90
${MPIFCCOMPILE_CHECK} -c -o $O/read_arrays_buffers_solver.o ${FCFLAGS_f90} $S/read_arrays_buffers_solver.f90
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -0,0 +1,273 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! create file OUTPUT_FILES/values_from_mesher.h based upon DATA/Par_file
+! in order to compile the solver with the right array sizes
+
+ subroutine create_header_file
+
+ implicit none
+
+ include "constants.h"
+
+! parameters read from parameter file
+ integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NSOURCES,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+
+ double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH
+
+ logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
+ TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D, &
+ RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
+ SAVE_MESH_FILES,ATTENUATION,CASE_3D, &
+ ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,first_time
+
+ character(len=150) LOCAL_PATH,MODEL
+
+! parameters deduced from parameters read from file
+ integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+
+! this for all the regions
+ integer, dimension(MAX_NUM_REGIONS) :: NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ nglob
+
+ double precision :: static_memory_size
+ character(len=150) HEADER_FILE
+
+ integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION
+
+ integer :: iregion
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
+ integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
+ integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
+ integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NGLOB1D_RADIAL_CORNER
+ integer, dimension(MAX_NUM_REGIONS) :: NGLOB1D_RADIAL_TEMP
+
+ integer :: c,compteur
+
+ double precision :: mem_per_core,percent
+
+ integer, parameter :: NB_COLONNES = 10
+ integer, dimension(NB_COLONNES) :: store_NEX_XI
+
+! maximum total number of processors we want to see in the table
+ integer, parameter :: MAX_NUMBER_OF_PROCS = 62976 !! Ranger !! 100000
+! integer, parameter :: MAX_NUMBER_OF_PROCS = 212992 ! current maximum on BlueGene at LLNL
+
+! amount of memory installed per processor core on the system (in gigabytes)
+ double precision, parameter :: MAX_MEMORY_PER_CORE = 2.d0 !! Ranger
+
+! base value depends if we implement three or four doublings (default is three)
+ integer, parameter :: NB_DOUBLING = 3
+ integer :: BASE_VALUE
+
+! output in LaTeX format or in regular ASCII format
+ logical, parameter :: OUTPUT_LATEX_FORMAT = .false.
+
+! ************** PROGRAM STARTS HERE **************
+
+ call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
+
+! count the total number of sources in the CMTSOLUTION file
+ NSOURCES = 1
+
+ BASE_VALUE = 2**NB_DOUBLING
+
+ first_time = .true.
+
+ print *
+ print *,'Number of GB of memory per core on the machine: ',MAX_MEMORY_PER_CORE
+ print *
+
+ print *,'total_proc, % of Ranger, NPROC_XI, NEX_XI, memory used per core, percentage:'
+ print *
+
+! total number of processors is 6 * NPROC_XI^2
+!!!! do NPROC_XI = 1,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
+!! start after 1000 processors
+ do NPROC_XI = 13,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
+
+ compteur = 1
+ c = 0
+
+ do while (compteur <= NB_COLONNES)
+ c = c + 1
+ NEX_XI = BASE_VALUE * c * NPROC_XI
+
+ if(NEX_XI >= 64 .and. mod(NEX_XI,2*BASE_VALUE) == 0) then
+
+ store_NEX_XI(compteur) = NEX_XI
+ compteur = compteur + 1
+
+! read the parameter file and compute additional parameters
+ call read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
+ ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
+ MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
+ PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
+ ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB, &
+ ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+
+ do iregion=1,MAX_NUM_REGIONS
+ NGLOB1D_RADIAL_CORNER(iregion,:) = NGLOB1D_RADIAL(iregion)
+ enddo
+
+ if (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA) then
+ NGLOB1D_RADIAL_CORNER(IREGION_OUTER_CORE,:) = NGLOB1D_RADIAL_CORNER(IREGION_OUTER_CORE,:) + &
+ maxval(DIFF_NSPEC1D_RADIAL(:,:))*(NGLLZ-1)
+ endif
+
+ if(first_time) then
+ first_time = .false.
+ if(ATTENUATION) then
+ print *,'ATTENUATION = .true.'
+ else
+ print *,'ATTENUATION = .false.'
+ endif
+ print *
+ endif
+
+!! compute memory usage per processor core
+! evaluate the amount of static memory needed by the solver
+ call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
+ ONE_CRUST,doubling_index,this_region_has_a_doubling,&
+ ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
+ NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
+ NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+
+ mem_per_core = static_memory_size/1073741824.d0
+
+ percent = 100.d0 * mem_per_core / MAX_MEMORY_PER_CORE
+
+ if(percent > 100.d0) goto 777
+
+ if(percent < 0.d0) then
+ goto 777
+! write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **mesher fails/temporary bug**')") &
+! 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+
+ else if(percent >= 93.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **too high**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+ goto 777
+
+ else if(percent >= 85.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **excellent**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+
+ else if(percent >= 80.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **very good**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+
+ else
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'%')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+ endif
+
+ endif
+ enddo
+
+ 777 continue
+
+ print *
+
+ enddo
+
+ end subroutine create_header_file
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -0,0 +1,2338 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
+ ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
+ MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
+ PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
+ ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ NGLOB, &
+ ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+
+
+ implicit none
+
+ include "constants.h"
+
+ integer NEX_ETA_dummy , NEX_XI_dummy, NPROC_ETA_dummy , NPROC_XI_dummy
+
+! parameters read from parameter file
+ integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+
+ double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
+ MOVIE_TOP_KM,MOVIE_TOP,MOVIE_BOTTOM_KM,MOVIE_BOTTOM, &
+ MOVIE_EAST_DEG,MOVIE_EAST,MOVIE_WEST_DEG,MOVIE_WEST,MOVIE_NORTH_DEG,MOVIE_NORTH,MOVIE_SOUTH_DEG,MOVIE_SOUTH
+
+ logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
+ TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D, &
+ RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
+ SAVE_MESH_FILES,ATTENUATION, &
+ ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,WRITE_SEISMOGRAMS_BY_MASTER,&
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE
+
+ character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
+
+! local variables
+ integer NEX_MAX
+
+ double precision RECORD_LENGTH_IN_MINUTES,ELEMENT_WIDTH
+
+ integer, external :: err_occurred
+
+! parameters to be computed based upon parameters above read from file
+ integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
+
+ integer, dimension(MAX_NUM_REGIONS) :: NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ NGLOB
+
+ integer nblocks_xi,nblocks_eta
+
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
+
+ integer :: ielem,elem_doubling_mantle,elem_doubling_middle_outer_core,elem_doubling_bottom_outer_core
+ double precision :: DEPTH_SECOND_DOUBLING_REAL,DEPTH_THIRD_DOUBLING_REAL, &
+ DEPTH_FOURTH_DOUBLING_REAL,distance,distance_min,zval
+
+! honor PREM Moho or not
+! doing so drastically reduces the stability condition and therefore the time step
+ logical :: HONOR_1D_SPHERICAL_MOHO,CASE_3D
+
+ integer :: ifirst_region, ilast_region, iter_region, iter_layer, doubling, padding, tmp_sum, tmp_sum_xi, tmp_sum_eta
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+ integer :: NUMBER_OF_MESH_LAYERS,layer_offset,nspec2D_xi_sb,nspec2D_eta_sb, &
+ nb_lay_sb, nspec_sb, nglob_vol, nglob_surf, nglob_edge
+
+ integer :: multiplication_factor
+
+! for the cut doublingbrick improvement
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
+ integer :: lastdoubling_layer, cut_doubling, nglob_int_surf_xi, nglob_int_surf_eta,nglob_ext_surf,&
+ normal_doubling, nglob_center_edge, nglob_corner_edge, nglob_border_edge
+ integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
+ integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
+
+! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
+ call open_parameter_file
+
+ call read_value_integer(SIMULATION_TYPE, 'solver.SIMULATION_TYPE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(SAVE_FORWARD, 'solver.SAVE_FORWARD')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ call read_value_integer(NCHUNKS, 'mesher.NCHUNKS')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
+ stop 'NCHUNKS must be either 1, 2, 3 or 6'
+
+ call read_value_double_precision(ANGULAR_WIDTH_XI_IN_DEGREES, 'mesher.ANGULAR_WIDTH_XI_IN_DEGREES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(ANGULAR_WIDTH_ETA_IN_DEGREES, 'mesher.ANGULAR_WIDTH_ETA_IN_DEGREES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(CENTER_LATITUDE_IN_DEGREES, 'mesher.CENTER_LATITUDE_IN_DEGREES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(CENTER_LONGITUDE_IN_DEGREES, 'mesher.CENTER_LONGITUDE_IN_DEGREES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(GAMMA_ROTATION_AZIMUTH, 'mesher.GAMMA_ROTATION_AZIMUTH')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+! this MUST be 90 degrees for two chunks or more to match geometrically
+ if(NCHUNKS > 1 .and. abs(ANGULAR_WIDTH_XI_IN_DEGREES - 90.d0) > 0.00000001d0) &
+ stop 'ANGULAR_WIDTH_XI_IN_DEGREES must be 90 for more than one chunk'
+
+! this can be any value in the case of two chunks
+ if(NCHUNKS > 2 .and. abs(ANGULAR_WIDTH_ETA_IN_DEGREES - 90.d0) > 0.00000001d0) &
+ stop 'ANGULAR_WIDTH_ETA_IN_DEGREES must be 90 for more than two chunks'
+
+! include central cube or not
+! use regular cubed sphere instead of cube for large distances
+ if(NCHUNKS == 6) then
+ INCLUDE_CENTRAL_CUBE = .true.
+ INFLATE_CENTRAL_CUBE = .false.
+ else
+ INCLUDE_CENTRAL_CUBE = .false.
+ INFLATE_CENTRAL_CUBE = .true.
+ endif
+
+! number of elements at the surface along the two sides of the first chunk
+ call read_value_integer(NEX_XI_dummy, 'mesher.NEX_XI')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NEX_ETA_dummy, 'mesher.NEX_ETA')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NPROC_XI_dummy, 'mesher.NPROC_XI')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NPROC_ETA_dummy, 'mesher.NPROC_ETA')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+!! DK DK
+ NEX_ETA = NEX_XI
+ NPROC_ETA = NPROC_XI
+
+! define the velocity model
+ call read_value_string(MODEL, 'model.name')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+! use PREM as the 1D reference model by default
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
+
+! HONOR_1D_SPHERICAL_MOHO: honor PREM Moho or not: doing so drastically reduces
+! the stability condition and therefore the time step, resulting in expensive
+! calculations. If not, honor a fictitious Moho at the depth of 40 km
+! in order to have even radial sampling from the d220 to the Earth surface.
+
+! ONE_CRUST: in order to increase stability and therefore to allow cheaper
+! simulations (larger time step), 1D models can be run with just one average crustal
+! layer instead of two.
+
+! CASE_3D : this flag allows the stretching of the elements in the crustal
+! layers in the case of 3D models. The purpose of this stretching is to squeeze more
+! GLL points per km in the upper part of the crust than in the lower part.
+ HONOR_1D_SPHERICAL_MOHO = .false.
+ ONE_CRUST = .false.
+ CASE_3D = .false.
+
+! default is no 3D model
+ THREE_D_MODEL = 0
+
+ if(MODEL == '1D_isotropic_prem') then
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+
+ else if(MODEL == '1D_transversely_isotropic_prem') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+
+ else if(MODEL == '1D_iasp91' .or. MODEL == '1D_1066a' .or. MODEL == '1D_ak135') then
+ if(MODEL == '1D_iasp91') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_IASP91
+ else if(MODEL == '1D_1066a') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_1066A
+ else if(MODEL == '1D_ak135') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
+ else
+ stop 'reference 1D Earth model unknown'
+ endif
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+
+ else if(MODEL == '1D_ref') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+
+ else if(MODEL == '1D_ref_iso') then
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+
+ else if(MODEL == '1D_isotropic_prem_onecrust') then
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ ONE_CRUST = .true.
+
+ else if(MODEL == '1D_transversely_isotropic_prem_onecrust') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ ONE_CRUST = .true.
+
+ else if(MODEL == '1D_iasp91_onecrust' .or. MODEL == '1D_1066a_onecrust' .or. MODEL == '1D_ak135_onecrust') then
+ if(MODEL == '1D_iasp91_onecrust') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_IASP91
+ else if(MODEL == '1D_1066a_onecrust') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_1066A
+ else if(MODEL == '1D_ak135_onecrust') then
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
+ else
+ stop 'reference 1D Earth model unknown'
+ endif
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ ONE_CRUST = .true.
+
+ else if(MODEL == 'transversely_isotropic_prem_plus_3D_crust_2.0') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+
+ else if(MODEL == 's20rts') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
+ THREE_D_MODEL = THREE_D_MODEL_S20RTS
+
+ else if(MODEL == 's362ani') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+ THREE_D_MODEL = THREE_D_MODEL_S362ANI
+
+ else if(MODEL == 's362iso') then
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+ THREE_D_MODEL = THREE_D_MODEL_S362ANI
+
+ else if(MODEL == 's362wmani') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+ THREE_D_MODEL = THREE_D_MODEL_S362WMANI
+
+ else if(MODEL == 's362ani_prem') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
+ THREE_D_MODEL = THREE_D_MODEL_S362ANI_PREM
+
+ else if(MODEL == 's29ea') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .true.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
+ THREE_D_MODEL = THREE_D_MODEL_S29EA
+
+ else if(MODEL == '3D_attenuation') then
+ TRANSVERSE_ISOTROPY = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .true.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+
+ else if(MODEL == '3D_anisotropic') then
+ TRANSVERSE_ISOTROPY = .true.
+ ISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_INNER_CORE = .false.
+ CRUSTAL = .false.
+ ATTENUATION_3D = .false.
+ ONE_CRUST = .true.
+ CASE_3D = .true.
+
+ else
+ stop 'model not implemented, edit read_compute_parameters.f90 and recompile'
+ endif
+
+! set time step, radial distribution of elements, and attenuation period range
+! right distribution is determined based upon maximum value of NEX
+ NEX_MAX = max(NEX_XI,NEX_ETA)
+
+!----
+!---- case prem_onecrust by default
+!----
+ if (SUPPRESS_CRUSTAL_MESH) then
+ multiplication_factor=2
+ else
+ multiplication_factor=1
+ endif
+
+ ! element width = 0.5625000 degrees = 62.54715 km
+ if(NEX_MAX*multiplication_factor <= 160) then
+ DT = 0.252d0
+
+ MIN_ATTENUATION_PERIOD = 30
+ MAX_ATTENUATION_PERIOD = 1500
+
+ NER_CRUST = 1
+ NER_80_MOHO = 1
+ NER_220_80 = 2
+ NER_400_220 = 2
+ NER_600_400 = 2
+ NER_670_600 = 1
+ NER_771_670 = 1
+ NER_TOPDDOUBLEPRIME_771 = 15
+ NER_CMB_TOPDDOUBLEPRIME = 1
+ NER_OUTER_CORE = 16
+ NER_TOP_CENTRAL_CUBE_ICB = 2
+ R_CENTRAL_CUBE = 950000.d0
+
+ ! element width = 0.3515625 degrees = 39.09196 km
+ else if(NEX_MAX*multiplication_factor <= 256) then
+ DT = 0.225d0
+
+ MIN_ATTENUATION_PERIOD = 20
+ MAX_ATTENUATION_PERIOD = 1000
+
+ NER_CRUST = 1
+ NER_80_MOHO = 1
+ NER_220_80 = 2
+ NER_400_220 = 3
+ NER_600_400 = 3
+ NER_670_600 = 1
+ NER_771_670 = 1
+ NER_TOPDDOUBLEPRIME_771 = 22
+ NER_CMB_TOPDDOUBLEPRIME = 2
+ NER_OUTER_CORE = 24
+ NER_TOP_CENTRAL_CUBE_ICB = 3
+ R_CENTRAL_CUBE = 965000.d0
+
+ ! element width = 0.2812500 degrees = 31.27357 km
+ else if(NEX_MAX*multiplication_factor <= 320) then
+ DT = 0.16d0
+
+ MIN_ATTENUATION_PERIOD = 15
+ MAX_ATTENUATION_PERIOD = 750
+
+ NER_CRUST = 1
+ NER_80_MOHO = 1
+ NER_220_80 = 3
+ NER_400_220 = 4
+ NER_600_400 = 4
+ NER_670_600 = 1
+ NER_771_670 = 2
+ NER_TOPDDOUBLEPRIME_771 = 29
+ NER_CMB_TOPDDOUBLEPRIME = 2
+ NER_OUTER_CORE = 32
+ NER_TOP_CENTRAL_CUBE_ICB = 4
+ R_CENTRAL_CUBE = 940000.d0
+
+ ! element width = 0.1875000 degrees = 20.84905 km
+ else if(NEX_MAX*multiplication_factor <= 480) then
+ DT = 0.11d0
+
+ MIN_ATTENUATION_PERIOD = 10
+ MAX_ATTENUATION_PERIOD = 500
+
+ NER_CRUST = 1
+ NER_80_MOHO = 2
+ NER_220_80 = 4
+ NER_400_220 = 5
+ NER_600_400 = 6
+ NER_670_600 = 2
+ NER_771_670 = 2
+ NER_TOPDDOUBLEPRIME_771 = 44
+ NER_CMB_TOPDDOUBLEPRIME = 3
+ NER_OUTER_CORE = 48
+ NER_TOP_CENTRAL_CUBE_ICB = 5
+ R_CENTRAL_CUBE = 988000.d0
+
+ ! element width = 0.1757812 degrees = 19.54598 km
+ else if(NEX_MAX*multiplication_factor <= 512) then
+ DT = 0.1125d0
+
+ MIN_ATTENUATION_PERIOD = 9
+ MAX_ATTENUATION_PERIOD = 500
+
+ NER_CRUST = 1
+ NER_80_MOHO = 2
+ NER_220_80 = 4
+ NER_400_220 = 6
+ NER_600_400 = 6
+ NER_670_600 = 2
+ NER_771_670 = 3
+ NER_TOPDDOUBLEPRIME_771 = 47
+ NER_CMB_TOPDDOUBLEPRIME = 3
+ NER_OUTER_CORE = 51
+ NER_TOP_CENTRAL_CUBE_ICB = 5
+ R_CENTRAL_CUBE = 1010000.d0
+
+ ! element width = 0.1406250 degrees = 15.63679 km
+ else if(NEX_MAX*multiplication_factor <= 640) then
+ DT = 0.09d0
+
+ MIN_ATTENUATION_PERIOD = 8
+ MAX_ATTENUATION_PERIOD = 400
+
+ NER_CRUST = 2
+ NER_80_MOHO = 3
+ NER_220_80 = 5
+ NER_400_220 = 7
+ NER_600_400 = 8
+ NER_670_600 = 3
+ NER_771_670 = 3
+ NER_TOPDDOUBLEPRIME_771 = 59
+ NER_CMB_TOPDDOUBLEPRIME = 4
+ NER_OUTER_CORE = 64
+ NER_TOP_CENTRAL_CUBE_ICB = 6
+ R_CENTRAL_CUBE = 1020000.d0
+
+ ! element width = 0.1041667 degrees = 11.58280 km
+ else if(NEX_MAX*multiplication_factor <= 864) then
+ DT = 0.0667d0
+
+ MIN_ATTENUATION_PERIOD = 6
+ MAX_ATTENUATION_PERIOD = 300
+
+ NER_CRUST = 2
+ NER_80_MOHO = 4
+ NER_220_80 = 6
+ NER_400_220 = 10
+ NER_600_400 = 10
+ NER_670_600 = 3
+ NER_771_670 = 4
+ NER_TOPDDOUBLEPRIME_771 = 79
+ NER_CMB_TOPDDOUBLEPRIME = 5
+ NER_OUTER_CORE = 86
+ NER_TOP_CENTRAL_CUBE_ICB = 9
+ R_CENTRAL_CUBE = 990000.d0
+
+ ! element width = 7.8125000E-02 degrees = 8.687103 km
+ else if(NEX_MAX*multiplication_factor <= 1152) then
+ DT = 0.05d0
+
+ MIN_ATTENUATION_PERIOD = 4
+ MAX_ATTENUATION_PERIOD = 200
+
+ NER_CRUST = 3
+ NER_80_MOHO = 6
+ NER_220_80 = 8
+ NER_400_220 = 13
+ NER_600_400 = 13
+ NER_670_600 = 4
+ NER_771_670 = 6
+ NER_TOPDDOUBLEPRIME_771 = 106
+ NER_CMB_TOPDDOUBLEPRIME = 7
+ NER_OUTER_CORE = 116
+ NER_TOP_CENTRAL_CUBE_ICB = 12
+ R_CENTRAL_CUBE = 985000.d0
+
+ ! element width = 7.2115384E-02 degrees = 8.018865 km
+ else if(NEX_MAX*multiplication_factor <= 1248) then
+ DT = 0.0462d0
+
+ MIN_ATTENUATION_PERIOD = 4
+ MAX_ATTENUATION_PERIOD = 200
+
+ NER_CRUST = 3
+ NER_80_MOHO = 6
+ NER_220_80 = 9
+ NER_400_220 = 14
+ NER_600_400 = 14
+ NER_670_600 = 5
+ NER_771_670 = 6
+ NER_TOPDDOUBLEPRIME_771 = 114
+ NER_CMB_TOPDDOUBLEPRIME = 8
+ NER_OUTER_CORE = 124
+ NER_TOP_CENTRAL_CUBE_ICB = 13
+ R_CENTRAL_CUBE = 985000.d0
+
+ else
+
+! scale with respect to 1248 if above that limit
+ DT = 0.0462d0 * 1248.d0 / (2.d0*NEX_MAX)
+
+ MIN_ATTENUATION_PERIOD = 4
+ MAX_ATTENUATION_PERIOD = 200
+
+ NER_CRUST = nint(3 * 2.d0*NEX_MAX / 1248.d0)
+ NER_80_MOHO = nint(6 * 2.d0*NEX_MAX / 1248.d0)
+ NER_220_80 = nint(9 * 2.d0*NEX_MAX / 1248.d0)
+ NER_400_220 = nint(14 * 2.d0*NEX_MAX / 1248.d0)
+ NER_600_400 = nint(14 * 2.d0*NEX_MAX / 1248.d0)
+ NER_670_600 = nint(5 * 2.d0*NEX_MAX / 1248.d0)
+ NER_771_670 = nint(6 * 2.d0*NEX_MAX / 1248.d0)
+ NER_TOPDDOUBLEPRIME_771 = nint(114 * 2.d0*NEX_MAX / 1248.d0)
+ NER_CMB_TOPDDOUBLEPRIME = nint(8 * 2.d0*NEX_MAX / 1248.d0)
+ NER_OUTER_CORE = nint(124 * 2.d0*NEX_MAX / 1248.d0)
+ NER_TOP_CENTRAL_CUBE_ICB = nint(13 * 2.d0*NEX_MAX / 1248.d0)
+ R_CENTRAL_CUBE = 985000.d0
+
+!! removed this limit else
+!! removed this limit stop 'problem with this value of NEX_MAX'
+ endif
+
+!----
+!---- change some values in the case of regular PREM with two crustal layers or of 3D models
+!----
+
+! case of regular PREM with two crustal layers: change the time step for small meshes
+! because of a different size of elements in the radial direction in the crust
+ if (HONOR_1D_SPHERICAL_MOHO) then
+ if (.not. ONE_CRUST) then
+ ! case 1D + two crustal layers
+ if (NER_CRUST<2) NER_CRUST=2
+ if(NEX_MAX*multiplication_factor <= 160) then
+ DT = 0.20d0
+ else if(NEX_MAX*multiplication_factor <= 256) then
+ DT = 0.20d0
+ endif
+ endif
+ else
+ ! case 3D
+ if (NER_CRUST<2) NER_CRUST=2
+ if(NEX_MAX*multiplication_factor <= 160) then
+ DT = 0.15d0
+ else if(NEX_MAX*multiplication_factor <= 256) then
+ DT = 0.17d0
+ else if(NEX_MAX*multiplication_factor <= 320) then
+ DT = 0.155d0
+ endif
+ endif
+
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
+ DT = DT*0.20d0
+ endif
+
+
+ if( .not. ATTENUATION_RANGE_PREDEFINED ) then
+ call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
+ MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD)
+ endif
+
+ if(ANGULAR_WIDTH_XI_IN_DEGREES < 90.0d0 .or. &
+ ANGULAR_WIDTH_ETA_IN_DEGREES < 90.0d0 .or. &
+ NEX_MAX > 1248) then
+
+ call auto_ner(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
+ NER_CRUST, NER_80_MOHO, NER_220_80, NER_400_220, NER_600_400, &
+ NER_670_600, NER_771_670, NER_TOPDDOUBLEPRIME_771, &
+ NER_CMB_TOPDDOUBLEPRIME, NER_OUTER_CORE, NER_TOP_CENTRAL_CUBE_ICB, &
+ R_CENTRAL_CUBE, CASE_3D)
+
+ call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
+ MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD)
+
+ call auto_time_stepping(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, DT)
+
+goto 777
+ write(*,*)'##############################################################'
+ write(*,*)
+ write(*,*)' Auto Radial Meshing Code '
+ write(*,*)' Consult read_compute_parameters.f90 and auto_ner.f90 '
+ write(*,*)' This should only be invoked for chunks less than 90 degrees'
+ write(*,*)' and for chunks greater than 1248 elements wide'
+ write(*,*)
+ write(*,*)'CHUNK WIDTH: ', ANGULAR_WIDTH_XI_IN_DEGREES
+ write(*,*)'NEX: ', NEX_MAX
+ write(*,*)'NER_CRUST: ', NER_CRUST
+ write(*,*)'NER_80_MOHO: ', NER_80_MOHO
+ write(*,*)'NER_220_80: ', NER_220_80
+ write(*,*)'NER_400_220: ', NER_400_220
+ write(*,*)'NER_600_400: ', NER_600_400
+ write(*,*)'NER_670_600: ', NER_670_600
+ write(*,*)'NER_771_670: ', NER_771_670
+ write(*,*)'NER_TOPDDOUBLEPRIME_771: ', NER_TOPDDOUBLEPRIME_771
+ write(*,*)'NER_CMB_TOPDDOUBLEPRIME: ', NER_CMB_TOPDDOUBLEPRIME
+ write(*,*)'NER_OUTER_CORE: ', NER_OUTER_CORE
+ write(*,*)'NER_TOP_CENTRAL_CUBE_ICB: ', NER_TOP_CENTRAL_CUBE_ICB
+ write(*,*)'R_CENTRAL_CUBE: ', R_CENTRAL_CUBE
+ write(*,*)'multiplication factor: ', multiplication_factor
+ write(*,*)
+ write(*,*)'DT: ',DT
+ write(*,*)'MIN_ATTENUATION_PERIOD ',MIN_ATTENUATION_PERIOD
+ write(*,*)'MAX_ATTENUATION_PERIOD ',MAX_ATTENUATION_PERIOD
+ write(*,*)
+ write(*,*)'##############################################################'
+777 continue
+
+ if (HONOR_1D_SPHERICAL_MOHO) then
+ if (.not. ONE_CRUST) then
+ ! case 1D + two crustal layers
+ if (NER_CRUST<2) NER_CRUST=2
+ endif
+ else
+ ! case 3D
+ if (NER_CRUST<2) NER_CRUST=2
+ endif
+ endif
+
+
+! take a 5% safety margin on the maximum stable time step
+! which was obtained by trial and error
+ DT = DT * (1.d0 - 0.05d0)
+
+ call read_value_logical(OCEANS, 'model.OCEANS')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(ELLIPTICITY, 'model.ELLIPTICITY')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(TOPOGRAPHY, 'model.TOPOGRAPHY')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(GRAVITY, 'model.GRAVITY')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(ROTATION, 'model.ROTATION')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(ATTENUATION, 'model.ATTENUATION')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ call read_value_logical(ABSORBING_CONDITIONS, 'solver.ABSORBING_CONDITIONS')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ if(ABSORBING_CONDITIONS .and. NCHUNKS == 6) stop 'cannot have absorbing conditions in the full Earth'
+
+ if(ABSORBING_CONDITIONS .and. NCHUNKS == 3) stop 'absorbing conditions not supported for three chunks yet'
+
+ if(ATTENUATION_3D .and. .not. ATTENUATION) stop 'need ATTENUATION to use ATTENUATION_3D'
+
+! radii in PREM or IASP91
+! and normalized density at fluid-solid interface on fluid size for coupling
+! ROCEAN: radius of the ocean (m)
+! RMIDDLE_CRUST: radius of the middle crust (m)
+! RMOHO: radius of the Moho (m)
+! R80: radius of 80 km discontinuity (m)
+! R120: radius of 120 km discontinuity (m) in IASP91
+! R220: radius of 220 km discontinuity (m)
+! R400: radius of 400 km discontinuity (m)
+! R600: radius of 600 km 2nd order discontinuity (m)
+! R670: radius of 670 km discontinuity (m)
+! R771: radius of 771 km 2nd order discontinuity (m)
+! RTOPDDOUBLEPRIME: radius of top of D" 2nd order discontinuity (m)
+! RCMB: radius of CMB (m)
+! RICB: radius of ICB (m)
+
+! by default there is no d120 discontinuity, except in IASP91, therefore set to fictitious value
+ R120 = -1.d0
+
+! value common to all models
+ RHO_OCEANS = 1020.0 / RHOAV
+
+ if(REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91) then
+
+! IASP91
+ ROCEAN = 6371000.d0
+ RMIDDLE_CRUST = 6351000.d0
+ RMOHO = 6336000.d0
+ R80 = 6291000.d0
+ R120 = 6251000.d0
+ R220 = 6161000.d0
+ R400 = 5961000.d0
+! there is no d600 discontinuity in IASP91 therefore this value is useless
+! but it needs to be there for compatibility with other subroutines
+ R600 = R_EARTH - 600000.d0
+ R670 = 5711000.d0
+ R771 = 5611000.d0
+ RTOPDDOUBLEPRIME = 3631000.d0
+ RCMB = 3482000.d0
+ RICB = 1217000.d0
+
+ RHO_TOP_OC = 9900.2379 / RHOAV
+ RHO_BOTTOM_OC = 12168.6383 / RHOAV
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
+
+! our implementation of AK135 has not been checked carefully yet
+! therefore let us doublecheck it carefully one day
+
+! values below corrected by Ying Zhou <yingz at gps.caltech.edu>
+
+! AK135 without the 300 meters of mud layer
+ ROCEAN = 6368000.d0
+ RMIDDLE_CRUST = 6361000.d0
+ RMOHO = 6353000.d0
+ R80 = 6291000.d0
+ R220 = 6161000.d0
+ R400 = 5961000.d0
+ R670 = 5711000.d0
+ RTOPDDOUBLEPRIME = 3631000.d0
+ RCMB = 3479500.d0
+ RICB = 1217500.d0
+
+! values for AK135 that are not discontinuities
+ R600 = 5771000.d0
+ R771 = 5611000.d0
+
+ RHO_TOP_OC = 9914.5000 / RHOAV
+ RHO_BOTTOM_OC = 12139.1000 / RHOAV
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
+
+! values below corrected by Ying Zhou <yingz at gps.caltech.edu>
+
+! 1066A
+ RMOHO = 6360000.d0
+ R400 = 5950000.d0
+ R600 = 5781000.d0
+ R670 = 5700000.d0
+ RCMB = 3484300.d0
+ RICB = 1229480.d0
+
+! values for 1066A that are not discontinuities
+ RTOPDDOUBLEPRIME = 3631000.d0
+ R220 = 6161000.d0
+ R771 = 5611000.d0
+! RMIDDLE_CRUST used only for high resolution FFSW1C model, with 3 elements crust simulations
+! mid_crust = 10 km
+ RMIDDLE_CRUST = 6361000.d0
+ R80 = 6291000.d0
+
+! model 1066A has no oceans, therefore we use the radius of the Earth instead
+ ROCEAN = R_EARTH
+
+ RHO_TOP_OC = 9917.4500 / RHOAV
+ RHO_BOTTOM_OC = 12160.6500 / RHOAV
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
+
+! REF
+ ROCEAN = 6368000.d0
+ RMIDDLE_CRUST = 6356000.d0
+ RMOHO = 6346600.d0
+ R80 = 6291000.d0
+ R220 = 6151000.d0
+ R400 = 5961000.d0
+ R600 = 5771000.d0
+ R670 = 5721000.d0
+ R771 = 5600000.d0
+ RTOPDDOUBLEPRIME = 3630000.d0
+ RCMB = 3479958.d0
+ RICB = 1221491.d0
+
+ RHO_TOP_OC = 9903.48 / RHOAV
+ RHO_BOTTOM_OC = 12166.35 / RHOAV
+
+ else
+
+! PREM
+ ROCEAN = 6368000.d0
+ RMIDDLE_CRUST = 6356000.d0
+ RMOHO = 6346600.d0
+ R80 = 6291000.d0
+ R220 = 6151000.d0
+ R400 = 5971000.d0
+ R600 = 5771000.d0
+ R670 = 5701000.d0
+ R771 = 5600000.d0
+ RTOPDDOUBLEPRIME = 3630000.d0
+ RCMB = 3480000.d0
+ RICB = 1221000.d0
+
+ RHO_TOP_OC = 9903.4384 / RHOAV
+ RHO_BOTTOM_OC = 12166.5885 / RHOAV
+
+ endif
+
+! honor the PREM Moho or define a fictitious Moho in order to have even radial sampling
+! from the d220 to the Earth surface
+ if(HONOR_1D_SPHERICAL_MOHO) then
+ RMOHO_FICTITIOUS_IN_MESHER = RMOHO
+ else
+ RMOHO_FICTITIOUS_IN_MESHER = (R80 + R_EARTH) / 2
+ endif
+
+ call read_value_double_precision(RECORD_LENGTH_IN_MINUTES, 'solver.RECORD_LENGTH_IN_MINUTES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+! compute total number of time steps, rounded to next multiple of 100
+ NSTEP = 100 * (int(RECORD_LENGTH_IN_MINUTES * 60.d0 / (100.d0*DT)) + 1)
+
+ call read_value_logical(MOVIE_SURFACE, 'solver.MOVIE_SURFACE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(MOVIE_VOLUME, 'solver.MOVIE_VOLUME')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NTSTEP_BETWEEN_FRAMES, 'solver.NTSTEP_BETWEEN_FRAMES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(HDUR_MOVIE, 'solver.HDUR_MOVIE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+! computes a default hdur_movie that creates nice looking movies.
+! Sets HDUR_MOVIE as the minimum period the mesh can resolve
+ if(HDUR_MOVIE <= TINYVAL) &
+ HDUR_MOVIE = 1.1d0*max(240.d0/NEX_XI*18.d0*ANGULAR_WIDTH_XI_IN_DEGREES/90.d0, &
+ 240.d0/NEX_ETA*18.d0*ANGULAR_WIDTH_ETA_IN_DEGREES/90.d0)
+
+ call read_value_integer(MOVIE_VOLUME_TYPE, 'solver.MOVIE_VOLUME_TYPE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(MOVIE_VOLUME_COARSE,'solver.MOVIE_VOLUME_COARSE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_TOP_KM, 'solver.MOVIE_TOP_KM')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_BOTTOM_KM, 'solver.MOVIE_BOTTOM_KM')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_WEST_DEG, 'solver.MOVIE_WEST_DEG')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_EAST_DEG, 'solver.MOVIE_EAST_DEG')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_NORTH_DEG, 'solver.MOVIE_NORTH_DEG')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_double_precision(MOVIE_SOUTH_DEG, 'solver.MOVIE_SOUTH_DEG')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(MOVIE_START, 'solver.MOVIE_START')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(MOVIE_STOP, 'solver.MOVIE_STOP')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ MOVIE_TOP = (R_EARTH_KM-MOVIE_TOP_KM)/R_EARTH_KM
+ MOVIE_BOTTOM = (R_EARTH_KM-MOVIE_BOTTOM_KM)/R_EARTH_KM
+ MOVIE_EAST = MOVIE_EAST_DEG * DEGREES_TO_RADIANS
+ MOVIE_WEST = MOVIE_WEST_DEG * DEGREES_TO_RADIANS
+ MOVIE_NORTH = (90.0d0 - MOVIE_NORTH_DEG) * DEGREES_TO_RADIANS ! converting from latitude to colatitude
+ MOVIE_SOUTH = (90.0d0 - MOVIE_SOUTH_DEG) * DEGREES_TO_RADIANS
+
+ call read_value_logical(SAVE_MESH_FILES, 'mesher.SAVE_MESH_FILES')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NUMBER_OF_RUNS, 'solver.NUMBER_OF_RUNS')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NUMBER_OF_THIS_RUN, 'solver.NUMBER_OF_THIS_RUN')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_string(LOCAL_PATH, 'LOCAL_PATH')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NTSTEP_BETWEEN_OUTPUT_INFO, 'solver.NTSTEP_BETWEEN_OUTPUT_INFO')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NTSTEP_BETWEEN_OUTPUT_SEISMOS, 'solver.NTSTEP_BETWEEN_OUTPUT_SEISMOS')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_integer(NTSTEP_BETWEEN_READ_ADJSRC, 'solver.NTSTEP_BETWEEN_READ_ADJSRC')
+ if(err_occurred() /= 0) return
+
+ call read_value_logical(OUTPUT_SEISMOS_ASCII_TEXT, 'solver.OUTPUT_SEISMOS_ASCII_TEXT')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(OUTPUT_SEISMOS_SAC_ALPHANUM, 'solver.OUTPUT_SEISMOS_SAC_ALPHANUM')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(OUTPUT_SEISMOS_SAC_BINARY, 'solver.OUTPUT_SEISMOS_SAC_BINARY')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(ROTATE_SEISMOGRAMS_RT, 'solver.ROTATE_SEISMOGRAMS_RT')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(WRITE_SEISMOGRAMS_BY_MASTER, 'solver.WRITE_SEISMOGRAMS_BY_MASTER')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(SAVE_ALL_SEISMOS_IN_ONE_FILE, 'solver.SAVE_ALL_SEISMOS_IN_ONE_FILE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(USE_BINARY_FOR_LARGE_FILE, 'solver.USE_BINARY_FOR_LARGE_FILE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ call read_value_logical(RECEIVERS_CAN_BE_BURIED, 'solver.RECEIVERS_CAN_BE_BURIED')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ call read_value_logical(PRINT_SOURCE_TIME_FUNCTION, 'solver.PRINT_SOURCE_TIME_FUNCTION')
+
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+! close parameter file
+ call close_parameter_file
+!--- check that parameters make sense
+
+ if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
+ stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5, check the Par_file'
+
+! subsets used to save seismograms must not be larger than the whole time series,
+! otherwise we waste memory
+ if(NTSTEP_BETWEEN_OUTPUT_SEISMOS > NSTEP) then
+ NTSTEP_BETWEEN_OUTPUT_SEISMOS = NSTEP
+ if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
+ stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then modified NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5'
+ endif
+
+! check that reals are either 4 or 8 bytes
+ if(CUSTOM_REAL /= SIZE_REAL .and. CUSTOM_REAL /= SIZE_DOUBLE) stop 'wrong size of CUSTOM_REAL for reals'
+
+! check that the parameter file is correct
+ if(NGNOD /= 27) stop 'number of control nodes must be 27'
+ if(NGNOD == 27 .and. NGNOD2D /= 9) stop 'elements with 27 points should have NGNOD2D = 9'
+
+! for the number of standard linear solids for attenuation
+ if(N_SLS /= 3) stop 'number of SLS must be 3'
+
+! check number of slices in each direction
+ if(NCHUNKS < 1) stop 'must have at least one chunk'
+ if(NPROC_XI < 1) stop 'NPROC_XI must be at least 1'
+ if(NPROC_ETA < 1) stop 'NPROC_ETA must be at least 1'
+
+! check number of chunks
+ if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
+ stop 'only one, two, three or six chunks can be meshed'
+
+! check that the central cube can be included
+ if(INCLUDE_CENTRAL_CUBE .and. NCHUNKS /= 6) stop 'need six chunks to include central cube'
+
+! check that sphere can be cut into slices without getting negative Jacobian
+ if(NEX_XI < 48) stop 'NEX_XI must be greater than 48 to cut the sphere into slices with positive Jacobian'
+ if(NEX_ETA < 48) stop 'NEX_ETA must be greater than 48 to cut the sphere into slices with positive Jacobian'
+
+! check that mesh can be coarsened in depth three or four times
+ CUT_SUPERBRICK_XI=.false.
+ CUT_SUPERBRICK_ETA=.false.
+
+ if (SUPPRESS_CRUSTAL_MESH .and. .not. ADD_4TH_DOUBLING) then
+ if(mod(NEX_XI,8) /= 0) stop 'NEX_XI must be a multiple of 8'
+ if(mod(NEX_ETA,8) /= 0) stop 'NEX_ETA must be a multiple of 8'
+ if(mod(NEX_XI/4,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 4*NPROC_XI'
+ if(mod(NEX_ETA/4,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 4*NPROC_ETA'
+ if(mod(NEX_XI/8,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
+ if(mod(NEX_ETA/8,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
+ elseif (SUPPRESS_CRUSTAL_MESH .or. .not. ADD_4TH_DOUBLING) then
+ if(mod(NEX_XI,16) /= 0) stop 'NEX_XI must be a multiple of 16'
+ if(mod(NEX_ETA,16) /= 0) stop 'NEX_ETA must be a multiple of 16'
+ if(mod(NEX_XI/8,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 8*NPROC_XI'
+ if(mod(NEX_ETA/8,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 8*NPROC_ETA'
+ if(mod(NEX_XI/16,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
+ if(mod(NEX_ETA/16,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
+ else
+ if(mod(NEX_XI,32) /= 0) stop 'NEX_XI must be a multiple of 32'
+ if(mod(NEX_ETA,32) /= 0) stop 'NEX_ETA must be a multiple of 32'
+ if(mod(NEX_XI/16,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 16*NPROC_XI'
+ if(mod(NEX_ETA/16,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 16*NPROC_ETA'
+ if(mod(NEX_XI/32,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
+ if(mod(NEX_ETA/32,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
+ endif
+
+! check that topology is correct if more than two chunks
+ if(NCHUNKS > 2 .and. NEX_XI /= NEX_ETA) stop 'must have NEX_XI = NEX_ETA for more than two chunks'
+ if(NCHUNKS > 2 .and. NPROC_XI /= NPROC_ETA) stop 'must have NPROC_XI = NPROC_ETA for more than two chunks'
+
+! check that IASP91, AK135, or 1066A is isotropic
+ if((REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91 .or. &
+ REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135 .or. &
+ REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) .and. TRANSVERSE_ISOTROPY) &
+ stop 'models IASP91, AK135 and 1066A are currently isotropic'
+
+ ELEMENT_WIDTH = ANGULAR_WIDTH_XI_IN_DEGREES/dble(NEX_MAX) * DEGREES_TO_RADIANS
+
+!
+!--- compute additional parameters
+!
+
+! 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
+ NPROC = NPROC_XI * NPROC_ETA
+
+! total number of processors in the full Earth composed of the six chunks
+ NPROCTOT = NCHUNKS * NPROC
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! definition of general mesh parameters below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! find element below top of which we should implement the second doubling in the mantle
+! locate element closest to optimal value
+ distance_min = HUGEVAL
+ do ielem = 2,NER_TOPDDOUBLEPRIME_771
+ zval = RTOPDDOUBLEPRIME + ielem * (R771 - RTOPDDOUBLEPRIME) / dble(NER_TOPDDOUBLEPRIME_771)
+ distance = abs(zval - (R_EARTH - DEPTH_SECOND_DOUBLING_OPTIMAL))
+ if(distance < distance_min) then
+ elem_doubling_mantle = ielem
+ distance_min = distance
+ DEPTH_SECOND_DOUBLING_REAL = R_EARTH - zval
+ endif
+ enddo
+
+! find element below top of which we should implement the third doubling in the middle of the outer core
+! locate element closest to optimal value
+ distance_min = HUGEVAL
+! start at element number 4 because we need at least two elements below for the fourth doubling
+! implemented at the bottom of the outer core
+ do ielem = 4,NER_OUTER_CORE
+ zval = RICB + ielem * (RCMB - RICB) / dble(NER_OUTER_CORE)
+ distance = abs(zval - (R_EARTH - DEPTH_THIRD_DOUBLING_OPTIMAL))
+ if(distance < distance_min) then
+ elem_doubling_middle_outer_core = ielem
+ distance_min = distance
+ DEPTH_THIRD_DOUBLING_REAL = R_EARTH - zval
+ endif
+ enddo
+
+ if (ADD_4TH_DOUBLING) then
+! find element below top of which we should implement the fourth doubling in the middle of the outer core
+! locate element closest to optimal value
+ distance_min = HUGEVAL
+! end two elements before the top because we need at least two elements above for the third doubling
+! implemented in the middle of the outer core
+ do ielem = 2,NER_OUTER_CORE-2
+ zval = RICB + ielem * (RCMB - RICB) / dble(NER_OUTER_CORE)
+ distance = abs(zval - (R_EARTH - DEPTH_FOURTH_DOUBLING_OPTIMAL))
+ if(distance < distance_min) then
+ elem_doubling_bottom_outer_core = ielem
+ distance_min = distance
+ DEPTH_FOURTH_DOUBLING_REAL = R_EARTH - zval
+ endif
+ enddo
+! make sure that the two doublings in the outer core are found in the right order
+ if(elem_doubling_bottom_outer_core >= elem_doubling_middle_outer_core) &
+ stop 'error in location of the two doublings in the outer core'
+ endif
+
+ ratio_sampling_array(15) = 0
+
+! define all the layers of the mesh
+ if (.not. ADD_4TH_DOUBLING) then
+
+ if (SUPPRESS_CRUSTAL_MESH) then
+
+ ONE_CRUST = .false.
+ OCEANS= .false.
+ TOPOGRAPHY = .false.
+ CRUSTAL = .false.
+
+ NUMBER_OF_MESH_LAYERS = 14
+ layer_offset = 1
+
+ ! now only one region
+ ner( 1) = NER_CRUST + NER_80_MOHO
+ ner( 2) = 0
+ ner( 3) = 0
+
+ ner( 4) = NER_220_80
+ ner( 5) = NER_400_220
+ ner( 6) = NER_600_400
+ ner( 7) = NER_670_600
+ ner( 8) = NER_771_670
+ ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner(10) = elem_doubling_mantle
+ ner(11) = NER_CMB_TOPDDOUBLEPRIME
+ ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(13) = elem_doubling_middle_outer_core
+ ner(14) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1:9) = 1
+ ratio_sampling_array(10:12) = 2
+ ratio_sampling_array(13:14) = 4
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1:3) = IFLAG_CRUST !!!!! IFLAG_80_MOHO
+ doubling_index(4) = IFLAG_220_80
+ doubling_index(5:7) = IFLAG_670_220
+ doubling_index(8:11) = IFLAG_MANTLE_NORMAL
+ doubling_index(12:13) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(14) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(10) = .true.
+ this_region_has_a_doubling(13) = .true.
+ lastdoubling_layer = 13
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = R80
+
+ r_top(2) = RMIDDLE_CRUST !!!! now fictitious
+ r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
+
+ r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
+ r_bottom(3) = R80 !!!! now fictitious
+
+ r_top(4) = R80
+ r_bottom(4) = R220
+
+ r_top(5) = R220
+ r_bottom(5) = R400
+
+ r_top(6) = R400
+ r_bottom(6) = R600
+
+ r_top(7) = R600
+ r_bottom(7) = R670
+
+ r_top(8) = R670
+ r_bottom(8) = R771
+
+ r_top(9) = R771
+ r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(10) = RTOPDDOUBLEPRIME
+
+ r_top(11) = RTOPDDOUBLEPRIME
+ r_bottom(11) = RCMB
+
+ r_top(12) = RCMB
+ r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(13) = RICB
+
+ r_top(15) = RICB
+ r_bottom(15) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = R80 / R_EARTH
+
+ rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
+ rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+
+ rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+ rmins(3) = R80 / R_EARTH !!!! now fictitious
+
+ rmaxs(4) = R80 / R_EARTH
+ rmins(4) = R220 / R_EARTH
+
+ rmaxs(5) = R220 / R_EARTH
+ rmins(5) = R400 / R_EARTH
+
+ rmaxs(6) = R400 / R_EARTH
+ rmins(6) = R600 / R_EARTH
+
+ rmaxs(7) = R600 / R_EARTH
+ rmins(7) = R670 / R_EARTH
+
+ rmaxs(8) = R670 / R_EARTH
+ rmins(8) = R771 / R_EARTH
+
+ rmaxs(9:10) = R771 / R_EARTH
+ rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(11) = RCMB / R_EARTH
+
+ rmaxs(12:13) = RCMB / R_EARTH
+ rmins(12:13) = RICB / R_EARTH
+
+ rmaxs(14) = RICB / R_EARTH
+ rmins(14) = R_CENTRAL_CUBE / R_EARTH
+
+ elseif (ONE_CRUST) then
+
+ NUMBER_OF_MESH_LAYERS = 13
+ layer_offset = 0
+
+ ner( 1) = NER_CRUST
+ ner( 2) = NER_80_MOHO
+ ner( 3) = NER_220_80
+ ner( 4) = NER_400_220
+ ner( 5) = NER_600_400
+ ner( 6) = NER_670_600
+ ner( 7) = NER_771_670
+ ner( 8) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner( 9) = elem_doubling_mantle
+ ner(10) = NER_CMB_TOPDDOUBLEPRIME
+ ner(11) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(12) = elem_doubling_middle_outer_core
+ ner(13) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1) = 1
+ ratio_sampling_array(2:8) = 2
+ ratio_sampling_array(9:11) = 4
+ ratio_sampling_array(12:13) = 8
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1) = IFLAG_CRUST
+ doubling_index(2) = IFLAG_80_MOHO
+ doubling_index(3) = IFLAG_220_80
+ doubling_index(4:6) = IFLAG_670_220
+ doubling_index(7:10) = IFLAG_MANTLE_NORMAL
+ doubling_index(11:12) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(13) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(2) = .true.
+ this_region_has_a_doubling(9) = .true.
+ this_region_has_a_doubling(12) = .true.
+ lastdoubling_layer = 12
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ !!!!!!!!!!! DK DK UGLY: beware, is there a bug when 3D crust crosses anisotropy in the mantle?
+ !!!!!!!!!!! DK DK UGLY: i.e. if there is no thick crust there, some elements above the Moho
+ !!!!!!!!!!! DK DK UGLY: should be anisotropic but anisotropy is currently only
+ !!!!!!!!!!! DK DK UGLY: stored between d220 and MOHO to save memory? Clarify this one day.
+ !!!!!!!!!!! DK DK UGLY: The Moho stretching and squishing that Jeroen added to V4.0
+ !!!!!!!!!!! DK DK UGLY: should partly deal with this problem.
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
+
+ r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
+ r_bottom(2) = R80
+
+ r_top(3) = R80
+ r_bottom(3) = R220
+
+ r_top(4) = R220
+ r_bottom(4) = R400
+
+ r_top(5) = R400
+ r_bottom(5) = R600
+
+ r_top(6) = R600
+ r_bottom(6) = R670
+
+ r_top(7) = R670
+ r_bottom(7) = R771
+
+ r_top(8) = R771
+ r_bottom(8) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(9) = RTOPDDOUBLEPRIME
+
+ r_top(10) = RTOPDDOUBLEPRIME
+ r_bottom(10) = RCMB
+
+ r_top(11) = RCMB
+ r_bottom(11) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(12) = RICB
+
+ r_top(13) = RICB
+ r_bottom(13) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+ rmins(2) = R80 / R_EARTH
+
+ rmaxs(3) = R80 / R_EARTH
+ rmins(3) = R220 / R_EARTH
+
+ rmaxs(4) = R220 / R_EARTH
+ rmins(4) = R400 / R_EARTH
+
+ rmaxs(5) = R400 / R_EARTH
+ rmins(5) = R600 / R_EARTH
+
+ rmaxs(6) = R600 / R_EARTH
+ rmins(6) = R670 / R_EARTH
+
+ rmaxs(7) = R670 / R_EARTH
+ rmins(7) = R771 / R_EARTH
+
+ rmaxs(8:9) = R771 / R_EARTH
+ rmins(8:9) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(10) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(10) = RCMB / R_EARTH
+
+ rmaxs(11:12) = RCMB / R_EARTH
+ rmins(11:12) = RICB / R_EARTH
+
+ rmaxs(13) = RICB / R_EARTH
+ rmins(13) = R_CENTRAL_CUBE / R_EARTH
+
+ else
+
+ NUMBER_OF_MESH_LAYERS = 14
+ layer_offset = 1
+ if ((RMIDDLE_CRUST-RMOHO_FICTITIOUS_IN_MESHER)<(R_EARTH-RMIDDLE_CRUST)) then
+ ner( 1) = ceiling (NER_CRUST / 2.d0)
+ ner( 2) = floor (NER_CRUST / 2.d0)
+ else
+ ner( 1) = floor (NER_CRUST / 2.d0)
+ ner( 2) = ceiling (NER_CRUST / 2.d0)
+ endif
+ ner( 3) = NER_80_MOHO
+ ner( 4) = NER_220_80
+ ner( 5) = NER_400_220
+ ner( 6) = NER_600_400
+ ner( 7) = NER_670_600
+ ner( 8) = NER_771_670
+ ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner(10) = elem_doubling_mantle
+ ner(11) = NER_CMB_TOPDDOUBLEPRIME
+ ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(13) = elem_doubling_middle_outer_core
+ ner(14) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1:2) = 1
+ ratio_sampling_array(3:9) = 2
+ ratio_sampling_array(10:12) = 4
+ ratio_sampling_array(13:14) = 8
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1:2) = IFLAG_CRUST
+ doubling_index(3) = IFLAG_80_MOHO
+ doubling_index(4) = IFLAG_220_80
+ doubling_index(5:7) = IFLAG_670_220
+ doubling_index(8:11) = IFLAG_MANTLE_NORMAL
+ doubling_index(12:13) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(14) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(3) = .true.
+ this_region_has_a_doubling(10) = .true.
+ this_region_has_a_doubling(13) = .true.
+ this_region_has_a_doubling(14) = .false.
+ lastdoubling_layer = 13
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = RMIDDLE_CRUST
+
+ r_top(2) = RMIDDLE_CRUST
+ r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
+
+ r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
+ r_bottom(3) = R80
+
+ r_top(4) = R80
+ r_bottom(4) = R220
+
+ r_top(5) = R220
+ r_bottom(5) = R400
+
+ r_top(6) = R400
+ r_bottom(6) = R600
+
+ r_top(7) = R600
+ r_bottom(7) = R670
+
+ r_top(8) = R670
+ r_bottom(8) = R771
+
+ r_top(9) = R771
+ r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(10) = RTOPDDOUBLEPRIME
+
+ r_top(11) = RTOPDDOUBLEPRIME
+ r_bottom(11) = RCMB
+
+ r_top(12) = RCMB
+ r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(13) = RICB
+
+ r_top(14) = RICB
+ r_bottom(14) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = RMIDDLE_CRUST / R_EARTH
+
+ rmaxs(2) = RMIDDLE_CRUST / R_EARTH
+ rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+ rmins(3) = R80 / R_EARTH
+
+ rmaxs(4) = R80 / R_EARTH
+ rmins(4) = R220 / R_EARTH
+
+ rmaxs(5) = R220 / R_EARTH
+ rmins(5) = R400 / R_EARTH
+
+ rmaxs(6) = R400 / R_EARTH
+ rmins(6) = R600 / R_EARTH
+
+ rmaxs(7) = R600 / R_EARTH
+ rmins(7) = R670 / R_EARTH
+
+ rmaxs(8) = R670 / R_EARTH
+ rmins(8) = R771 / R_EARTH
+
+ rmaxs(9:10) = R771 / R_EARTH
+ rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(11) = RCMB / R_EARTH
+
+ rmaxs(12:13) = RCMB / R_EARTH
+ rmins(12:13) = RICB / R_EARTH
+
+ rmaxs(14) = RICB / R_EARTH
+ rmins(14) = R_CENTRAL_CUBE / R_EARTH
+
+ endif
+ else
+ if (SUPPRESS_CRUSTAL_MESH) then
+
+ ONE_CRUST = .false.
+ OCEANS= .false.
+ TOPOGRAPHY = .false.
+ CRUSTAL = .false.
+
+ NUMBER_OF_MESH_LAYERS = 15
+ layer_offset = 1
+
+ ! now only one region
+ ner( 1) = NER_CRUST + NER_80_MOHO
+ ner( 2) = 0
+ ner( 3) = 0
+
+ ner( 4) = NER_220_80
+ ner( 5) = NER_400_220
+ ner( 6) = NER_600_400
+ ner( 7) = NER_670_600
+ ner( 8) = NER_771_670
+ ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner(10) = elem_doubling_mantle
+ ner(11) = NER_CMB_TOPDDOUBLEPRIME
+ ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(13) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
+ ner(14) = elem_doubling_bottom_outer_core
+ ner(15) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1:9) = 1
+ ratio_sampling_array(10:12) = 2
+ ratio_sampling_array(13) = 4
+ ratio_sampling_array(14:15) = 8
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1:3) = IFLAG_CRUST !!!!! IFLAG_80_MOHO
+ doubling_index(4) = IFLAG_220_80
+ doubling_index(5:7) = IFLAG_670_220
+ doubling_index(8:11) = IFLAG_MANTLE_NORMAL
+ doubling_index(12:14) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(15) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(10) = .true.
+ this_region_has_a_doubling(13) = .true.
+ this_region_has_a_doubling(14) = .true.
+ lastdoubling_layer = 14
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = R80
+
+ r_top(2) = RMIDDLE_CRUST !!!! now fictitious
+ r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
+
+ r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
+ r_bottom(3) = R80 !!!! now fictitious
+
+ r_top(4) = R80
+ r_bottom(4) = R220
+
+ r_top(5) = R220
+ r_bottom(5) = R400
+
+ r_top(6) = R400
+ r_bottom(6) = R600
+
+ r_top(7) = R600
+ r_bottom(7) = R670
+
+ r_top(8) = R670
+ r_bottom(8) = R771
+
+ r_top(9) = R771
+ r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(10) = RTOPDDOUBLEPRIME
+
+ r_top(11) = RTOPDDOUBLEPRIME
+ r_bottom(11) = RCMB
+
+ r_top(12) = RCMB
+ r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+
+ r_top(14) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+ r_bottom(14) = RICB
+
+ r_top(15) = RICB
+ r_bottom(15) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = R80 / R_EARTH
+
+ rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
+ rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+
+ rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+ rmins(3) = R80 / R_EARTH !!!! now fictitious
+
+ rmaxs(4) = R80 / R_EARTH
+ rmins(4) = R220 / R_EARTH
+
+ rmaxs(5) = R220 / R_EARTH
+ rmins(5) = R400 / R_EARTH
+
+ rmaxs(6) = R400 / R_EARTH
+ rmins(6) = R600 / R_EARTH
+
+ rmaxs(7) = R600 / R_EARTH
+ rmins(7) = R670 / R_EARTH
+
+ rmaxs(8) = R670 / R_EARTH
+ rmins(8) = R771 / R_EARTH
+
+ rmaxs(9:10) = R771 / R_EARTH
+ rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(11) = RCMB / R_EARTH
+
+ rmaxs(12:14) = RCMB / R_EARTH
+ rmins(12:14) = RICB / R_EARTH
+
+ rmaxs(15) = RICB / R_EARTH
+ rmins(15) = R_CENTRAL_CUBE / R_EARTH
+
+ elseif (ONE_CRUST) then
+
+ NUMBER_OF_MESH_LAYERS = 14
+ layer_offset = 0
+
+ ner( 1) = NER_CRUST
+ ner( 2) = NER_80_MOHO
+ ner( 3) = NER_220_80
+ ner( 4) = NER_400_220
+ ner( 5) = NER_600_400
+ ner( 6) = NER_670_600
+ ner( 7) = NER_771_670
+ ner( 8) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner( 9) = elem_doubling_mantle
+ ner(10) = NER_CMB_TOPDDOUBLEPRIME
+ ner(11) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(12) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
+ ner(13) = elem_doubling_bottom_outer_core
+ ner(14) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1) = 1
+ ratio_sampling_array(2:8) = 2
+ ratio_sampling_array(9:11) = 4
+ ratio_sampling_array(12) = 8
+ ratio_sampling_array(13:14) = 16
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1) = IFLAG_CRUST
+ doubling_index(2) = IFLAG_80_MOHO
+ doubling_index(3) = IFLAG_220_80
+ doubling_index(4:6) = IFLAG_670_220
+ doubling_index(7:10) = IFLAG_MANTLE_NORMAL
+ doubling_index(11:13) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(14) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(2) = .true.
+ this_region_has_a_doubling(9) = .true.
+ this_region_has_a_doubling(12) = .true.
+ this_region_has_a_doubling(13) = .true.
+ lastdoubling_layer = 13
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ !!!!!!!!!!! DK DK UGLY: beware, is there a bug when 3D crust crosses anisotropy in the mantle?
+ !!!!!!!!!!! DK DK UGLY: i.e. if there is no thick crust there, some elements above the Moho
+ !!!!!!!!!!! DK DK UGLY: should be anisotropic but anisotropy is currently only
+ !!!!!!!!!!! DK DK UGLY: stored between d220 and MOHO to save memory? Clarify this one day.
+ !!!!!!!!!!! DK DK UGLY: The Moho stretching and squishing that Jeroen added to V4.0
+ !!!!!!!!!!! DK DK UGLY: should partly deal with this problem.
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
+
+ r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
+ r_bottom(2) = R80
+
+ r_top(3) = R80
+ r_bottom(3) = R220
+
+ r_top(4) = R220
+ r_bottom(4) = R400
+
+ r_top(5) = R400
+ r_bottom(5) = R600
+
+ r_top(6) = R600
+ r_bottom(6) = R670
+
+ r_top(7) = R670
+ r_bottom(7) = R771
+
+ r_top(8) = R771
+ r_bottom(8) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(9) = RTOPDDOUBLEPRIME
+
+ r_top(10) = RTOPDDOUBLEPRIME
+ r_bottom(10) = RCMB
+
+ r_top(11) = RCMB
+ r_bottom(11) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(12) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+
+ r_top(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+ r_bottom(13) = RICB
+
+ r_top(14) = RICB
+ r_bottom(14) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+ rmins(2) = R80 / R_EARTH
+
+ rmaxs(3) = R80 / R_EARTH
+ rmins(3) = R220 / R_EARTH
+
+ rmaxs(4) = R220 / R_EARTH
+ rmins(4) = R400 / R_EARTH
+
+ rmaxs(5) = R400 / R_EARTH
+ rmins(5) = R600 / R_EARTH
+
+ rmaxs(6) = R600 / R_EARTH
+ rmins(6) = R670 / R_EARTH
+
+ rmaxs(7) = R670 / R_EARTH
+ rmins(7) = R771 / R_EARTH
+
+ rmaxs(8:9) = R771 / R_EARTH
+ rmins(8:9) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(10) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(10) = RCMB / R_EARTH
+
+ rmaxs(11:13) = RCMB / R_EARTH
+ rmins(11:13) = RICB / R_EARTH
+
+ rmaxs(14) = RICB / R_EARTH
+ rmins(14) = R_CENTRAL_CUBE / R_EARTH
+
+ else
+
+ NUMBER_OF_MESH_LAYERS = 15
+ layer_offset = 1
+ if ((RMIDDLE_CRUST-RMOHO_FICTITIOUS_IN_MESHER)<(R_EARTH-RMIDDLE_CRUST)) then
+ ner( 1) = ceiling (NER_CRUST / 2.d0)
+ ner( 2) = floor (NER_CRUST / 2.d0)
+ else
+ ner( 1) = floor (NER_CRUST / 2.d0)
+ ner( 2) = ceiling (NER_CRUST / 2.d0)
+ endif
+ ner( 3) = NER_80_MOHO
+ ner( 4) = NER_220_80
+ ner( 5) = NER_400_220
+ ner( 6) = NER_600_400
+ ner( 7) = NER_670_600
+ ner( 8) = NER_771_670
+ ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
+ ner(10) = elem_doubling_mantle
+ ner(11) = NER_CMB_TOPDDOUBLEPRIME
+ ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
+ ner(13) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
+ ner(14) = elem_doubling_bottom_outer_core
+ ner(15) = NER_TOP_CENTRAL_CUBE_ICB
+
+ ! value of the doubling ratio in each radial region of the mesh
+ ratio_sampling_array(1:2) = 1
+ ratio_sampling_array(3:9) = 2
+ ratio_sampling_array(10:12) = 4
+ ratio_sampling_array(13) = 8
+ ratio_sampling_array(14:15) = 16
+
+ ! value of the doubling index flag in each radial region of the mesh
+ doubling_index(1:2) = IFLAG_CRUST
+ doubling_index(3) = IFLAG_80_MOHO
+ doubling_index(4) = IFLAG_220_80
+ doubling_index(5:7) = IFLAG_670_220
+ doubling_index(8:11) = IFLAG_MANTLE_NORMAL
+ doubling_index(12:14) = IFLAG_OUTER_CORE_NORMAL
+ doubling_index(15) = IFLAG_INNER_CORE_NORMAL
+
+ ! define the three regions in which we implement a mesh doubling at the top of that region
+ this_region_has_a_doubling(:) = .false.
+ this_region_has_a_doubling(3) = .true.
+ this_region_has_a_doubling(10) = .true.
+ this_region_has_a_doubling(13) = .true.
+ this_region_has_a_doubling(14) = .true.
+ lastdoubling_layer = 14
+
+ ! define the top and bottom radii of all the regions of the mesh in the radial direction
+ ! the first region is the crust at the surface of the Earth
+ ! the last region is in the inner core near the center of the Earth
+
+ r_top(1) = R_EARTH
+ r_bottom(1) = RMIDDLE_CRUST
+
+ r_top(2) = RMIDDLE_CRUST
+ r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
+
+ r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
+ r_bottom(3) = R80
+
+ r_top(4) = R80
+ r_bottom(4) = R220
+
+ r_top(5) = R220
+ r_bottom(5) = R400
+
+ r_top(6) = R400
+ r_bottom(6) = R600
+
+ r_top(7) = R600
+ r_bottom(7) = R670
+
+ r_top(8) = R670
+ r_bottom(8) = R771
+
+ r_top(9) = R771
+ r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+
+ r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
+ r_bottom(10) = RTOPDDOUBLEPRIME
+
+ r_top(11) = RTOPDDOUBLEPRIME
+ r_bottom(11) = RCMB
+
+ r_top(12) = RCMB
+ r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+
+ r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
+ r_bottom(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+
+ r_top(14) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
+ r_bottom(14) = RICB
+
+ r_top(15) = RICB
+ r_bottom(15) = R_CENTRAL_CUBE
+
+ ! new definition of rmins & rmaxs
+ rmaxs(1) = ONE
+ rmins(1) = RMIDDLE_CRUST / R_EARTH
+
+ rmaxs(2) = RMIDDLE_CRUST / R_EARTH
+ rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+ rmins(3) = R80 / R_EARTH
+
+ rmaxs(4) = R80 / R_EARTH
+ rmins(4) = R220 / R_EARTH
+
+ rmaxs(5) = R220 / R_EARTH
+ rmins(5) = R400 / R_EARTH
+
+ rmaxs(6) = R400 / R_EARTH
+ rmins(6) = R600 / R_EARTH
+
+ rmaxs(7) = R600 / R_EARTH
+ rmins(7) = R670 / R_EARTH
+
+ rmaxs(8) = R670 / R_EARTH
+ rmins(8) = R771 / R_EARTH
+
+ rmaxs(9:10) = R771 / R_EARTH
+ rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
+
+ rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
+ rmins(11) = RCMB / R_EARTH
+
+ rmaxs(12:14) = RCMB / R_EARTH
+ rmins(12:14) = RICB / R_EARTH
+
+ rmaxs(15) = RICB / R_EARTH
+ rmins(15) = R_CENTRAL_CUBE / R_EARTH
+ endif
+ endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! calculation of number of elements (NSPEC) below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ ratio_divide_central_cube = maxval(ratio_sampling_array)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 1D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! theoretical number of spectral elements in radial direction
+do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
+ if(iter_region == IREGION_CRUST_MANTLE) then
+ ifirst_region = 1
+ ilast_region = 10 + layer_offset
+ else if(iter_region == IREGION_OUTER_CORE) then
+ ifirst_region = 11 + layer_offset
+ ilast_region = NUMBER_OF_MESH_LAYERS - 1
+ else if(iter_region == IREGION_INNER_CORE) then
+ ifirst_region = NUMBER_OF_MESH_LAYERS
+ ilast_region = NUMBER_OF_MESH_LAYERS
+ else
+ stop 'incorrect region code detected'
+ endif
+ NSPEC1D_RADIAL(iter_region) = sum(ner(ifirst_region:ilast_region))
+enddo
+
+! difference of radial number of element for outer core if the superbrick is cut
+ DIFF_NSPEC1D_RADIAL(:,:) = 0
+ if (CUT_SUPERBRICK_XI) then
+ if (CUT_SUPERBRICK_ETA) then
+ DIFF_NSPEC1D_RADIAL(2,1) = 1
+ DIFF_NSPEC1D_RADIAL(3,1) = 2
+ DIFF_NSPEC1D_RADIAL(4,1) = 1
+
+ DIFF_NSPEC1D_RADIAL(1,2) = 1
+ DIFF_NSPEC1D_RADIAL(2,2) = 2
+ DIFF_NSPEC1D_RADIAL(3,2) = 1
+
+ DIFF_NSPEC1D_RADIAL(1,3) = 1
+ DIFF_NSPEC1D_RADIAL(3,3) = 1
+ DIFF_NSPEC1D_RADIAL(4,3) = 2
+
+ DIFF_NSPEC1D_RADIAL(1,4) = 2
+ DIFF_NSPEC1D_RADIAL(2,4) = 1
+ DIFF_NSPEC1D_RADIAL(4,4) = 1
+ else
+ DIFF_NSPEC1D_RADIAL(2,1) = 1
+ DIFF_NSPEC1D_RADIAL(3,1) = 1
+
+ DIFF_NSPEC1D_RADIAL(1,2) = 1
+ DIFF_NSPEC1D_RADIAL(4,2) = 1
+ endif
+ else
+ if (CUT_SUPERBRICK_ETA) then
+ DIFF_NSPEC1D_RADIAL(3,1) = 1
+ DIFF_NSPEC1D_RADIAL(4,1) = 1
+
+ DIFF_NSPEC1D_RADIAL(1,2) = 1
+ DIFF_NSPEC1D_RADIAL(2,2) = 1
+ endif
+ endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 2D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! exact number of surface elements for faces along XI and ETA
+
+do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
+ if(iter_region == IREGION_CRUST_MANTLE) then
+ ifirst_region = 1
+ ilast_region = 10 + layer_offset
+ else if(iter_region == IREGION_OUTER_CORE) then
+ ifirst_region = 11 + layer_offset
+ ilast_region = NUMBER_OF_MESH_LAYERS - 1
+ else if(iter_region == IREGION_INNER_CORE) then
+ ifirst_region = NUMBER_OF_MESH_LAYERS
+ ilast_region = NUMBER_OF_MESH_LAYERS
+ else
+ stop 'incorrect region code detected'
+ endif
+ tmp_sum_xi = 0
+ tmp_sum_eta = 0
+ do iter_layer = ifirst_region, ilast_region
+ if (this_region_has_a_doubling(iter_layer)) then
+ if (ner(iter_layer) == 1) then
+ nb_lay_sb = 1
+ nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK_1L
+ nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK_1L
+ else
+ nb_lay_sb = 2
+ nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK
+ nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK
+ endif
+ doubling = 1
+ else
+ doubling = 0
+ nb_lay_sb = 0
+ nspec2D_xi_sb = 0
+ nspec2D_eta_sb = 0
+ endif
+
+ tmp_sum_xi = tmp_sum_xi + ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)) * &
+ (ner(iter_layer) - doubling*nb_lay_sb)) + &
+ doubling * ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)) * (nspec2D_xi_sb/2))
+
+ tmp_sum_eta = tmp_sum_eta + ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)) * &
+ (ner(iter_layer) - doubling*nb_lay_sb)) + &
+ doubling * ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)) * (nspec2D_eta_sb/2))
+ enddo
+ NSPEC2D_XI(iter_region) = tmp_sum_xi
+ NSPEC2D_ETA(iter_region) = tmp_sum_eta
+ if (iter_region == IREGION_INNER_CORE .and. INCLUDE_CENTRAL_CUBE) then
+ NSPEC2D_XI(iter_region) = NSPEC2D_XI(iter_region) + &
+ ((NEX_PER_PROC_XI / ratio_divide_central_cube)*(NEX_XI / ratio_divide_central_cube))
+ NSPEC2D_ETA(iter_region) = NSPEC2D_ETA(iter_region) + &
+ ((NEX_PER_PROC_ETA / ratio_divide_central_cube)*(NEX_XI / ratio_divide_central_cube))
+ endif
+enddo
+
+! difference of number of surface elements along xi or eta for outer core if the superbrick is cut
+ DIFF_NSPEC2D_XI(:,:) = 0
+ DIFF_NSPEC2D_ETA(:,:) = 0
+ if (CUT_SUPERBRICK_XI) then
+ if (CUT_SUPERBRICK_ETA) then
+ DIFF_NSPEC2D_XI(2,1) = 2
+ DIFF_NSPEC2D_XI(1,2) = 2
+ DIFF_NSPEC2D_XI(2,3) = 2
+ DIFF_NSPEC2D_XI(1,4) = 2
+
+ DIFF_NSPEC2D_ETA(2,1) = 1
+ DIFF_NSPEC2D_ETA(2,2) = 1
+ DIFF_NSPEC2D_ETA(1,3) = 1
+ DIFF_NSPEC2D_ETA(1,4) = 1
+ else
+ DIFF_NSPEC2D_ETA(2,1) = 1
+ DIFF_NSPEC2D_ETA(1,2) = 1
+ endif
+ else
+ if (CUT_SUPERBRICK_ETA) then
+ DIFF_NSPEC2D_XI(2,1) = 2
+ DIFF_NSPEC2D_XI(1,2) = 2
+ endif
+ endif
+ DIFF_NSPEC2D_XI(:,:) = DIFF_NSPEC2D_XI(:,:) * (NEX_PER_PROC_XI / ratio_divide_central_cube)
+ DIFF_NSPEC2D_ETA(:,:) = DIFF_NSPEC2D_ETA(:,:) * (NEX_PER_PROC_ETA / ratio_divide_central_cube)
+
+! exact number of surface elements on the bottom and top boundaries
+
+! in the crust and mantle
+ NSPEC2D_TOP(IREGION_CRUST_MANTLE) = (NEX_XI/ratio_sampling_array(1))*(NEX_ETA/ratio_sampling_array(1))/NPROC
+ NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE) = (NEX_XI/ratio_sampling_array(10+layer_offset))*&
+ (NEX_ETA/ratio_sampling_array(10+layer_offset))/NPROC
+
+! in the outer core with mesh doubling
+ if (ADD_4TH_DOUBLING) then
+ NSPEC2D_TOP(IREGION_OUTER_CORE) = (NEX_XI/(ratio_divide_central_cube/4))*(NEX_ETA/(ratio_divide_central_cube/4))/NPROC
+ NSPEC2D_BOTTOM(IREGION_OUTER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
+ else
+ NSPEC2D_TOP(IREGION_OUTER_CORE) = (NEX_XI/(ratio_divide_central_cube/2))*(NEX_ETA/(ratio_divide_central_cube/2))/NPROC
+ NSPEC2D_BOTTOM(IREGION_OUTER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
+ endif
+
+! in the top of the inner core
+ NSPEC2D_TOP(IREGION_INNER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
+ NSPEC2D_BOTTOM(IREGION_INNER_CORE) = NSPEC2D_TOP(IREGION_INNER_CORE)
+
+! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX(:) = NSPEC2D_ETA(:)
+ NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE) = NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE) + maxval(DIFF_NSPEC2D_ETA(:,:))
+ NSPEC2DMAX_YMIN_YMAX(:) = NSPEC2D_XI(:)
+ NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE) = NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE) + maxval(DIFF_NSPEC2D_XI(:,:))
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 3D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! exact number of spectral elements in each region
+
+do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
+ if(iter_region == IREGION_CRUST_MANTLE) then
+ ifirst_region = 1
+ ilast_region = 10 + layer_offset
+ else if(iter_region == IREGION_OUTER_CORE) then
+ ifirst_region = 11 + layer_offset
+ ilast_region = NUMBER_OF_MESH_LAYERS - 1
+ else if(iter_region == IREGION_INNER_CORE) then
+ ifirst_region = NUMBER_OF_MESH_LAYERS
+ ilast_region = NUMBER_OF_MESH_LAYERS
+ else
+ stop 'incorrect region code detected'
+ endif
+ tmp_sum = 0;
+ do iter_layer = ifirst_region, ilast_region
+ if (this_region_has_a_doubling(iter_layer)) then
+ if (ner(iter_layer) == 1) then
+ nb_lay_sb = 1
+ nspec_sb = NSPEC_SUPERBRICK_1L
+ else
+ nb_lay_sb = 2
+ nspec_sb = NSPEC_DOUBLING_SUPERBRICK
+ endif
+ doubling = 1
+ else
+ doubling = 0
+ nb_lay_sb = 0
+ nspec_sb = 0
+ endif
+ tmp_sum = tmp_sum + ((NEX_XI / ratio_sampling_array(iter_layer)) * (NEX_ETA / ratio_sampling_array(iter_layer)) * &
+ (ner(iter_layer) - doubling*nb_lay_sb)) + &
+ doubling * ((NEX_XI / ratio_sampling_array(iter_layer)) * (NEX_ETA / ratio_sampling_array(iter_layer)) * &
+ (nspec_sb/4))
+ enddo
+ NSPEC(iter_region) = tmp_sum / NPROC
+enddo
+
+ if(INCLUDE_CENTRAL_CUBE) NSPEC(IREGION_INNER_CORE) = NSPEC(IREGION_INNER_CORE) + &
+ (NEX_PER_PROC_XI / ratio_divide_central_cube) * &
+ (NEX_PER_PROC_ETA / ratio_divide_central_cube) * &
+ (NEX_XI / ratio_divide_central_cube)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! calculation of number of points (NGLOB) below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 1D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! theoretical number of Gauss-Lobatto points in radial direction
+ NGLOB1D_RADIAL(:) = NSPEC1D_RADIAL(:)*(NGLLZ-1)+1
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 2D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! 2-D addressing and buffers for summation between slices
+! we add one to number of points because of the flag after the last point
+ NGLOB2DMAX_XMIN_XMAX(:) = NSPEC2DMAX_XMIN_XMAX(:)*NGLLY*NGLLZ + 1
+ NGLOB2DMAX_YMIN_YMAX(:) = NSPEC2DMAX_YMIN_YMAX(:)*NGLLX*NGLLZ + 1
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 3D case
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! exact number of global points in each region
+
+! initialize array
+ NGLOB(:) = 0
+
+! in the inner core (no doubling region + eventually central cube)
+ if(INCLUDE_CENTRAL_CUBE) then
+ NGLOB(IREGION_INNER_CORE) = ((NEX_PER_PROC_XI/ratio_divide_central_cube) &
+ *(NGLLX-1)+1)*((NEX_PER_PROC_ETA/ratio_divide_central_cube) &
+ *(NGLLY-1)+1)*((NER_TOP_CENTRAL_CUBE_ICB + NEX_XI / ratio_divide_central_cube)*(NGLLZ-1)+1)
+ else
+ NGLOB(IREGION_INNER_CORE) = ((NEX_PER_PROC_XI/ratio_divide_central_cube) &
+ *(NGLLX-1)+1)*((NEX_PER_PROC_ETA/ratio_divide_central_cube) &
+ *(NGLLY-1)+1)*((NER_TOP_CENTRAL_CUBE_ICB)*(NGLLZ-1)+1)
+ endif
+
+! in the crust-mantle and outercore
+ do iter_region = IREGION_CRUST_MANTLE,IREGION_OUTER_CORE
+ if(iter_region == IREGION_CRUST_MANTLE) then
+ ifirst_region = 1
+ ilast_region = 10 + layer_offset
+ else if(iter_region == IREGION_OUTER_CORE) then
+ ifirst_region = 11 + layer_offset
+ ilast_region = NUMBER_OF_MESH_LAYERS - 1
+ else
+ stop 'incorrect region code detected'
+ endif
+ tmp_sum = 0;
+ do iter_layer = ifirst_region, ilast_region
+ nglob_int_surf_eta=0
+ nglob_int_surf_xi=0
+ nglob_ext_surf = 0
+ nglob_center_edge = 0
+ nglob_corner_edge = 0
+ nglob_border_edge = 0
+ if (this_region_has_a_doubling(iter_layer)) then
+ if (iter_region == IREGION_OUTER_CORE .and. iter_layer == lastdoubling_layer .and. &
+ (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA)) then
+ doubling = 1
+ normal_doubling = 0
+ cut_doubling = 1
+ nb_lay_sb = 2
+ nglob_edge = 0
+ nglob_surf = 0
+ nglob_vol = 8*NGLLX**3 - 12*NGLLX**2 + 6*NGLLX - 1
+ nglob_int_surf_eta = 6*NGLLX**2 - 7*NGLLX + 2
+ nglob_int_surf_xi = 5*NGLLX**2 - 5*NGLLX + 1
+ nglob_ext_surf = 4*NGLLX**2-4*NGLLX+1
+ nglob_center_edge = 4*(NGLLX-1)+1
+ nglob_corner_edge = 2*(NGLLX-1)+1
+ nglob_border_edge = 3*(NGLLX-1)+1
+ else
+ if (ner(iter_layer) == 1) then
+ nb_lay_sb = 1
+ nglob_vol = 28*NGLLX**3 - 62*NGLLX**2 + 47*NGLLX - 12
+ nglob_surf = 6*NGLLX**2-8*NGLLX+3
+ nglob_edge = NGLLX
+ else
+ nb_lay_sb = 2
+ nglob_vol = 32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13
+ nglob_surf = 8*NGLLX**2-11*NGLLX+4
+ nglob_edge = 2*NGLLX-1
+ endif
+ doubling = 1
+ normal_doubling = 1
+ cut_doubling = 0
+ endif
+ padding = -1
+ else
+ doubling = 0
+ normal_doubling = 0
+ cut_doubling = 0
+ padding = 0
+ nb_lay_sb = 0
+ nglob_vol = 0
+ nglob_surf = 0
+ nglob_edge = 0
+ endif
+ if (iter_layer == ilast_region) padding = padding +1
+ nblocks_xi = NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)
+ nblocks_eta = NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)
+
+ tmp_sum = tmp_sum + &
+ ((nblocks_xi)*(NGLLX-1)+1) * ((nblocks_eta)*(NGLLX-1)+1) * ((ner(iter_layer) - doubling*nb_lay_sb)*(NGLLX-1)+padding)+&
+ normal_doubling * ((((nblocks_xi*nblocks_eta)/4)*nglob_vol) - &
+ (((nblocks_eta/2-1)*nblocks_xi/2+(nblocks_xi/2-1)*nblocks_eta/2)*nglob_surf) + &
+ ((nblocks_eta/2-1)*(nblocks_xi/2-1)*nglob_edge)) + &
+ cut_doubling*(nglob_vol*(nblocks_xi*nblocks_eta) - &
+ ( nblocks_eta*(int(nblocks_xi/2)*nglob_int_surf_xi + int((nblocks_xi-1)/2)*nglob_ext_surf) + &
+ nblocks_xi*(int(nblocks_eta/2)*nglob_int_surf_eta + int((nblocks_eta-1)/2)*nglob_ext_surf)&
+ ) + &
+ ( int(nblocks_xi/2)*int(nblocks_eta/2)*nglob_center_edge + &
+ int((nblocks_xi-1)/2)*int((nblocks_eta-1)/2)*nglob_corner_edge + &
+ ((int(nblocks_eta/2)*int((nblocks_xi-1)/2))+(int((nblocks_eta-1)/2)*int(nblocks_xi/2)))*nglob_border_edge&
+ ))
+ enddo
+ NGLOB(iter_region) = tmp_sum
+ enddo
+
+!!! example :
+!!! nblocks_xi/2=5
+!!! ____________________________________
+!!! I I I I I I
+!!! I I I I I I
+!!! I I I I I I
+!!! nblocks_eta/2=3 I______+______+______+______+______I
+!!! I I I I I I
+!!! I I I I I I
+!!! I I I I I I
+!!! I______+______+______+______+______I
+!!! I I I I I I
+!!! I I I I I I
+!!! I I I I I I
+!!! I______I______I______I______I______I
+!!!
+!!! NGLOB for this doubling layer = 3*5*Volume - ((3-1)*5+(5-1)*3)*Surface + (3-1)*(5-1)*Edge
+!!!
+!!! 32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13 -> nb GLL points in a superbrick (Volume)
+!!! 8*NGLLX**2-11*NGLLX+4 -> nb GLL points on a superbrick side (Surface)
+!!! 2*NGLLX-1 -> nb GLL points on a corner edge of a superbrick (Edge)
+
+!!! for the one layer superbrick :
+!!! NGLOB = 28.NGLL^3 - 62.NGLL^2 + 47.NGLL - 12 (Volume)
+!!! NGLOB = 6.NGLL^2 - 8.NGLL + 3 (Surface)
+!!! NGLOB = NGLL (Edge)
+!!!
+!!! those results were obtained by using the script UTILS/doubling_brick/count_nglob_analytical.pl
+!!! with an opendx file of the superbrick's geometry
+
+!!! for the basic doubling bricks (two layers)
+!!! NGLOB = 8.NGLL^3 - 12.NGLL^2 + 6.NGLL - 1 (VOLUME)
+!!! NGLOB = 5.NGLL^2 - 5.NGLL + 1 (SURFACE 1)
+!!! NGLOB = 6.NGLL^2 - 7.NGLL + 2 (SURFACE 2)
+
+ end subroutine read_compute_parameters
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2008-03-22 15:45:42 UTC (rev 11502)
@@ -490,3 +490,9 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
+! for Cuthill McKee permutation
+ logical, parameter :: PERFORM_CUTHILL_MCKEE = .true.
+ logical, parameter :: PERMUTE_INNER_CORE = .false.
+ integer, parameter :: NGNOD_HEXAHEDRA = 8
+ logical, parameter :: CMcK_MULTI = .true.
+ integer, parameter :: LIMIT_MULTI_CUTHILL=50
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -36,7 +36,7 @@
ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY, &
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE, &
- NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
+ NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,NGLOB2DMAX_XY, &
myrank,LOCAL_PATH,OCEANS,ibathy_topo, &
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD,&
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
@@ -256,7 +256,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
@@ -484,6 +484,8 @@
integer NUMBER_OF_MESH_LAYERS,layer_shift,cpt,first_layer_aniso,last_layer_aniso,FIRST_ELT_NON_ANISO
double precision, dimension(:,:), allocatable :: stretch_tab
+ integer :: NGLOB2DMAX_XY
+
integer :: nb_layer_above_aniso,FIRST_ELT_ABOVE_ANISO
integer, parameter :: maxker=200
@@ -535,6 +537,14 @@
double precision r_moho,r_400,r_670
logical :: is_superbrick
+! added for Cuthill McKee permutation
+ integer, dimension(:), allocatable :: perm,perm_tmp,temp_array_1D_int
+ logical, dimension(:,:), allocatable :: temp_array_2D_log
+ integer, dimension(:,:,:,:), allocatable :: temp_array_int
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: temp_array_real
+ double precision, dimension(:,:,:,:), allocatable :: temp_array_dble
+ double precision, dimension(:,:,:,:,:), allocatable :: temp_array_dble_5dim
+
! the height at which the central cube is cut
integer :: nz_inf_limit
@@ -1434,13 +1444,13 @@
! 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_ETA_FACE,iregion_code)
+ NSPEC2D_ETA_FACE,iregion_code,NGLOB2DMAX_XY,nglob)
call get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_XI_FACE,iregion_code)
+ NSPEC2D_XI_FACE,iregion_code,NGLOB2DMAX_XY,nglob)
call get_MPI_1D_buffers(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool,idoubling, &
xstore,ystore,zstore,ifseg,npointot, &
- NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion_code)
+ NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion_code,nglob)
! Stacey
if(NCHUNKS /= 6) &
@@ -1475,6 +1485,141 @@
! should be zero in all the regions except in the mantle
nspec_tiso = count(idoubling(1:nspec) == IFLAG_220_80) + count(idoubling(1:nspec) == IFLAG_80_MOHO)
+! ***************************************************
+! Cuthill McKee permutation
+! ***************************************************
+
+ if (iregion_code /= IREGION_INNER_CORE .or. PERMUTE_INNER_CORE) then
+ allocate(perm(nspec))
+ if(iregion_code == IREGION_CRUST_MANTLE) then
+ ! do not permute anisotropic elements
+ perm(1:FIRST_ELT_NON_ANISO-1) = (/ (i,i=1,FIRST_ELT_NON_ANISO-1) /)
+
+ ! no more connectivity between layers below and above the anisotropic layers => 2 permutations
+ ! permute the bottom of the region : below the aniso layers
+ allocate(perm_tmp(FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO))
+ call get_perm(ibool(:,:,:,FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1),perm_tmp,LIMIT_MULTI_CUTHILL,&
+(FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO),nglob,.true.,.false.)
+ perm(FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1) = perm_tmp(:)+(FIRST_ELT_NON_ANISO-1)
+ deallocate(perm_tmp)
+
+ ! permute the top of the region : above the aniso layers
+ allocate(perm_tmp(nspec-FIRST_ELT_ABOVE_ANISO+1))
+ call get_perm(ibool(:,:,:,FIRST_ELT_ABOVE_ANISO:nspec),perm_tmp,LIMIT_MULTI_CUTHILL,&
+(nspec-FIRST_ELT_ABOVE_ANISO+1),nglob,.true.,.false.)
+ perm(FIRST_ELT_ABOVE_ANISO:nspec) = perm_tmp(:)+(FIRST_ELT_ABOVE_ANISO-1)
+ deallocate(perm_tmp)
+ else
+ ! the 3 last parameters are : PERFORM_CUTHILL_MCKEE,INVERSE,FACE
+ call get_perm(ibool,perm,LIMIT_MULTI_CUTHILL,nspec,nglob,.true.,.false.)
+ endif
+
+ ! permutation of xstore, ystore, zstore, rhostore, kappavstore, kappahstore,
+ ! muvstore, muhstore, eta_anisostore, xixstore, xiystore, xizstore, etaxstore,
+ ! etaystore, etazstore, gammaxstore, gammaystore, gammazstore, no more jacobianstore
+
+ allocate(temp_array_dble(NGLLX,NGLLY,NGLLZ,nspec))
+ if(ATTENUATION .and. ATTENUATION_3D) then
+ call permute_elements_dble(Qmu_store,temp_array_dble,perm,nspec)
+ allocate(temp_array_dble_5dim(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
+ temp_array_dble_5dim(:,:,:,:,:) = tau_e_store(:,:,:,:,:)
+ do i = 1,nspec
+ tau_e_store(:,:,:,:,perm(i)) = temp_array_dble_5dim(:,:,:,:,i)
+ enddo
+ deallocate(temp_array_dble_5dim)
+ endif
+ call permute_elements_dble(xstore,temp_array_dble,perm,nspec)
+ call permute_elements_dble(ystore,temp_array_dble,perm,nspec)
+ call permute_elements_dble(zstore,temp_array_dble,perm,nspec)
+ deallocate(temp_array_dble)
+
+
+ allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
+ if(NCHUNKS /= 6) then
+ call permute_elements_real(rho_vp,temp_array_real,perm,nspec)
+ call permute_elements_real(rho_vs,temp_array_real,perm,nspec)
+ endif
+ if((ANISOTROPIC_INNER_CORE .and. iregion_code == IREGION_INNER_CORE) .or. &
+ (ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) then
+ call permute_elements_real(c11store,temp_array_real,perm,nspec)
+ call permute_elements_real(c12store,temp_array_real,perm,nspec)
+ call permute_elements_real(c13store,temp_array_real,perm,nspec)
+ call permute_elements_real(c14store,temp_array_real,perm,nspec)
+ call permute_elements_real(c15store,temp_array_real,perm,nspec)
+ call permute_elements_real(c16store,temp_array_real,perm,nspec)
+ call permute_elements_real(c22store,temp_array_real,perm,nspec)
+ call permute_elements_real(c23store,temp_array_real,perm,nspec)
+ call permute_elements_real(c24store,temp_array_real,perm,nspec)
+ call permute_elements_real(c25store,temp_array_real,perm,nspec)
+ call permute_elements_real(c26store,temp_array_real,perm,nspec)
+ call permute_elements_real(c33store,temp_array_real,perm,nspec)
+ call permute_elements_real(c34store,temp_array_real,perm,nspec)
+ call permute_elements_real(c35store,temp_array_real,perm,nspec)
+ call permute_elements_real(c36store,temp_array_real,perm,nspec)
+ call permute_elements_real(c44store,temp_array_real,perm,nspec)
+ call permute_elements_real(c45store,temp_array_real,perm,nspec)
+ call permute_elements_real(c46store,temp_array_real,perm,nspec)
+ call permute_elements_real(c55store,temp_array_real,perm,nspec)
+ call permute_elements_real(c56store,temp_array_real,perm,nspec)
+ call permute_elements_real(c66store,temp_array_real,perm,nspec)
+ endif
+ call permute_elements_real(rhostore,temp_array_real,perm,nspec)
+ call permute_elements_real(kappavstore,temp_array_real,perm,nspec)
+ call permute_elements_real(kappahstore,temp_array_real,perm,nspec)
+ call permute_elements_real(muvstore,temp_array_real,perm,nspec)
+ call permute_elements_real(muhstore,temp_array_real,perm,nspec)
+ call permute_elements_real(eta_anisostore,temp_array_real,perm,nspec)
+ call permute_elements_real(xixstore,temp_array_real,perm,nspec)
+ call permute_elements_real(xiystore,temp_array_real,perm,nspec)
+ call permute_elements_real(xizstore,temp_array_real,perm,nspec)
+ call permute_elements_real(etaxstore,temp_array_real,perm,nspec)
+ call permute_elements_real(etaystore,temp_array_real,perm,nspec)
+ call permute_elements_real(etazstore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammaxstore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammaystore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammazstore,temp_array_real,perm,nspec)
+ deallocate(temp_array_real)
+
+ ! permutation of ibool
+ allocate(temp_array_int(NGLLX,NGLLY,NGLLZ,nspec))
+ call permute_elements_integer(ibool,temp_array_int,perm,nspec)
+ deallocate(temp_array_int)
+
+ ! permutation of iMPIcut_*
+ allocate(temp_array_2D_log(2,nspec))
+ temp_array_2D_log(:,:) = iMPIcut_xi(:,:)
+ do i = 1,nspec
+ iMPIcut_xi(:,perm(i)) = temp_array_2D_log(:,i)
+ enddo
+ temp_array_2D_log(:,:) = iMPIcut_eta(:,:)
+ do i = 1,nspec
+ iMPIcut_eta(:,perm(i)) = temp_array_2D_log(:,i)
+ enddo
+ deallocate(temp_array_2D_log)
+
+ ! permutation of iboun
+ allocate(temp_array_2D_log(6,nspec))
+ temp_array_2D_log(:,:) = iboun(:,:)
+ do i = 1,nspec
+ iboun(:,perm(i)) = temp_array_2D_log(:,i)
+ enddo
+ deallocate(temp_array_2D_log)
+
+ ! permutation of idoubling
+ allocate(temp_array_1D_int(nspec))
+ temp_array_1D_int(:) = idoubling(:)
+ do i = 1,nspec
+ idoubling(perm(i)) = temp_array_1D_int(i)
+ enddo
+ deallocate(temp_array_1D_int)
+
+ deallocate(perm)
+ endif
+
+! ***************************************************
+! end of Cuthill McKee permutation
+! ***************************************************
+
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, &
@@ -1658,7 +1803,7 @@
jacobian2D_xmin,jacobian2D_xmax, &
jacobian2D_ymin,jacobian2D_ymax, &
jacobian2D_bottom,jacobian2D_top, &
- iMPIcut_xi,iMPIcut_eta,nspec,nglob, &
+ nspec,nglob, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -6,7 +6,8 @@
! Main authors: Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory, California Institute of Technology, USA
! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau, April 2007
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
@@ -26,7 +27,7 @@
subroutine get_MPI_1D_buffers(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
idoubling,xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion)
+ NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion,nglob_ori)
! routine to create the MPI 1D chunk buffers for edges
@@ -34,7 +35,7 @@
include "constants.h"
- integer nspec,myrank,iregion
+ integer nspec,myrank,nglob_ori,nglob,ipoin1D,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER
logical iMPIcut_xi(2,nspec)
@@ -61,6 +62,27 @@
! processor identification
character(len=150) prname
+! arrays for sorting routine
+ integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: work
+ integer, dimension(:), allocatable :: ibool_selected
+ double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
+
+! allocate arrays for message buffers with maximum size
+! define maximum size for message buffers
+ allocate(ibool_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(xstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(ystore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(zstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(ind(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(ninseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(iglob(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(locval(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(ifseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(iwork(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+ allocate(work(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
+
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -83,37 +105,40 @@
ispeccount=0
do ispec=1,nspec
+ ! remove central cube for chunk buffers
+ if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+ ! corner detection here
+ if(iMPIcut_xi(1,ispec) .and. iMPIcut_eta(1,ispec)) then
+ ispeccount=ispeccount+1
+ ! loop on all the points
+ ix = 1
+ iy = 1
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin1D = npoin1D + 1
+ ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
+ endif
+ enddo
-! remove central cube for chunk buffers
- if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-! corner detection here
- if(iMPIcut_xi(1,ispec) .and. iMPIcut_eta(1,ispec)) then
-
- ispeccount=ispeccount+1
-
-! loop on all the points
- ix = 1
- iy = 1
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin1D = npoin1D + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ do ipoin1D=1,npoin1D
+ write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
+ ystore_selected(ipoin1D),zstore_selected(ipoin1D)
enddo
- endif
- enddo
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -139,39 +164,41 @@
! nb of elements in this 1D buffer
ispeccount=0
-
do ispec=1,nspec
+ ! remove central cube for chunk buffers
+ if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
+ idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+ ! corner detection here
+ if(iMPIcut_xi(2,ispec) .and. iMPIcut_eta(1,ispec)) then
+ ispeccount=ispeccount+1
+ ! loop on all the points
+ ix = NGLLX
+ iy = 1
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin1D = npoin1D + 1
+ ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
+ endif
+ enddo
-! remove central cube for chunk buffers
- if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
- idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-! corner detection here
- if(iMPIcut_xi(2,ispec) .and. iMPIcut_eta(1,ispec)) then
-
- ispeccount=ispeccount+1
-
-! loop on all the points
- ix = NGLLX
- iy = 1
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin1D = npoin1D + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ do ipoin1D=1,npoin1D
+ write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
+ ystore_selected(ipoin1D),zstore_selected(ipoin1D)
enddo
- endif
- enddo
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -220,20 +247,28 @@
iy = NGLLY
do iz=1,NGLLZ
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin1D = npoin1D + 1
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin1D = npoin1D + 1
+ ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
+ endif
+ enddo
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+ do ipoin1D=1,npoin1D
+ write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
+ ystore_selected(ipoin1D),zstore_selected(ipoin1D)
enddo
- endif
- enddo
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -278,20 +313,28 @@
iy = NGLLY
do iz=1,NGLLZ
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin1D = npoin1D + 1
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin1D = npoin1D + 1
+ ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
+ endif
+ enddo
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+ do ipoin1D=1,npoin1D
+ write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
+ ystore_selected(ipoin1D),zstore_selected(ipoin1D)
enddo
- endif
- enddo
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -304,5 +347,17 @@
if(ispeccount /= NSPEC1D_RADIAL_CORNER(iregion,3) .or. npoin1D /= NGLOB1D_RADIAL_CORNER(iregion,3)) &
call exit_MPI(myrank,'error MPI 1D buffer detection in xi=right')
+ deallocate(ibool_selected)
+ deallocate(xstore_selected)
+ deallocate(ystore_selected)
+ deallocate(zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
end subroutine get_MPI_1D_buffers
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -6,7 +6,8 @@
! Main authors: Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory, California Institute of Technology, USA
! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau, April 2007
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
@@ -26,7 +27,7 @@
subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_XI_FACE,iregion)
+ NSPEC2D_XI_FACE,iregion,NGLOB2DMAX_XY,nglob_ori)
! this routine detects cut planes along eta
! In principle the left cut plane of the first slice
@@ -37,9 +38,10 @@
include "constants.h"
- integer nspec,myrank,iregion
+ integer nspec,myrank,nglob_ori,nglob,ipoin2D,NGLOB2DMAX_XY,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE
+
logical iMPIcut_eta(2,nspec)
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
@@ -62,9 +64,31 @@
! processor identification
character(len=150) prname
+! arrays for sorting routine
+ integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: work
+ integer, dimension(:), allocatable :: ibool_selected
+ double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
+
+
+! allocate arrays for message buffers with maximum size
+! define maximum size for message buffers
+ allocate(ibool_selected(NGLOB2DMAX_XY))
+ allocate(xstore_selected(NGLOB2DMAX_XY))
+ allocate(ystore_selected(NGLOB2DMAX_XY))
+ allocate(zstore_selected(NGLOB2DMAX_XY))
+ allocate(ind(NGLOB2DMAX_XY))
+ allocate(ninseg(NGLOB2DMAX_XY))
+ allocate(iglob(NGLOB2DMAX_XY))
+ allocate(locval(NGLOB2DMAX_XY))
+ allocate(ifseg(NGLOB2DMAX_XY))
+ allocate(iwork(NGLOB2DMAX_XY))
+ allocate(work(NGLOB2DMAX_XY))
+
! theoretical number of surface elements in the buffers
! cut planes along eta=constant correspond to XI faces
- nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
+ nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -86,28 +110,33 @@
ispecc1=0
do ispec=1,nspec
- if(iMPIcut_eta(1,ispec)) then
-
- ispecc1=ispecc1+1
-
-! loop on all the points in that 2-D element, including edges
- iy = 1
- do ix=1,NGLLX
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin2D_eta = npoin2D_eta + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ if(iMPIcut_eta(1,ispec)) then
+ ispecc1=ispecc1+1
+ ! loop on all the points in that 2-D element, including edges
+ iy = 1
+ do ix=1,NGLLX
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin2D_eta = npoin2D_eta + 1
+ ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
enddo
+ endif
enddo
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+
+ do ipoin2D=1,npoin2D_eta
+ write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
+ ystore_selected(ipoin2D),zstore_selected(ipoin2D)
enddo
! put flag to indicate end of the list of points
@@ -124,7 +153,7 @@
!
! determine if the element falls on the right MPI cut plane
!
- nspec2Dtheor = NSPEC2D_XI_FACE(iregion,2)
+ nspec2Dtheor = NSPEC2D_XI_FACE(iregion,2)
! global point number and coordinates right MPI cut-plane
open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown')
@@ -139,28 +168,33 @@
ispecc2=0
do ispec=1,nspec
- if(iMPIcut_eta(2,ispec)) then
-
- ispecc2=ispecc2+1
-
-! loop on all the points in that 2-D element, including edges
- iy = NGLLY
- do ix=1,NGLLX
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin2D_eta = npoin2D_eta + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ if(iMPIcut_eta(2,ispec)) then
+ ispecc2=ispecc2+1
+ ! loop on all the points in that 2-D element, including edges
+ iy = NGLLY
+ do ix=1,NGLLX
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin2D_eta = npoin2D_eta + 1
+ ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
enddo
+ endif
enddo
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+
+ do ipoin2D=1,npoin2D_eta
+ write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
+ ystore_selected(ipoin2D),zstore_selected(ipoin2D)
enddo
! put flag to indicate end of the list of points
@@ -174,5 +208,17 @@
! compare number of surface elements detected to analytical value
if(ispecc2 /= nspec2Dtheor) call exit_MPI(myrank,'error MPI cut-planes detection in eta=right')
+ deallocate(ibool_selected)
+ deallocate(xstore_selected)
+ deallocate(ystore_selected)
+ deallocate(zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
end subroutine get_MPI_cutplanes_eta
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -6,7 +6,8 @@
! Main authors: Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory, California Institute of Technology, USA
! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau, April 2007
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
@@ -26,7 +27,7 @@
subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_ETA_FACE,iregion)
+ NSPEC2D_ETA_FACE,iregion,NGLOB2DMAX_XY,nglob_ori)
! this routine detects cut planes along xi
! In principle the left cut plane of the first slice
@@ -37,7 +38,7 @@
include "constants.h"
- integer nspec,myrank,iregion
+ integer nspec,myrank,nglob_ori,nglob,ipoin2D,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_ETA_FACE
logical iMPIcut_xi(2,nspec)
@@ -62,10 +63,32 @@
! processor identification
character(len=150) prname,errmsg
+! arrays for sorting routine
+ integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: work
+ integer NGLOB2DMAX_XY
+ integer, dimension(:), allocatable :: ibool_selected
+ double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
+
+! allocate arrays for message buffers with maximum size
+! define maximum size for message buffers
+ allocate(ibool_selected(NGLOB2DMAX_XY))
+ allocate(xstore_selected(NGLOB2DMAX_XY))
+ allocate(ystore_selected(NGLOB2DMAX_XY))
+ allocate(zstore_selected(NGLOB2DMAX_XY))
+ allocate(ind(NGLOB2DMAX_XY))
+ allocate(ninseg(NGLOB2DMAX_XY))
+ allocate(iglob(NGLOB2DMAX_XY))
+ allocate(locval(NGLOB2DMAX_XY))
+ allocate(ifseg(NGLOB2DMAX_XY))
+ allocate(iwork(NGLOB2DMAX_XY))
+ allocate(work(NGLOB2DMAX_XY))
+
+
! theoretical number of surface elements in the buffers
! cut planes along xi=constant correspond to ETA faces
- nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
-
+ nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -86,28 +109,33 @@
ispecc1=0
do ispec=1,nspec
- if(iMPIcut_xi(1,ispec)) then
-
- ispecc1=ispecc1+1
-
-! loop on all the points in that 2-D element, including edges
- ix = 1
- do iy=1,NGLLY
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin2D_xi = npoin2D_xi + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ if(iMPIcut_xi(1,ispec)) then
+ ispecc1=ispecc1+1
+ ! loop on all the points in that 2-D element, including edges
+ ix = 1
+ do iy=1,NGLLY
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin2D_xi = npoin2D_xi + 1
+ ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
enddo
+ endif
enddo
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+
+ do ipoin2D=1,npoin2D_xi
+ write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
+ ystore_selected(ipoin2D),zstore_selected(ipoin2D)
enddo
! put flag to indicate end of the list of points
@@ -126,7 +154,7 @@
!
! determine if the element falls on the right MPI cut plane
!
- nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,2)
+ nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,2)
! global point number and coordinates right MPI cut-plane
open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='unknown')
@@ -141,30 +169,36 @@
ispecc2=0
do ispec=1,nspec
- if(iMPIcut_xi(2,ispec)) then
-
- ispecc2=ispecc2+1
-
-! loop on all the points in that 2-D element, including edges
- ix = NGLLX
- do iy=1,NGLLY
- do iz=1,NGLLZ
-
-! select point, if not already selected
- if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
- mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
- npoin2D_xi = npoin2D_xi + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
+ if(iMPIcut_xi(2,ispec)) then
+ ispecc2=ispecc2+1
+ ! loop on all the points in that 2-D element, including edges
+ ix = NGLLX
+ do iy=1,NGLLY
+ do iz=1,NGLLZ
+ ! select point, if not already selected
+ if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
+ mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
+ npoin2D_xi = npoin2D_xi + 1
+ ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
+ xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
+ ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
+ zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
+ endif
+ enddo
enddo
+ endif
enddo
- endif
+ nglob=nglob_ori
+ call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
+ ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
+
+ do ipoin2D=1,npoin2D_xi
+ write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
+ ystore_selected(ipoin2D),zstore_selected(ipoin2D)
enddo
+
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -179,5 +213,17 @@
call exit_MPI(myrank,errmsg)
endif
+ deallocate(ibool_selected)
+ deallocate(xstore_selected)
+ deallocate(ystore_selected)
+ deallocate(zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
end subroutine get_MPI_cutplanes_xi
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -42,7 +42,7 @@
include "constants.h"
- integer nspec,myrank
+ integer nspec,myrank,dummy_var,ispec_tmp
integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
@@ -83,6 +83,12 @@
double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+! arrays for sorting routine
+ integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: work
+ double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
+
! check that the parameter file is correct
if(NGNOD /= 27) call exit_MPI(myrank,'elements should have 27 control nodes')
if(NGNOD2D /= 9) call exit_MPI(myrank,'surface elements should have 9 control nodes')
@@ -136,6 +142,38 @@
call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmin,normal_xmin,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+
+ allocate (xstore_selected(ispecb1))
+ allocate (ystore_selected(ispecb1))
+ allocate (zstore_selected(ispecb1))
+ allocate(ind(ispecb1))
+ allocate(ninseg(ispecb1))
+ allocate(iglob(ispecb1))
+ allocate(locval(ispecb1))
+ allocate(ifseg(ispecb1))
+ allocate(iwork(ispecb1))
+ allocate(work(ispecb1))
+
+ do ispec_tmp=1,ispecb1
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_xmin(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_xmin(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_xmin(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb1,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_xmin,iglob,normal_xmin,jacobian2D_xmin,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLY,NGLLZ)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
endif
! on boundary: xmax
@@ -176,6 +214,38 @@
call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmax,normal_xmax,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+
+ allocate (xstore_selected(ispecb2))
+ allocate (ystore_selected(ispecb2))
+ allocate (zstore_selected(ispecb2))
+ allocate(ind(ispecb2))
+ allocate(ninseg(ispecb2))
+ allocate(iglob(ispecb2))
+ allocate(locval(ispecb2))
+ allocate(ifseg(ispecb2))
+ allocate(iwork(ispecb2))
+ allocate(work(ispecb2))
+
+ do ispec_tmp=1,ispecb2
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_xmax(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_xmax(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_xmax(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb2,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_xmax,iglob,normal_xmax,jacobian2D_xmax,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLY,NGLLZ)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
endif
! on boundary: ymin
@@ -216,6 +286,38 @@
call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymin,normal_ymin,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+
+ allocate (xstore_selected(ispecb3))
+ allocate (ystore_selected(ispecb3))
+ allocate (zstore_selected(ispecb3))
+ allocate(ind(ispecb3))
+ allocate(ninseg(ispecb3))
+ allocate(iglob(ispecb3))
+ allocate(locval(ispecb3))
+ allocate(ifseg(ispecb3))
+ allocate(iwork(ispecb3))
+ allocate(work(ispecb3))
+
+ do ispec_tmp=1,ispecb3
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_ymin(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_ymin(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_ymin(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb3,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_ymin,iglob,normal_ymin,jacobian2D_ymin,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLX,NGLLZ)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
endif
! on boundary: ymax
@@ -256,6 +358,38 @@
call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymax,normal_ymax,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+
+ allocate (xstore_selected(ispecb4))
+ allocate (ystore_selected(ispecb4))
+ allocate (zstore_selected(ispecb4))
+ allocate(ind(ispecb4))
+ allocate(ninseg(ispecb4))
+ allocate(iglob(ispecb4))
+ allocate(locval(ispecb4))
+ allocate(ifseg(ispecb4))
+ allocate(iwork(ispecb4))
+ allocate(work(ispecb4))
+
+ do ispec_tmp=1,ispecb4
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_ymax(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_ymax(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_ymax(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb4,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_ymax,iglob,normal_ymax,jacobian2D_ymax,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLX,NGLLZ)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
endif
! on boundary: bottom
@@ -295,6 +429,38 @@
call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, &
jacobian2D_bottom,normal_bottom,NGLLX,NGLLY,NSPEC2D_BOTTOM)
+
+ allocate (xstore_selected(ispecb5))
+ allocate (ystore_selected(ispecb5))
+ allocate (zstore_selected(ispecb5))
+ allocate(ind(ispecb5))
+ allocate(ninseg(ispecb5))
+ allocate(iglob(ispecb5))
+ allocate(locval(ispecb5))
+ allocate(ifseg(ispecb5))
+ allocate(iwork(ispecb5))
+ allocate(work(ispecb5))
+
+ do ispec_tmp=1,ispecb5
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_bottom(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_bottom(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_bottom(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb5,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_bottom,iglob,normal_bottom,jacobian2D_bottom,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLX,NGLLY)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
+
endif
! on boundary: top
@@ -334,10 +500,42 @@
call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
+
+ allocate (xstore_selected(ispecb6))
+ allocate (ystore_selected(ispecb6))
+ allocate (zstore_selected(ispecb6))
+ allocate(ind(ispecb6))
+ allocate(ninseg(ispecb6))
+ allocate(iglob(ispecb6))
+ allocate(locval(ispecb6))
+ allocate(ifseg(ispecb6))
+ allocate(iwork(ispecb6))
+ allocate(work(ispecb6))
+
+ do ispec_tmp=1,ispecb6
+ xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm_top(ispec_tmp))
+ ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm_top(ispec_tmp))
+ zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm_top(ispec_tmp))
+ enddo
+
+ call sort_array_coordinates_gjb(ispecb6,xstore_selected,ystore_selected,zstore_selected, &
+ ibelm_top,iglob,normal_top,jacobian2D_top,locval,ifseg,dummy_var,ind,ninseg,iwork,work,NGLLX,NGLLY)
+
+ deallocate (xstore_selected)
+ deallocate (ystore_selected)
+ deallocate (zstore_selected)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iglob)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(iwork)
+ deallocate(work)
endif
enddo
+
! check theoretical value of elements at the bottom
if(ispecb5 /= NSPEC2D_BOTTOM) then
call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM')
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -0,0 +1,773 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_perm (ibool,perm,limit,nspec,nglob,INVERSE,FACE)
+
+ implicit none
+
+ include "constants.h"
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+! include "OUTPUT_FILES/values_from_mesher.h"
+
+ logical :: INVERSE,FACE
+
+! input
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+! output
+ integer, dimension(nspec) :: perm
+
+! local variables
+ integer nspec,nglob_GLL_full
+
+! a neighbor of a hexahedral node is a hexahedra which share a face whith it -> max degre of a node = 6
+ integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 100
+
+! global corner numbers that need to be created
+ integer, dimension(nglob) :: global_corner_number
+
+ integer mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+ integer, dimension(:), allocatable :: ne,np,adj
+ integer xadj(nspec+1)
+
+! arrays to store the permutation and inverse permutation of the Cuthill-McKee algorithm
+ integer, dimension(nspec) :: invperm
+
+ logical maskel(nspec)
+
+ integer i,istart,istop,number_of_neighbors
+
+ integer nglob_eight_corners_only,nglob
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only
+ integer total_size_ne,total_size_adj,limit
+
+!
+!-----------------------------------------------------------------------
+!
+ if(PERFORM_CUTHILL_MCKEE) then
+ ! total number of points in the mesh
+ nglob_GLL_full = nglob
+
+ !---- call Charbel Farhat's routines
+ write(IMAIN,*) 'form_elt_connectivity_foelco'
+ call form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,nglob_GLL_full,ibool,nglob_eight_corners_only)
+ do i=1,nspec
+ istart = mp(i)
+ istop = mp(i+1) - 1
+ enddo
+
+ allocate(np(nglob_eight_corners_only+1))
+ count_only = .true.
+ total_size_ne = 1
+ write(IMAIN,*) 'form_node_connectivity_fonoco 1'
+ allocate(ne(total_size_ne))
+ call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
+ deallocate(ne)
+ allocate(ne(total_size_ne))
+ count_only = .false.
+ write(IMAIN,*) 'form_node_connectivity_fonoco 2'
+ call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
+ do i=1,nglob_eight_corners_only
+ istart = np(i)
+ istop = np(i+1) - 1
+ enddo
+
+ count_only = .true.
+ total_size_adj = 1
+ write(IMAIN,*) 'create_adjacency_table_adjncy 1'
+ allocate(adj(total_size_adj))
+ call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,FACE)
+ deallocate(adj)
+ allocate(adj(total_size_adj))
+ count_only = .false.
+ write(IMAIN,*) 'create_adjacency_table_adjncy 2'
+ call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,FACE)
+ do i=1,nspec
+ istart = xadj(i)
+ istop = xadj(i+1) - 1
+ number_of_neighbors = istop-istart+1
+ if(number_of_neighbors < 1 .or. number_of_neighbors > MAX_NUMBER_OF_NEIGHBORS) stop 'incorrect number of neighbors'
+ enddo
+ deallocate(ne,np)
+
+! call the Cuthill-McKee sorting algorithm
+ write(IMAIN,*) 'launch cuthill_mckee'
+ call cuthill_mckee(adj,xadj,perm,invperm,nspec,total_size_adj,limit,INVERSE)
+ deallocate(adj)
+ write(IMAIN,*) 'permutation achieved'
+ else
+! create identity permutation in order to do nothing
+ do i=1,nspec
+ perm(i) = i
+ enddo
+ endif
+
+ end subroutine get_perm
+
+!
+!----------------------------------------------------
+!
+
+!=======================================================================
+!
+! Charbel Farhat's FEM topology routines
+!
+! Dimitri Komatitsch, February 1996 - Code based on Farhat's original version from 1987
+!
+! modified and adapted by Dimitri Komatitsch, May 2006
+!
+!=======================================================================
+
+ subroutine form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,&
+nglob_GLL_full,ibool,nglob_eight_corners_only)
+
+!-----------------------------------------------------------------------
+!
+! Forms the MN and MP arrays
+!
+! Input :
+! -------
+! ibool Array needed to build the element connectivity table
+! nspec Number of elements in the domain
+! NGNOD_HEXAHEDRA number of nodes per hexahedron (brick with 8 corners)
+!
+! Output :
+! --------
+! MN, MP This is the element connectivity array pair.
+! Array MN contains the list of the element
+! connectivity, that is, the nodes contained in each
+! element. They are stored in a stacked fashion.
+!
+! Pointer array MP stores the location of each
+! element list. Its length is equal to the number
+! of elements plus one.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,nglob_GLL_full
+
+! arrays with mesh parameters per slice
+ integer, intent(in), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! global corner numbers that need to be created
+ integer, intent(out), dimension(nglob_GLL_full) :: global_corner_number
+ integer, intent(out) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+ integer, intent(out) :: nglob_eight_corners_only
+
+ integer ninter,nsum,ispec,node,k,inumcorner,ix,iy,iz
+
+ ninter = 1
+ nsum = 1
+ mp(1) = 1
+
+!---- define topology of the elements in the mesh
+!---- we need to define adjacent numbers from the sub-mesh consisting of the corners only
+ nglob_eight_corners_only = 0
+ global_corner_number(:) = -1
+
+ do ispec=1,nspec
+
+ inumcorner = 0
+ do iz = 1,NGLLZ,NGLLZ-1
+ do iy = 1,NGLLY,NGLLY-1
+ do ix = 1,NGLLX,NGLLX-1
+
+ inumcorner = inumcorner + 1
+ if(inumcorner > NGNOD_HEXAHEDRA) stop 'corner number too large'
+
+! check if this point was already assigned a number previously, otherwise create one and store it
+ if(global_corner_number(ibool(ix,iy,iz,ispec)) == -1) then
+ nglob_eight_corners_only = nglob_eight_corners_only + 1
+ global_corner_number(ibool(ix,iy,iz,ispec)) = nglob_eight_corners_only
+ endif
+
+ node = global_corner_number(ibool(ix,iy,iz,ispec))
+ do k=nsum,ninter-1
+ if(node == mn(k)) goto 200
+ enddo
+
+ mn(ninter) = node
+ ninter = ninter + 1
+ 200 continue
+
+ enddo
+ enddo
+ enddo
+
+ nsum = ninter
+ mp(ispec + 1) = nsum
+
+ enddo
+
+ end subroutine form_elt_connectivity_foelco
+
+!
+!----------------------------------------------------
+!
+
+ subroutine form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,&
+nspec,count_only,total_size_ne)
+
+!-----------------------------------------------------------------------
+!
+! Forms the NE and NP arrays
+!
+! Input :
+! -------
+! MN, MP, nspec
+! nglob_eight_corners_only Number of nodes in the domain
+!
+! Output :
+! --------
+! NE, NP This is the node-connected element array pair.
+! Integer array NE contains a list of the
+! elements connected to each node, stored in stacked fashion.
+!
+! Array NP is the pointer array for the
+! location of a node's element list in the NE array.
+! Its length is equal to the number of points plus one.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only
+ integer total_size_ne
+
+ integer nglob_eight_corners_only,nspec
+
+ integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+
+ integer, intent(out) :: ne(total_size_ne),np(nglob_eight_corners_only+1)
+
+ integer nsum,inode,ispec,j
+
+ nsum = 1
+ np(1) = 1
+
+ do inode=1,nglob_eight_corners_only
+ do 200 ispec=1,nspec
+
+ do j=mp(ispec),mp(ispec + 1) - 1
+ if (mn(j) == inode) then
+ if(count_only) then
+ total_size_ne = nsum
+ else
+ ne(nsum) = ispec
+ endif
+ nsum = nsum + 1
+ goto 200
+ endif
+ enddo
+ 200 continue
+
+ np(inode + 1) = nsum
+
+ enddo
+
+ end subroutine form_node_connectivity_fonoco
+
+!
+!----------------------------------------------------
+!
+
+ subroutine create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,face)
+
+!-----------------------------------------------------------------------
+!
+! Establishes the element adjacency information of the mesh
+! Two elements are considered adjacent if they share a face.
+!
+! Input :
+! -------
+! MN, MP, NE, NP, nspec
+! MASKEL logical mask (length = nspec)
+!
+! Output :
+! --------
+! ADJ, XADJ This is the element adjacency array pair. Array
+! ADJ contains the list of the elements adjacent to
+! element i. They are stored in a stacked fashion.
+! Pointer array XADJ stores the location of each element list.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only,face
+ integer total_size_ne,total_size_adj
+
+ integer nglob_eight_corners_only
+
+ integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1),ne(total_size_ne),np(nglob_eight_corners_only+1)
+
+ integer, intent(out) :: adj(total_size_adj),xadj(nspec+1)
+
+ logical maskel(nspec)
+ integer countel(nspec)
+
+ integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
+
+ xadj(1) = 1
+ iad = 1
+
+ do ispec=1,nspec
+
+! reset mask
+ maskel(:) = .false.
+
+! mask current element
+ maskel(ispec) = .true.
+ if (face) countel(:) = 0
+
+ istart = mp(ispec)
+ istop = mp(ispec+1) - 1
+ do ino=istart,istop
+ node = mn(ino)
+ jstart = np(node)
+ jstop = np(node + 1) - 1
+ do 120 jel=jstart,jstop
+ nelem = ne(jel)
+ if(maskel(nelem)) goto 120
+ if (face) then
+ ! if 2 elements share at least 3 corners, therefore they share a face
+ countel(nelem) = countel(nelem) + 1
+ if (countel(nelem)>=3) then
+ if(count_only) then
+ total_size_adj = iad
+ else
+ adj(iad) = nelem
+ endif
+ maskel(nelem) = .true.
+ iad = iad + 1
+ endif
+ else
+ if(count_only) then
+ total_size_adj = iad
+ else
+ adj(iad) = nelem
+ endif
+ maskel(nelem) = .true.
+ iad = iad + 1
+ endif
+ 120 continue
+ enddo
+
+ xadj(ispec+1) = iad
+
+ enddo
+
+ end subroutine create_adjacency_table_adjncy
+
+ subroutine cuthill_mckee(adj,xadj,mask,invperm_all,nspec,total_size_adj,limit,INVERSE)
+
+ implicit none
+ include "constants.h"
+
+ integer, intent(in) :: nspec,total_size_adj, limit
+ integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
+
+ integer, intent(out), dimension(nspec) :: mask,invperm_all
+! integer, dimension(limit) :: invperm_sub
+ integer, dimension(nspec) :: invperm_sub
+ logical :: INVERSE
+ integer ispec,gsize,counter,nspec_sub,root,total_ordered_elts, next_root
+
+
+! fill the mask with ones
+ mask(:) = 1
+ invperm_all(:) = 0
+ counter = 0
+ nspec_sub = limit
+ root = 1
+ total_ordered_elts = 0
+
+ do while(total_ordered_elts < nspec)
+ ! création of a sublist of sorted elements which fit in the cache (the criterion of size is limit)
+ ! limit = nb of element that can fit in the L2 cache
+ call Cut_McK( root, nspec, total_size_adj, xadj, adj, mask, gsize, invperm_sub, limit, nspec_sub, next_root)
+ ! add the sublist in the main permutation list
+ invperm_all(total_ordered_elts+1:total_ordered_elts+nspec_sub) = invperm_sub(1:nspec_sub)
+ total_ordered_elts = total_ordered_elts + nspec_sub
+ ! seek for a new root to build the new sublist
+ if (next_root > 0) then
+ root = next_root
+ else
+! if (total_ordered_elts /= nspec) call find_next_root(next_root,xadj,adj,nspec,total_size_adj,mask)
+ if (total_ordered_elts /= nspec) &
+ call find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
+ root = next_root
+ endif
+ enddo
+
+ if (INVERSE) then
+ do ispec=1,nspec
+ mask(invperm_all(ispec)) = ispec
+ enddo
+ else
+ mask(:) = invperm_all(:)
+ endif
+
+ end subroutine cuthill_mckee
+
+
+!*******************************************************************************
+! Objective: Cuthill McKee ordering
+! The algorithm is:
+!
+! X(1) = ROOT.
+! for ( I = 1 to N-1)
+! Find all unlabeled neighbors of X(I),
+! assign them the next available labels, in order of increasing degree.
+!
+! Parameters:
+! root the starting point for the cm ordering.
+! nbnodes the number of nodes.
+! nnz the number of adjacency entries.
+!
+! xadj/adj the graph
+! mask only those nodes with nonzero mask are considered
+!
+! gsize the number of the connected component
+! invp Inverse invputation (from new order to old order)
+!*******************************************************************************
+
+
+subroutine find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
+ implicit none
+ include "constants.h"
+
+! input
+ integer, intent(in) :: total_size_adj,total_ordered_elts,nspec
+ integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
+ integer, intent(in), dimension(nspec) :: mask,invperm_all
+! output
+ integer, intent(out) :: next_root
+! variables
+ integer :: cur_node,neighbor_node,i,j
+
+ do i=total_ordered_elts, 1, -1
+ cur_node = invperm_all(i)
+ do j= xadj(cur_node), xadj(cur_node+1)-1
+ neighbor_node = adj(j)
+ if (mask(neighbor_node)/=0) then
+ next_root=neighbor_node
+ return
+ endif
+ enddo
+ enddo
+
+end subroutine find_next_root
+
+!*******************************************************************************
+! Objective: Cuthill McKee ordering
+! The algorithm is:
+!
+! X(1) = ROOT.
+! for ( I = 1 to N-1)
+! Find all unlabeled neighbors of X(I),
+! assign them the next available labels, in order of increasing degree.
+!
+! Parameters:
+! root the starting point for the cm ordering.
+! nbnodes the number of nodes.
+! nnz the number of adjacency entries.
+!
+! xadj/adj the graph
+! mask only those nodes with nonzero mask are considered
+!
+! gsize the number of the connected component
+! invp Inverse invputation (from new order to old order)
+!*******************************************************************************
+
+subroutine Cut_McK( root, nbnodes, nnz, xadj, adj, mask, gsize, invp, limit, nspec_sub, next_root)
+ implicit none
+ include "constants.h"
+!--------------------------------------------------------------- Input Variables
+ integer root, nnz, nbnodes, limit, nspec_sub, next_root
+
+ integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
+
+!-------------------------------------------------------------- Output Variables
+ integer gsize
+ integer invp(nbnodes)
+
+!--------------------------------------------------------------- Local Variables
+ integer i, j, k, l, lbegin, lnbr, linvp, lvlend, nbr, node, fnbr
+ integer deg(nbnodes)
+
+!------------------------------------------------------------------------- BEGIN
+! Find the degrees of the nodes in the subgraph specified by mask and root
+! Here invp is used to store a levelization of the subgraph
+ invp(:)=0
+ deg(:)=0
+ call degree ( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, invp)
+
+ mask(root) = 0
+
+ IF (gsize > 1) THEN
+ !If there is at least 2 nodes in the subgraph
+ lvlend = 0
+ lnbr = 1
+
+ DO while (lvlend < lnbr)
+ !lbegin/lvlend point to the begin/end of the present level
+ lbegin = lvlend + 1
+ lvlend = lnbr
+
+ do i= lbegin, lvlend
+ node = invp(i)
+
+ !Find the unnumbered neighbours of node.
+ !fnbr/lnbr point to the first/last neighbors of node
+ fnbr = lnbr + 1
+ do j= xadj(node), xadj(node+1)-1
+ nbr = adj(j)
+
+ if (mask(nbr) /= 0) then
+ lnbr = lnbr + 1
+ mask(nbr) = 0
+ invp(lnbr) = nbr
+ endif
+ enddo
+
+ !If no neighbors, go to next node in this level.
+ IF (lnbr > fnbr) THEN
+ !Sort the neighbors of NODE in increasing order by degree.
+ !Linear insertion is used.
+ k = fnbr
+ do while (k < lnbr)
+ l = k
+ k = k + 1
+ nbr = invp(k)
+
+ DO WHILE (fnbr < l)
+ linvp = invp(l)
+
+ if (deg(linvp) <= deg(nbr)) then
+ exit
+ endif
+
+ invp(l+1) = linvp
+ l = l-1
+ ENDDO
+
+ invp(l+1) = nbr
+ enddo
+ ENDIF
+ enddo
+ ENDDO
+
+ ENDIF
+
+ if (gsize > limit) then
+ do i = limit + 1 , nbnodes
+ node=invp(i)
+ if (node /=0) mask(node) = 1
+ enddo
+ next_root = invp(limit +1)
+ nspec_sub = limit
+ else
+ next_root = -1
+ nspec_sub = gsize
+ endif
+
+!--------------------------------------------------------------------------- END
+END subroutine Cut_McK
+
+
+!*******************************************************************************
+! Objective: computes the degrees of the nodes in the connected graph
+!
+! Parameters:
+! root the root node
+! nbnodes the number of nodes in the graph
+! nnz the graph size
+! xadj/adj the whole graph
+! mask Only nodes with mask == 0 are considered
+!
+! gsize the number of nodes in the connected graph
+! deg degree for all the nodes in the connected graph
+! level levelization of the connected graph
+!
+!*******************************************************************************
+
+subroutine degree( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, level )
+ implicit none
+!--------------------------------------------------------------- Input Variables
+ integer root, nbnodes, nnz
+ integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
+
+!-------------------------------------------------------------- Output Variables
+ integer gsize
+ integer deg(nbnodes), level(nbnodes)
+
+!--------------------------------------------------------------- Local Variables
+ integer i, j, ideg, lbegin, lvlend, lvsize, nxt, nbr, node
+
+!------------------------------------------------------------------------- BEGIN
+! The sign of xadj(I) is used to indicate if node i has been considered
+ xadj(root) = -xadj(root)
+ level(1) = root
+ nxt = 1
+ lvlend = 0
+ lvsize = 1
+
+ DO WHILE (lvsize > 0)
+ ! lbegin/lvlend points the begin/end of the present level
+ lbegin = lvlend + 1
+ lvlend = nxt
+
+ ! Find the degrees of nodes in the present level and generate the next level
+ DO i= lbegin, lvlend
+ node = level(i)
+ ideg = 0
+ do j= ABS( xadj(node) ), ABS( xadj(node+1) )-1
+ nbr = adj(j)
+
+ if (mask(nbr) /= 0) then
+ ideg = ideg + 1
+
+ if (xadj(nbr) >= 0) then
+ xadj(nbr) = -xadj(nbr)
+ nxt = nxt + 1
+ level(nxt) = nbr
+ endif
+ endif
+ enddo
+
+ deg(node) = ideg
+ ENDDO
+
+ !Compute the level size of the next level
+ lvsize = nxt - lvlend
+ ENDDO
+
+ !Reset xadj to its correct sign
+ do i = 1, nxt
+ node = level(i)
+ xadj(node) = -xadj(node)
+ enddo
+
+ gsize = nxt
+!--------------------------------------------------------------------------- END
+END subroutine degree
+
+
+ subroutine permute_elements_real(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ real(kind=CUSTOM_REAL), intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_real
+
+!
+!-----------------------------------------------------------------------
+!
+
+! implement permutation of elements for arrays of integer type
+ subroutine permute_elements_integer(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ integer, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_integer
+
+!
+!-----------------------------------------------------------------------
+!
+
+! implement permutation of elements for arrays of double precision type
+ subroutine permute_elements_dble(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ double precision, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_dble
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -366,7 +366,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
@@ -1098,7 +1098,7 @@
call MPI_BCAST(D3MM_V%qq,3*(NK+1)*(NK+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
elseif(THREE_D_MODEL == THREE_D_MODEL_SEA99_JP3D) then
! the variables read are declared and stored in structure SEA99M_V and JP3DM_V
- if(myrank == 0) then
+ if(myrank == 0) then
call read_sea99_s_model(SEA99M_V)
call read_iso3d_dpzhao_model(JP3DM_V)
endif
@@ -1406,6 +1406,7 @@
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE, &
NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
+ max(NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code)), &
myrank,LOCAL_PATH,OCEANS,ibathy_topo, &
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -42,7 +42,7 @@
jacobian2D_xmin,jacobian2D_xmax, &
jacobian2D_ymin,jacobian2D_ymax, &
jacobian2D_bottom,jacobian2D_top, &
- iMPIcut_xi,iMPIcut_eta,nspec,nglob, &
+ nspec,nglob, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source, &
@@ -144,9 +144,6 @@
! number of elements on the boundaries
integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
-! MPI cut-planes parameters along xi and along eta
- logical iMPIcut_xi(2,nspec),iMPIcut_eta(2,nspec)
-
integer i,j,k,ispec,iglob,nspec1, nglob1
! attenuation
@@ -362,15 +359,6 @@
close(27)
-! MPI cut-planes parameters along xi and along eta
- open(unit=27,file=prname(1:len_trim(prname))//'iMPIcut_xi.bin',status='unknown',form='unformatted')
- write(27) iMPIcut_xi
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'iMPIcut_eta.bin',status='unknown',form='unformatted')
- write(27) iMPIcut_eta
- close(27)
-
if(ATTENUATION .and. ATTENUATION_3D) then
open(unit=27, file=prname(1:len_trim(prname))//'attenuation3D.bin', status='unknown', form='unformatted')
write(27) tau_s
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/sort_array_coordinates.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/sort_array_coordinates.f90 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/sort_array_coordinates.f90 2008-03-22 15:45:42 UTC (rev 11502)
@@ -232,3 +232,164 @@
end subroutine swap_all_buffers
+ subroutine sort_array_coordinates_gjb(npointot,x,y,z,ibool,iglob,normal,jacobian_2D,loc,ifseg,nglob,ind,ninseg,iwork,work,dim1,dim2)
+
+! this routine MUST be in double precision to avoid sensitivity
+! to roundoff errors in the coordinates of the points
+
+ implicit none
+
+ include "constants.h"
+
+ integer npointot,nglob,dim1,dim2
+
+ real(kind=CUSTOM_REAL) jacobian_2D(dim1,dim2,npointot), normal(NDIM,dim1,dim2,npointot)
+ integer ibool(npointot),iglob(npointot),loc(npointot)
+ integer ind(npointot),ninseg(npointot)
+ logical ifseg(npointot)
+ double precision x(npointot),y(npointot),z(npointot)
+ integer iwork(npointot)
+ double precision work(npointot)
+
+ integer ipoin,i,j
+ integer nseg,ioff,iseg,ig
+ double precision xtol
+
+! establish initial pointers
+ do ipoin=1,npointot
+ loc(ipoin)=ipoin
+ enddo
+
+! define a tolerance, normalized radius is 1., so let's use a small value
+ xtol = SMALLVALTOL
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=npointot
+
+ do j=1,NDIM
+
+! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
+
+ call rank_buffers(x(ioff),ind,ninseg(iseg))
+
+ else if(j == 2) then
+
+ call rank_buffers(y(ioff),ind,ninseg(iseg))
+
+ else
+
+ call rank_buffers(z(ioff),ind,ninseg(iseg))
+
+ endif
+
+ call swap_all_buffers_gjb(ibool(ioff),loc(ioff), &
+ x(ioff),y(ioff),z(ioff),normal(:,:,:,ioff),jacobian_2D(:,:,ioff),iwork,work,ind,ninseg(iseg),dim1,dim2)
+
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+! check for jumps in current coordinate
+ if(j == 1) then
+ do i=2,npointot
+ if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,npointot
+ if(dabs(y(i)-y(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,npointot
+ if(dabs(z(i)-z(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ endif
+
+! count up number of different segments
+ nseg=0
+ do i=1,npointot
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+! assign global node numbers (now sorted lexicographically)
+ ig=0
+ do i=1,npointot
+ if(ifseg(i)) ig=ig+1
+ iglob(loc(i))=ig
+ enddo
+
+ nglob=ig
+
+ end subroutine sort_array_coordinates_gjb
+
+ subroutine swap_all_buffers_gjb(IA,IB,A,B,C,NORM,JAC,IW,W,ind,n,dim1,dim2)
+!
+! swap arrays IA, IB, A, B and C according to addressing in array IND
+!
+ implicit none
+ include "constants.h"
+
+ integer n,dim1,dim2
+
+ integer IND(n)
+ integer IA(n),IB(n),IW(n)
+ double precision A(n),B(n),C(n),W(n)
+ real(kind=CUSTOM_REAL) NORM(NDIM,dim1,dim2,n), JAC(dim1,dim2,n), WN(NDIM,dim1,dim2,n), WJ(dim1,dim2,n)
+
+ integer i
+
+ do i=1,n
+ W(i)=A(i)
+ IW(i)=IA(i)
+ enddo
+
+ do i=1,n
+ A(i)=W(ind(i))
+ IA(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=B(i)
+ IW(i)=IB(i)
+ enddo
+
+ do i=1,n
+ B(i)=W(ind(i))
+ IB(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=C(i)
+ enddo
+
+ do i=1,n
+ C(i)=W(ind(i))
+ enddo
+
+ do i=1,n
+ WN(:,:,:,i)=NORM(:,:,:,i)
+ enddo
+ do i=1,n
+ NORM(:,:,:,i)=WN(:,:,:,ind(i))
+ enddo
+
+ do i=1,n
+ WJ(:,:,i)=JAC(:,:,i)
+ enddo
+ do i=1,n
+ JAC(:,:,i)=WJ(:,:,ind(i))
+ enddo
+
+ end subroutine swap_all_buffers_gjb
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt 2008-03-22 15:41:50 UTC (rev 11501)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt 2008-03-22 15:45:42 UTC (rev 11502)
@@ -4,10 +4,8 @@
Things that could be done in a future version 4.1:
-- Cuthill-McKee sorting to reduce cache misses
+- partial ALTIVEC support (will be done by Nicolas Le Goff)
-- maybe partial ALTIVEC support (if useful)
-
- use a potential of (rho * u) instead of u in the fluid, in case of
fluid-fluid discontinuities
@@ -150,3 +148,13 @@
- supprimer sections qui decrivent write_AVS_mesh_quality, check_buffers*.f90, check_mesh_quality*.f90 etc du manuel une fois que le mesher et le solver auront ete fusionnes
+---------------------------------------------------------------------------------
+---------------------------------------------------------------------------------
+---------------------------------------------------------------------------------
+---------------------------------------------------------------------------------
+
+Done:
+-----
+
+- Cuthill-McKee sorting to reduce cache misses: finished by David Michea in March 2008
+
More information about the cig-commits
mailing list