[cig-commits] r22467 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared trunk/src/shared
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Jun 30 18:22:14 PDT 2013
Author: dkomati1
Date: 2013-06-30 18:22:13 -0700 (Sun, 30 Jun 2013)
New Revision: 22467
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/Makefile
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_helpers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_manager.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_elements.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_points.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/define_all_layers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_timestep_and_layers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_c_binary.c
Log:
done merging new files in "shared"
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/Makefile 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,56 @@
+#=====================================================================
+#
+# S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+# --------------------------------------------------
+#
+# Main authors: Dimitri Komatitsch and Jeroen Tromp
+# Princeton University, USA
+# and University of Pau / CNRS / INRIA, France
+# (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+# April 2011
+#
+# 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.
+#
+#=====================================================================
+
+DIR = shared
+
+# The rest of this file is generic
+#######################################
+
+####
+#### targets
+####
+
+default:
+ $(MAKE) -C ../.. $(DIR)
+
+all:
+ $(MAKE) -C ../.. all
+
+clean:
+ $(MAKE) -C ../.. CLEAN=$(DIR) clean
+
+cleanall:
+ $(MAKE) -C ../.. clean
+
+backup:
+ mkdir -p bak
+ cp *f90 *h Makefile bak
+
+bak: backup
+
+.PHONY: default all clean cleanall backup bak
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,170 @@
+#=====================================================================
+#
+# S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+# --------------------------------------------------
+#
+# Main authors: Dimitri Komatitsch and Jeroen Tromp
+# Princeton University, USA
+# and University of Pau / CNRS / INRIA, France
+# (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+# April 2011
+#
+# 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.
+#
+#=====================================================================
+
+#######################################
+
+shared_TARGETS = \
+ $(shared_OBJECTS) \
+ $(EMPTY_MACRO)
+
+shared_OBJECTS = \
+ $O/auto_ner.o \
+ $O/broadcast_compute_parameters.o \
+ $O/calendar.o \
+ $O/count_number_of_sources.o \
+ $O/create_name_database.o \
+ $O/create_serial_name_database.o \
+ $O/euler_angles.o \
+ $O/exit_mpi.o \
+ $O/force_ftz.o \
+ $O/get_model_parameters.o \
+ $O/get_value_parameters.o \
+ $O/gll_library.o \
+ $O/hex_nodes.o \
+ $O/intgrl.o \
+ $O/lagrange_poly.o \
+ $O/make_ellipticity.o \
+ $O/make_gravity.o \
+ $O/memory_eval.o \
+ $O/model_prem.o \
+ $O/model_topo_bathy.o \
+ $O/param_reader.o \
+ $O/read_compute_parameters.o \
+ $O/read_parameter_file.o \
+ $O/read_value_parameters.o \
+ $O/reduce.o \
+ $O/rthetaphi_xyz.o \
+ $O/save_header_file.o \
+ $O/spline_routines.o \
+ $O/write_c_binary.o \
+ $(EMPTY_MACRO)
+
+#######################################
+
+## compilation directories
+S := ${S_TOP}/src/shared
+$(shared_OBJECTS): S = ${S_TOP}/src/shared
+
+####
+#### rule for each .o file below
+####
+
+##
+## shared
+##
+$O/auto_ner.o: ${SETUP}/constants.h $S/auto_ner.f90
+ ${FCCOMPILE_CHECK} -c -o $O/auto_ner.o ${FCFLAGS_f90} $S/auto_ner.f90
+
+$O/broadcast_compute_parameters.o: ${SETUP}/constants.h $S/broadcast_compute_parameters.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/broadcast_compute_parameters.o ${FCFLAGS_f90} $S/broadcast_compute_parameters.f90
+
+$O/calendar.o: $S/calendar.f90
+ ${FCCOMPILE_CHECK} -c -o $O/calendar.o ${FCFLAGS_f90} $S/calendar.f90
+
+$O/count_number_of_sources.o: ${SETUP}/constants.h $S/count_number_of_sources.f90
+ ${FCCOMPILE_CHECK} -c -o $O/count_number_of_sources.o ${FCFLAGS_f90} $S/count_number_of_sources.f90
+
+$O/create_name_database.o: ${SETUP}/constants.h $S/create_name_database.f90
+ ${FCCOMPILE_CHECK} -c -o $O/create_name_database.o ${FCFLAGS_f90} $S/create_name_database.f90
+
+$O/create_serial_name_database.o: ${SETUP}/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/euler_angles.o: ${SETUP}/constants.h $S/euler_angles.f90
+ ${FCCOMPILE_CHECK} -c -o $O/euler_angles.o ${FCFLAGS_f90} $S/euler_angles.f90
+
+### C compilation
+$O/force_ftz.o: $S/force_ftz.c ${SETUP}/config.h
+ ${CC} -c $(CPPFLAGS) $(CFLAGS) -o $O/force_ftz.o $S/force_ftz.c
+
+$O/get_model_parameters.o: ${SETUP}/constants.h $S/get_model_parameters.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_model_parameters.o ${FCFLAGS_f90} $S/get_model_parameters.f90
+
+$O/get_value_parameters.o: ${SETUP}/constants.h $S/get_value_parameters.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_value_parameters.o ${FCFLAGS_f90} $S/get_value_parameters.f90
+
+$O/gll_library.o: ${SETUP}/constants.h $S/gll_library.f90
+ ${FCCOMPILE_CHECK} -c -o $O/gll_library.o ${FCFLAGS_f90} $S/gll_library.f90
+
+$O/hex_nodes.o: ${SETUP}/constants.h $S/hex_nodes.f90
+ ${FCCOMPILE_CHECK} -c -o $O/hex_nodes.o ${FCFLAGS_f90} $S/hex_nodes.f90
+
+$O/intgrl.o: ${SETUP}/constants.h $S/intgrl.f90
+ ${FCCOMPILE_CHECK} -c -o $O/intgrl.o ${FCFLAGS_f90} $S/intgrl.f90
+
+$O/lagrange_poly.o: ${SETUP}/constants.h $S/lagrange_poly.f90
+ ${FCCOMPILE_CHECK} -c -o $O/lagrange_poly.o ${FCFLAGS_f90} $S/lagrange_poly.f90
+
+$O/make_ellipticity.o: ${SETUP}/constants.h $S/make_ellipticity.f90
+ ${FCCOMPILE_CHECK} -c -o $O/make_ellipticity.o ${FCFLAGS_f90} $S/make_ellipticity.f90
+
+$O/make_gravity.o: ${SETUP}/constants.h $S/make_gravity.f90
+ ${FCCOMPILE_CHECK} -c -o $O/make_gravity.o ${FCFLAGS_f90} $S/make_gravity.f90
+
+$O/memory_eval.o: ${SETUP}/constants.h $S/memory_eval.f90
+ ${FCCOMPILE_CHECK} -c -o $O/memory_eval.o ${FCFLAGS_f90} $S/memory_eval.f90
+
+$O/model_prem.o: ${SETUP}/constants.h $S/model_prem.f90
+ ${FCCOMPILE_CHECK} -c -o $O/model_prem.o ${FCFLAGS_f90} $S/model_prem.f90
+
+$O/param_reader.o: $S/param_reader.c ${SETUP}/config.h
+ ${CC} -c $(CPPFLAGS) $(CFLAGS) -o $O/param_reader.o $S/param_reader.c
+
+$O/read_compute_parameters.o: ${SETUP}/constants.h $S/read_compute_parameters.f90
+ ${FCCOMPILE_CHECK} -c -o $O/read_compute_parameters.o ${FCFLAGS_f90} $S/read_compute_parameters.f90
+
+$O/read_parameter_file.o: ${SETUP}/constants.h $S/read_parameter_file.f90
+ ${FCCOMPILE_CHECK} -c -o $O/read_parameter_file.o ${FCFLAGS_f90} $S/read_parameter_file.f90
+
+$O/read_value_parameters.o: ${SETUP}/constants.h $S/read_value_parameters.f90
+ ${FCCOMPILE_CHECK} -c -o $O/read_value_parameters.o ${FCFLAGS_f90} $S/read_value_parameters.f90
+
+$O/reduce.o: ${SETUP}/constants.h $S/reduce.f90
+ ${FCCOMPILE_CHECK} -c -o $O/reduce.o ${FCFLAGS_f90} $S/reduce.f90
+
+$O/rthetaphi_xyz.o: ${SETUP}/constants.h $S/rthetaphi_xyz.f90
+ ${FCCOMPILE_CHECK} -c -o $O/rthetaphi_xyz.o ${FCFLAGS_f90} $S/rthetaphi_xyz.f90
+
+$O/save_header_file.o: ${SETUP}/constants.h $S/save_header_file.f90
+ ${FCCOMPILE_CHECK} -c -o $O/save_header_file.o ${FCFLAGS_f90} $S/save_header_file.f90
+
+$O/spline_routines.o: ${SETUP}/constants.h $S/spline_routines.f90
+ ${FCCOMPILE_CHECK} -c -o $O/spline_routines.o ${FCFLAGS_f90} $S/spline_routines.f90
+
+### C compilation
+$O/write_c_binary.o: $S/write_c_binary.c ${SETUP}/config.h
+ $(CC) $(CPPFLAGS) $(CFLAGS) -c -o $O/write_c_binary.o $S/write_c_binary.c
+
+##
+## shared objects with mpi compilation
+##
+$O/exit_mpi.o: ${SETUP}/constants.h $S/exit_mpi.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/exit_mpi.o ${FCFLAGS_f90} $S/exit_mpi.f90
+
+$O/model_topo_bathy.o: ${SETUP}/constants.h $S/model_topo_bathy.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/model_topo_bathy.o ${FCFLAGS_f90} $S/model_topo_bathy.f90
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_c_binary.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_c_binary.c 2013-07-01 01:17:07 UTC (rev 22466)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_c_binary.c 2013-07-01 01:22:13 UTC (rev 22467)
@@ -6,8 +6,8 @@
!
! Main authors: Dimitri Komatitsch and Jeroen Tromp
! Princeton University, USA
- ! and University of Pau / CNRS / INRIA, France
- ! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+ ! and CNRS / INRIA / University of Pau, France
+ ! (c) Princeton University and CNRS / INRIA / University of Pau
! April 2011
!
! This program is free software; you can redistribute it and/or modify
@@ -280,7 +280,10 @@
FILE *ft;
int itemlen,remlen,donelen,ret;
- char *buf;
+// DK DK fixed the warning we got when compiling with Intel icc
+// DK DK solution found at http://osdir.com/ml/network.quagga.devel/2004-09/msg00090.html
+// DK DK void *buf;
+ char *buf = NULL;
// file pointer
ft = fp_abs[*fid];
@@ -317,7 +320,10 @@
FILE *ft;
int ret,itemlen,remlen,donelen;
long long pos;
- char *buf;
+// DK DK fixed the warning we got when compiling with Intel icc
+// DK DK solution found at http://osdir.com/ml/network.quagga.devel/2004-09/msg00090.html
+// DK DK void *buf;
+ char *buf = NULL;
// file pointer
ft = fp_abs[*fid];
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_helpers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_helpers.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_helpers.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,356 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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.
+!
+!=====================================================================
+
+
+!-------------------------------------------------------------------------------
+!> \file adios_helpers.f90
+!! \brief Helpers to set up adios features.
+!! \author MPBL
+!-------------------------------------------------------------------------------
+
+!===============================================================================
+!> Get the ADIOS error message from an adios error number if there is an error.
+!! \param adios_err The error code considered.
+subroutine check_adios_err(myrank, adios_err)
+ use adios_read_mod
+ implicit none
+ integer, intent(in) :: myrank, adios_err
+ character(len=1024) :: msg
+
+ if (adios_err /= 0) then
+ call adios_errmsg(msg)
+ print *, "process: ", myrank, ", error: ", msg
+ stop
+ endif
+end subroutine check_adios_err
+
+
+!===============================================================================
+!> Define an ADIOS scalar double precision variable and autoincrement
+!! the adios group size by (8).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+subroutine define_adios_double_scalar (adios_group, name, path, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 6 == real(kind=8)
+ call adios_define_var (adios_group, name, path, 6, "", "", "", varid)
+ group_size_inc = group_size_inc + 8
+end subroutine define_adios_double_scalar
+
+!===============================================================================
+!> Define an ADIOS scalar integer variable and autoincrement the adios
+!! group size by (4).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+subroutine define_adios_integer_scalar (adios_group, name, path, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 2 == integer(kind=4)
+ call adios_define_var (adios_group, name, path, adios_integer, "", "", "", varid)
+ group_size_inc = group_size_inc + 4
+end subroutine define_adios_integer_scalar
+
+!===============================================================================
+!> Define an ADIOS scalar byte variable and autoincrement the adios
+!! group size by (1).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+subroutine define_adios_byte_scalar (adios_group, name, path, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 0 == byte == any_data_type(kind=1)
+ call adios_define_var (adios_group, name, path, 0, "", "", "", varid)
+ group_size_inc = group_size_inc + 1
+end subroutine define_adios_byte_scalar
+
+!===============================================================================
+!> Define a local ADIOS array of integers and autoincrement the adios
+!! group size by (4 * number of elements).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param dim The number of elements in the 1D array. Required to
+!! correctly increment adios group size.
+!! \param dim_str The "stringified" version of dim. Needed by adios
+!! to define variables
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+subroutine define_adios_integer_local_array1D (adios_group, name, path, dim, dim_str, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path, dim_str
+ integer(kind=8), intent(inout) :: group_size_inc
+ integer, intent(in) :: dim
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 2 == integer
+ call adios_define_var (adios_group, name, path, 2, dim_str, "", "", varid)
+ group_size_inc = group_size_inc + 4*dim
+end subroutine define_adios_integer_local_array1D
+
+!===============================================================================
+!> Define a local ADIOS array of doubles and autoincrement the adios
+!! group size by (8 * number of elements).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param dim The number of elements in the 1D array. Required to
+!! correctly increment adios group size.
+!! \param dim_str The "stringified" version of dim. Needed by adios
+!! to define variables
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+subroutine define_adios_double_local_array1D (adios_group, name, path, dim, dim_str, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path, dim_str
+ integer(kind=8), intent(inout) :: group_size_inc
+ integer, intent(in) :: dim
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 6 == real(kind=8)
+ call adios_define_var (adios_group, name, path, 6, dim_str, "", "", varid)
+ group_size_inc = group_size_inc + 8*dim
+end subroutine define_adios_double_local_array1D
+
+!===============================================================================
+!> Define a local ADIOS string and autoincrement the adios
+!! group size by (1 * string's length).
+!! \param adios_group The adios group where the variables belongs
+!! \param name The variable to be defined
+!! \param path The logical path structuring the data and containing
+!! the variable
+!! \param len The length of the string(number of character. in Fortran
+!! it does not include a final '\0' -- null -- character)
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+!! \note Adios string are scalar values counting for (1) byte. It is
+!! mandatory to increase the group size by the length of the
+!! string in order not to overlap 'data regions'.
+subroutine define_adios_string (adios_group, name, path, length, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Arguments
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: name, path
+ integer(kind=8), intent(inout) :: group_size_inc
+ integer :: length
+ ! Local Variables
+ integer(kind=8) :: varid ! dummy variable, adios use var name
+
+ ! adios: 9 == string
+ call adios_define_var (adios_group, name, path, 9, "", "", "", varid)
+ group_size_inc = group_size_inc + 1*length
+end subroutine define_adios_string
+
+!===============================================================================
+!> Define a global ADIOS 1D real array and autoincrement the adios
+!! group size.
+!! \param adios_group The adios group where the variables belongs
+!! \param array_name The variable to be defined. This is actually the path for
+!! ADIOS. The values are stored in array_name/array
+!! \param len The local dimension of the array.
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+!! \note This function define local, global and offset sizes as well as the
+!! array to store the values in.
+subroutine define_adios_global_real_1d_array(adios_group, array_name, &
+ local_dim, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Parameters
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: array_name
+ integer, intent(in) :: local_dim
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Variables
+ integer(kind=8) :: var_id
+
+ call define_adios_integer_scalar (adios_group, "local_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "global_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "offset", array_name, &
+ group_size_inc)
+! call adios_define_var(adios_group, "array", array_name, 6, &
+! "local_dim", "global_dim", "offset", var_id)
+ call adios_define_var(adios_group, "array", array_name, 5, &
+ array_name // "/local_dim", array_name // "/global_dim", &
+ array_name // "/offset", var_id)
+ group_size_inc = group_size_inc + local_dim*4
+end subroutine define_adios_global_real_1d_array
+
+!===============================================================================
+!> Define a global ADIOS 1D integer array and autoincrement the adios
+!! group size.
+!! \param adios_group The adios group where the variables belongs
+!! \param array_name The variable to be defined. This is actually the path for
+!! ADIOS. The values are stored in array_name/array
+!! \param len The local dimension of the array.
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+!! \note This function define local, global and offset sizes as well as the
+!! array to store the values in.
+subroutine define_adios_global_integer_1d_array(adios_group, array_name, &
+ local_dim, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Parameters
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: array_name
+ integer, intent(in) :: local_dim
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Variables
+ integer(kind=8) :: var_id
+
+ call define_adios_integer_scalar (adios_group, "local_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "global_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "offset", array_name, &
+ group_size_inc)
+! call adios_define_var(adios_group, "array", array_name, 6, &
+! "local_dim", "global_dim", "offset", var_id)
+ call adios_define_var(adios_group, "array", array_name, 2, &
+ array_name // "/local_dim", array_name // "/global_dim", &
+ array_name // "/offset", var_id)
+ group_size_inc = group_size_inc + local_dim*4
+end subroutine define_adios_global_integer_1d_array
+
+!===============================================================================
+!> Define a global ADIOS 1D logical array and autoincrement the adios
+!! group size.
+!! \param adios_group The adios group where the variables belongs
+!! \param array_name The variable to be defined. This is actually the path for
+!! ADIOS. The values are stored in array_name/array
+!! \param len The local dimension of the array.
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+!! \note This function define local, global and offset sizes as well as the
+!! array to store the values in.
+subroutine define_adios_global_logical_1d_array(adios_group, array_name, &
+ local_dim, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Parameters
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: array_name
+ integer, intent(in) :: local_dim
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Variables
+ integer(kind=8) :: var_id
+
+ call define_adios_integer_scalar (adios_group, "local_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "global_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "offset", array_name, &
+ group_size_inc)
+! call adios_define_var(adios_group, "array", array_name, 6, &
+! "local_dim", "global_dim", "offset", var_id)
+ call adios_define_var(adios_group, "array", array_name, 1, &
+ array_name // "/local_dim", array_name // "/global_dim", &
+ array_name // "/offset", var_id)
+ group_size_inc = group_size_inc + local_dim*1
+end subroutine define_adios_global_logical_1d_array
+
+!===============================================================================
+!> Define a global ADIOS 1D real array and autoincrement the adios
+!! group size.
+!! \param adios_group The adios group where the variables belongs
+!! \param array_name The variable to be defined. This is actually the path for
+!! ADIOS. The values are stored in array_name/array
+!! \param len The local dimension of the array.
+!! \param group_size_inc The inout adios group size to increment
+!! with the size of the variable
+!! \note This function define local, global and offset sizes as well as the
+!! array to store the values in.
+subroutine define_adios_global_double_1d_array(adios_group, array_name, &
+ local_dim, group_size_inc)
+ use adios_write_mod
+ implicit none
+ ! Parameters
+ integer(kind=8), intent(in) :: adios_group
+ character(len=*), intent(in) :: array_name
+ integer, intent(in) :: local_dim
+ integer(kind=8), intent(inout) :: group_size_inc
+ ! Variables
+ integer(kind=8) :: var_id
+
+ call define_adios_integer_scalar (adios_group, "local_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "global_dim", array_name, &
+ group_size_inc)
+ call define_adios_integer_scalar (adios_group, "offset", array_name, &
+ group_size_inc)
+ call adios_define_var(adios_group, "array", array_name, 6, &
+ array_name // "/local_dim", array_name // "/global_dim", &
+ array_name // "/offset", var_id)
+ group_size_inc = group_size_inc + local_dim*8
+end subroutine define_adios_global_double_1d_array
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_manager.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_manager.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_manager.F90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,51 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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.
+!
+!=====================================================================
+
+!> @brief Initialize ADIOS and setup the xml output file
+subroutine adios_setup()
+ use adios_write_mod, only: adios_init
+
+ implicit none
+ integer :: adios_err, sizeMB
+
+ call adios_init_noxml (adios_err);
+ sizeMB = 200 ! TODO 200MB is surely not the right size for the adios buffer
+ call adios_allocate_buffer (sizeMB , adios_err)
+end subroutine adios_setup
+
+!> @brief Finalize ADIOS. Must be called once everything is written down.
+subroutine adios_cleanup()
+ use mpi
+ use adios_write_mod, only: adios_finalize
+
+ implicit none
+ integer :: myrank
+ integer :: adios_err, ierr
+
+ call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)
+ call adios_finalize (myrank, adios_err)
+end subroutine adios_cleanup
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_elements.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_elements.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_elements.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,375 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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 count_elements(NEX_XI,NEX_ETA,NEX_PER_PROC_XI,NPROC,&
+ NEX_PER_PROC_ETA,ratio_divide_central_cube,&
+ NSPEC,NSPEC2D_XI,NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ ner,ratio_sampling_array,this_region_has_a_doubling, &
+ ifirst_region,ilast_region,iter_region,iter_layer, &
+ doubling,tmp_sum,tmp_sum_xi,tmp_sum_eta, &
+ NUMBER_OF_MESH_LAYERS,layer_offset,nspec2D_xi_sb,nspec2D_eta_sb, &
+ nb_lay_sb, nspec_sb, nglob_surf, &
+ CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA, INCLUDE_CENTRAL_CUBE, &
+ last_doubling_layer, &
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ tmp_sum_nglob2D_xi, tmp_sum_nglob2D_eta,divider,nglob_edges_h,&
+ nglob_edge_v,to_remove)
+
+
+ implicit none
+
+ include "constants.h"
+
+
+ ! parameters to be computed based upon parameters above read from file
+ integer NPROC,NEX_XI,NEX_ETA,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,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX
+
+
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
+
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+
+
+ integer :: ifirst_region, ilast_region, iter_region, iter_layer, doubling, tmp_sum, tmp_sum_xi, tmp_sum_eta
+ integer :: NUMBER_OF_MESH_LAYERS,layer_offset,nspec2D_xi_sb,nspec2D_eta_sb, &
+ nb_lay_sb, nspec_sb, nglob_surf
+
+
+ ! for the cut doublingbrick improvement
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
+ logical :: INCLUDE_CENTRAL_CUBE
+ integer :: last_doubling_layer
+ 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 :: tmp_sum_nglob2D_xi, tmp_sum_nglob2D_eta,divider,nglob_edges_h,nglob_edge_v,to_remove
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! calculation of number of elements (NSPEC) below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ ratio_divide_central_cube = maxval(ratio_sampling_array(1:NUMBER_OF_MESH_LAYERS))
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 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
+ tmp_sum_nglob2D_xi = 0
+ tmp_sum_nglob2D_eta = 0
+ do iter_layer = ifirst_region, ilast_region
+ if (this_region_has_a_doubling(iter_layer)) then
+ if (iter_region == IREGION_OUTER_CORE .and. iter_layer == last_doubling_layer) then
+ ! simple brick
+ divider = 1
+ nglob_surf = 6*NGLLX**2 - 7*NGLLX + 2
+ nglob_edges_h = 2*(NGLLX-1)+1 + NGLLX
+ ! minimum value to be safe
+ nglob_edge_v = NGLLX-2
+ nb_lay_sb = 2
+ nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK
+ nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK
+ else
+ ! double brick
+ divider = 2
+ if (ner(iter_layer) == 1) then
+ nglob_surf = 6*NGLLX**2 - 8*NGLLX + 3
+ nglob_edges_h = 4*(NGLLX-1)+1 + 2*(NGLLX-1)+1
+ nglob_edge_v = NGLLX-2
+ nb_lay_sb = 1
+ nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK_1L
+ nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK_1L
+ else
+ nglob_surf = 8*NGLLX**2 - 11*NGLLX + 4
+ nglob_edges_h = 4*(NGLLX-1)+1 + 2*(NGLLX-1)+1
+ nglob_edge_v = 2*(NGLLX-1)+1 -2
+ nb_lay_sb = 2
+ nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK
+ nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK
+ divider = 2
+ endif
+ endif
+ doubling = 1
+ to_remove = 1
+ else
+ if (iter_layer /= ifirst_region) then
+ to_remove = 0
+ else
+ to_remove = 1
+ endif
+ ! dummy values to avoid a warning
+ nglob_surf = 0
+ nglob_edges_h = 0
+ nglob_edge_v = 0
+ divider = 1
+ 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))
+
+ tmp_sum_nglob2D_xi = tmp_sum_nglob2D_xi + (((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)) * &
+ (ner(iter_layer) - doubling*nb_lay_sb))*NGLLX*NGLLX) - &
+ ((((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer))-1)*(ner(iter_layer) - doubling*nb_lay_sb)) + &
+ ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer))*(ner(iter_layer) - to_remove - doubling*nb_lay_sb))*NGLLX) + &
+ (((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer))-1)*(ner(iter_layer) - to_remove - doubling*nb_lay_sb)) + &
+ doubling * (((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer))/divider) * (nglob_surf-nglob_edges_h) - &
+ ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer))/divider -1) * nglob_edge_v)
+
+ tmp_sum_nglob2D_eta = tmp_sum_nglob2D_eta + (((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)) * &
+ (ner(iter_layer) - doubling*nb_lay_sb))*NGLLX*NGLLX) - &
+ ((((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer))-1)*(ner(iter_layer) - doubling*nb_lay_sb)) + &
+ ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer))* &
+ (ner(iter_layer) - to_remove - doubling*nb_lay_sb))*NGLLX) + &
+ (((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer))-1)*(ner(iter_layer) - to_remove - doubling*nb_lay_sb)) + &
+ doubling * (((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer))/divider) * (nglob_surf-nglob_edges_h) - &
+ ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer))/divider -1) * nglob_edge_v)
+
+ enddo ! iter_layer
+
+ NSPEC2D_XI(iter_region) = tmp_sum_xi
+ NSPEC2D_ETA(iter_region) = tmp_sum_eta
+
+ NGLOB2DMAX_YMIN_YMAX(iter_region) = tmp_sum_nglob2D_xi
+ NGLOB2DMAX_XMIN_XMAX(iter_region) = tmp_sum_nglob2D_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))
+
+ NGLOB2DMAX_YMIN_YMAX(iter_region) = NGLOB2DMAX_YMIN_YMAX(iter_region) + &
+ (((NEX_PER_PROC_XI / ratio_divide_central_cube)*(NGLLX-1)+1)*((NEX_XI / ratio_divide_central_cube)*(NGLLX-1)+1))
+
+ NGLOB2DMAX_XMIN_XMAX(iter_region) = NGLOB2DMAX_XMIN_XMAX(iter_region) + &
+ (((NEX_PER_PROC_ETA / ratio_divide_central_cube)*(NGLLX-1)+1)*((NEX_XI / ratio_divide_central_cube)*(NGLLX-1)+1))
+ endif
+ enddo ! iter_region
+
+ ! 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))) / NPROC
+ enddo
+ NSPEC(iter_region) = tmp_sum
+ 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)
+
+ if(minval(NSPEC) <= 0) stop 'negative NSPEC, there is a problem somewhere, try to recompile :) '
+
+ end subroutine count_elements
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_points.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_points.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/count_points.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,242 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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 count_points(NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube,&
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB,&
+ nblocks_xi,nblocks_eta,ner,ratio_sampling_array,&
+ this_region_has_a_doubling,&
+ ifirst_region, ilast_region, iter_region, iter_layer, &
+ doubling, padding, tmp_sum, &
+ INCLUDE_CENTRAL_CUBE,NER_TOP_CENTRAL_CUBE_ICB,NEX_XI, &
+ NUMBER_OF_MESH_LAYERS,layer_offset, &
+ nb_lay_sb, nglob_vol, nglob_surf, nglob_edge, &
+ CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA, &
+ last_doubling_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)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! calculation of number of points (NGLOB) below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ implicit none
+
+ include "constants.h"
+
+! parameters read from parameter file
+
+! parameters to be computed based upon parameters above read from file
+ integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
+
+ integer, dimension(MAX_NUM_REGIONS) :: &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ NGLOB
+
+ integer NER_TOP_CENTRAL_CUBE_ICB,NEX_XI
+ integer nblocks_xi,nblocks_eta
+
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
+
+ integer :: ifirst_region, ilast_region, iter_region, iter_layer, doubling, padding, tmp_sum
+ integer :: NUMBER_OF_MESH_LAYERS,layer_offset, &
+ nb_lay_sb, nglob_vol, nglob_surf, nglob_edge
+
+! for the cut doublingbrick improvement
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,INCLUDE_CENTRAL_CUBE
+ integer :: last_doubling_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
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! 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(:) = NGLOB2DMAX_XMIN_XMAX(:) + 1
+ NGLOB2DMAX_YMIN_YMAX(:) = NGLOB2DMAX_YMIN_YMAX(:) + 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 == last_doubling_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 count_points
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/define_all_layers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/define_all_layers.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/define_all_layers.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,949 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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 define_all_layers(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,&
+ RMIDDLE_CRUST,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RMOHO_FICTITIOUS_IN_MESHER,R80_FICTITIOUS_IN_MESHER,&
+ ONE_CRUST, &
+ ner,ratio_sampling_array,&
+ NUMBER_OF_MESH_LAYERS,layer_offset,last_doubling_layer, &
+ r_bottom,r_top,this_region_has_a_doubling,&
+ ielem,elem_doubling_mantle,elem_doubling_middle_outer_core,&
+ elem_doubling_bottom_outer_core,&
+ DEPTH_SECOND_DOUBLING_REAL,DEPTH_THIRD_DOUBLING_REAL, &
+ DEPTH_FOURTH_DOUBLING_REAL,distance,distance_min,zval,&
+ doubling_index,rmins,rmaxs)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!
+!!!!!! definition of general mesh parameters below
+!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ include "constants.h"
+
+ ! parameters read from parameter file
+ integer 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
+
+ ! radii
+ double precision RMIDDLE_CRUST,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RMOHO_FICTITIOUS_IN_MESHER,R80_FICTITIOUS_IN_MESHER
+
+ logical ONE_CRUST
+
+ ! layers
+ integer :: NUMBER_OF_MESH_LAYERS,layer_offset,last_doubling_layer
+ 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
+
+ ! doubling elements
+ 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
+
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+
+
+! 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
+
+ ! default case:
+ ! no fourth doubling at the bottom of the outer core
+
+ if (SUPPRESS_CRUSTAL_MESH) then
+
+ ! suppress the crustal layers
+ ! will be replaced by an extension of the mantle: R_EARTH is not modified,
+ ! but no more crustal doubling
+
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ 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_FICTITIOUS_IN_MESHER !!!! now fictitious
+
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
+ 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) = R80_FICTITIOUS_IN_MESHER / 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_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / 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
+
+ ! 1D models:
+ ! 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.
+
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ r_top(3) = R80_FICTITIOUS_IN_MESHER
+ 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_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(3) = R80_FICTITIOUS_IN_MESHER / 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
+
+ ! default case for 3D models:
+ ! contains the crustal layers
+ ! doubling at the base of the crust
+
+ 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) ! regional mesh: ner(1) = 1 since NER_CRUST=3
+ ner( 2) = ceiling (NER_CRUST / 2.d0) ! ner(2) = 2
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
+ 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_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / 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
+
+ ! 4th doubling case:
+ ! includes fourth doubling at the bottom of the outer core
+
+ if (SUPPRESS_CRUSTAL_MESH) then
+
+ ! suppress the crustal layers
+ ! will be replaced by an extension of the mantle: R_EARTH is not modified,
+ ! but no more crustal doubling
+
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ 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_FICTITIOUS_IN_MESHER !!!! now fictitious
+
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
+ 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_FICTITIOUS_IN_MESHER / 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_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
+
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / 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
+
+ ! 1D models:
+ ! 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.
+
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ r_top(3) = R80_FICTITIOUS_IN_MESHER
+ 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_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(3) = R80_FICTITIOUS_IN_MESHER / 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
+
+ ! for 3D models:
+ ! contains the crustal layers
+ ! doubling at the base of the crust
+
+ 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.
+ last_doubling_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_FICTITIOUS_IN_MESHER
+
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
+ 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_FICTITIOUS_IN_MESHER / R_EARTH
+
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / 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
+
+
+ end subroutine define_all_layers
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_timestep_and_layers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_timestep_and_layers.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_timestep_and_layers.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,448 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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_timestep_and_layers(DT,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,R_CENTRAL_CUBE, &
+ NEX_MAX,NCHUNKS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,&
+ ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL, &
+ ANISOTROPIC_INNER_CORE)
+
+
+ implicit none
+
+ include "constants.h"
+
+ ! parameters read from parameter file
+ integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD
+
+ integer 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
+
+ integer NEX_MAX,NCHUNKS,REFERENCE_1D_MODEL,THREE_D_MODEL
+
+ double precision DT
+ double precision R_CENTRAL_CUBE
+ double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES
+
+ logical ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL,ANISOTROPIC_INNER_CORE
+
+ ! local variables
+ integer multiplication_factor
+
+ !----
+ !---- 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
+ ! time step
+ DT = 0.252d0
+
+ ! attenuation period range
+ MIN_ATTENUATION_PERIOD = 30
+ MAX_ATTENUATION_PERIOD = 1500
+
+ ! number of element layers in each mesh region
+ 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
+
+ ! radius of central cube
+ 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
+
+ !> Hejun
+ ! avoids elongated elements below the 670-discontinuity,
+ ! since for model REFERENCE_MODEL_1DREF,
+ ! the 670-discontinuity is moved up to 650 km depth.
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1DREF) then
+ NER_771_670 = NER_771_670 + 1
+ end if
+
+ !----
+ !---- 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
+ ! 1D models honor 1D spherical moho
+ if (.not. ONE_CRUST) then
+ ! case 1D + two crustal layers
+ if (NER_CRUST < 2 ) NER_CRUST = 2
+ ! makes time step smaller
+ if(NEX_MAX*multiplication_factor <= 160) then
+ DT = 0.20d0
+ else if(NEX_MAX*multiplication_factor <= 256) then
+ DT = 0.20d0
+ endif
+ endif
+ else
+ ! 3D models: must have two element layers for crust
+ if (NER_CRUST < 2 ) NER_CRUST = 2
+ ! makes time step smaller
+ if(NEX_MAX*multiplication_factor <= 80) then
+ DT = 0.125d0
+ else 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( .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, CRUSTAL, &
+ HONOR_1D_SPHERICAL_MOHO, REFERENCE_1D_MODEL)
+
+ 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)
+
+ !! DK DK suppressed because this routine should not write anything to the screen
+ ! 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(*,*)'##############################################################'
+
+ 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
+
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+
+
+ ! time step reductions are based on empirical values (..somehow)
+
+ ! following models need special attention, at least for global simulations:
+ if( NCHUNKS == 6 ) then
+ ! makes time step smaller for this ref model, otherwise becomes unstable in fluid
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91) &
+ DT = DT*(1.d0 - 0.3d0)
+
+ ! using inner core anisotropy, simulations might become unstable in solid
+ if( ANISOTROPIC_INNER_CORE ) then
+ ! DT = DT*(1.d0 - 0.1d0) not working yet...
+ stop 'anisotropic inner core - unstable feature, uncomment this line in get_timestep_and_layers.f90'
+ endif
+ endif
+
+ ! following models need special attention, regardless of number of chunks:
+ ! it makes the time step smaller for this ref model, otherwise becomes unstable in fluid
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) &
+ DT = DT*(1.d0 - 0.8d0) ! *0.20d0
+
+
+ if( ITYPE_CRUSTAL_MODEL == ICRUST_CRUSTMAPS ) &
+ DT = DT*(1.d0 - 0.3d0)
+
+ ! decreases time step as otherwise the solution might become unstable for rougher/unsmoothed models
+ if( .false. ) then
+ if( THREE_D_MODEL == THREE_D_MODEL_PPM ) DT = DT * (1.d0 - 0.2d0)
+ endif
+
+ ! takes a 5% safety margin on the maximum stable time step
+ ! which was obtained by trial and error
+ DT = DT * (1.d0 - 0.05d0)
+
+ ! adapts number of element layers in crust and time step for regional simulations
+ if( REGIONAL_MOHO_MESH ) then
+ ! hard coded number of crustal element layers and time step
+
+ ! checks
+ if( NCHUNKS > 1 ) stop 'regional moho mesh: NCHUNKS error in rcp_set_timestep_and_layers'
+ if( HONOR_1D_SPHERICAL_MOHO ) return
+
+ ! original values
+ !print*,'NER:',NER_CRUST
+ !print*,'DT:',DT
+
+ ! enforce 3 element layers
+ NER_CRUST = 3
+
+ ! increased stability, empirical
+ DT = DT*(1.d0 + 0.5d0)
+
+ if( REGIONAL_MOHO_MESH_EUROPE ) DT = 0.17 ! europe
+ if( REGIONAL_MOHO_MESH_ASIA ) DT = 0.15 ! asia & middle east
+
+ endif
+
+ end subroutine get_timestep_and_layers
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90 2013-07-01 01:22:13 UTC (rev 22467)
@@ -0,0 +1,608 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! 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 write_VTK_data_points(nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ points_globalindices,num_points_globalindices, &
+ prname_file)
+
+! external mesh routine for saving vtk files for points locations
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nglob
+
+ ! global coordinates
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! gll data values array
+ integer :: num_points_globalindices
+ integer, dimension(num_points_globalindices) :: points_globalindices
+
+ ! file name
+ character(len=150) prname_file
+
+ integer :: i,iglob
+
+ ! write source and receiver VTK files for Paraview
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', num_points_globalindices, ' float'
+ do i=1,num_points_globalindices
+ iglob = points_globalindices(i)
+ if( iglob <= 0 .or. iglob > nglob ) then
+ print*,'error: '//prname_file(1:len_trim(prname_file))//'.vtk'
+ print*,'error global index: ',iglob,i
+ close(IOVTK)
+ stop 'error vtk points file'
+ endif
+
+ write(IOVTK,'(3e18.6)') xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
+ enddo
+ write(IOVTK,*) ""
+
+ close(IOVTK)
+
+ end subroutine write_VTK_data_points
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine write_VTK_glob_points(nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ glob_values, &
+ prname_file)
+
+! external mesh routine for saving vtk files for points locations
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nglob
+
+ ! global coordinates
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! gll data values array
+ real(kind=CUSTOM_REAL), dimension(nglob) :: glob_values
+
+ ! file name
+ character(len=150) prname_file
+
+ integer :: iglob
+
+ ! write source and receiver VTK files for Paraview
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ do iglob=1,nglob
+ write(IOVTK,*) xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
+ enddo
+ write(IOVTK,*) ""
+
+ ! writes out gll-data (velocity) for each element point
+ write(IOVTK,'(a,i12)') "POINT_DATA ",nglob
+ write(IOVTK,'(a)') "SCALARS glob_data float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do iglob=1,nglob
+ write(IOVTK,*) glob_values(iglob)
+ enddo
+ write(IOVTK,*) ""
+
+ close(IOVTK)
+
+ end subroutine write_VTK_glob_points
+
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine write_VTK_data_elem_l(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ elem_flag,prname_file)
+
+! routine for saving vtk file holding logical flag on each spectral element
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nspec,nglob
+
+ ! global coordinates
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! element flag array
+ logical, dimension(nspec) :: elem_flag
+ integer :: ispec,i
+
+ ! file name
+ character(len=150) prname_file
+
+ ! write source and receiver VTK files for Paraview
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ do i=1,nglob
+ write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ enddo
+ write(IOVTK,*) ""
+
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ do ispec=1,nspec
+ write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
+ ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1,ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
+ enddo
+ write(IOVTK,*) ""
+
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOVTK,*) (12,ispec=1,nspec)
+ write(IOVTK,*) ""
+
+ write(IOVTK,'(a,i12)') "CELL_DATA ",nspec
+ write(IOVTK,'(a)') "SCALARS elem_flag integer"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do ispec = 1,nspec
+ if( elem_flag(ispec) .eqv. .true. ) then
+ write(IOVTK,*) 1
+ else
+ write(IOVTK,*) 0
+ endif
+ enddo
+ write(IOVTK,*) ""
+ close(IOVTK)
+
+
+ end subroutine write_VTK_data_elem_l
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine write_VTK_data_elem_i(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ elem_flag,prname_file)
+
+
+! routine for saving vtk file holding integer value on each spectral element
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nspec,nglob
+
+ ! global coordinates
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! element flag array
+ integer, dimension(nspec) :: elem_flag
+ integer :: ispec,i
+
+ ! file name
+ character(len=150) prname_file
+
+ ! write source and receiver VTK files for Paraview
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ do i=1,nglob
+ write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ enddo
+ write(IOVTK,*) ""
+
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ do ispec=1,nspec
+ write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
+ ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1,ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
+ enddo
+ write(IOVTK,*) ""
+
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOVTK,*) (12,ispec=1,nspec)
+ write(IOVTK,*) ""
+
+ write(IOVTK,'(a,i12)') "CELL_DATA ",nspec
+ write(IOVTK,'(a)') "SCALARS elem_val integer"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do ispec = 1,nspec
+ write(IOVTK,*) elem_flag(ispec)
+ enddo
+ write(IOVTK,*) ""
+ close(IOVTK)
+
+ end subroutine write_VTK_data_elem_i
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! external mesh routine for saving vtk files for custom_real values on global points
+
+ subroutine write_VTK_data_cr(idoubling,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ glob_data,prname_file)
+
+! outputs single file for each process
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nspec,nglob
+
+ integer, dimension(nspec):: idoubling
+
+ ! global coordinates
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! global data values array
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: glob_data
+
+ ! file name
+ character(len=256) prname_file
+
+ ! local parameters
+ integer :: ispec,i
+ real(kind=CUSTOM_REAL) :: rval,thetaval,phival,xval,yval,zval
+
+ ! write source and receiver VTK files for Paraview
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ do i=1,nglob
+
+ !x,y,z store have been converted to r theta phi already, need to revert back for xyz output
+ rval = xstore_dummy(i)
+ thetaval = ystore_dummy(i)
+ phival = zstore_dummy(i)
+ call rthetaphi_2_xyz(xval,yval,zval,rval,thetaval,phival)
+
+ !write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOVTK,'(3e18.6)') xval,yval,zval
+ enddo
+ write(IOVTK,*) ""
+
+ ! defines cell on coarse corner points
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ do ispec=1,nspec
+
+ ! specific to inner core elements
+ ! exclude fictitious elements in central cube
+ if(idoubling(ispec) /= IFLAG_IN_FICTITIOUS_CUBE) then
+ ! valid cell
+ write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1, &
+ ibool(NGLLX,1,1,ispec)-1, &
+ ibool(NGLLX,NGLLY,1,ispec)-1, &
+ ibool(1,NGLLY,1,ispec)-1, &
+ ibool(1,1,NGLLZ,ispec)-1, &
+ ibool(NGLLX,1,NGLLZ,ispec)-1, &
+ ibool(NGLLX,NGLLY,NGLLZ,ispec)-1, &
+ ibool(1,NGLLY,NGLLZ,ispec)-1
+ else
+ ! fictitious elements in central cube
+ ! maps cell onto a randomly chosen point
+ write(IOVTK,'(9i12)') 8,ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1, &
+ ibool(1,1,1,1)-1
+ endif
+
+ enddo
+ write(IOVTK,*) ""
+
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOVTK,*) (12,ispec=1,nspec)
+ write(IOVTK,*) ""
+
+ ! x components
+ write(IOVTK,'(a,i12)') "POINT_DATA ",nglob
+ write(IOVTK,'(a)') "SCALARS x_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do i = 1,nglob
+ write(IOVTK,*) glob_data(1,i)
+ enddo
+ ! y components
+ write(IOVTK,'(a)') "SCALARS y_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do i = 1,nglob
+ write(IOVTK,*) glob_data(2,i)
+ enddo
+ ! z components
+ write(IOVTK,'(a)') "SCALARS z_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do i = 1,nglob
+ write(IOVTK,*) glob_data(3,i)
+ enddo
+ ! norm
+ write(IOVTK,'(a)') "SCALARS norm float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do i = 1,nglob
+ write(IOVTK,*) sqrt( glob_data(1,i)*glob_data(1,i) &
+ + glob_data(2,i)*glob_data(2,i) &
+ + glob_data(3,i)*glob_data(3,i))
+ enddo
+ write(IOVTK,*) ""
+
+ close(IOVTK)
+
+
+ end subroutine write_VTK_data_cr
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! external mesh routine for saving vtk files for custom_real values on global points
+
+ subroutine write_VTK_data_cr_all(myrank,NPROCTOT,idoubling, &
+ nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ glob_data,prname_file)
+
+! outputs single file for all processes
+
+ implicit none
+
+ include "constants.h"
+
+ include 'mpif.h'
+ include "precision.h"
+
+ integer :: myrank,NPROCTOT
+
+ integer ::nspec,nglob
+
+ integer, dimension(nspec):: idoubling
+
+ ! global coordinates
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! global data values array
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: glob_data
+
+ ! file name
+ character(len=256) prname_file
+
+ ! local parameters
+ integer :: ispec,i,iproc,ier
+ real(kind=CUSTOM_REAL) :: rval,thetaval,phival,xval,yval,zval
+
+ real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: &
+ store_val_x_all,store_val_y_all,store_val_z_all, &
+ store_val_ux_all,store_val_uy_all,store_val_uz_all
+ integer, dimension(:,:,:,:,:),allocatable :: ibool_all
+ integer, dimension(:,:),allocatable :: idoubling_all
+
+ ! master collect arrays
+ if( myrank == 0 ) then
+ allocate(store_val_x_all(nglob,0:NPROCTOT-1), &
+ store_val_y_all(nglob,0:NPROCTOT-1), &
+ store_val_z_all(nglob,0:NPROCTOT-1), &
+ store_val_ux_all(nglob,0:NPROCTOT-1), &
+ store_val_uy_all(nglob,0:NPROCTOT-1), &
+ store_val_uz_all(nglob,0:NPROCTOT-1), &
+ idoubling_all(nspec,0:NPROCTOT-1), &
+ ibool_all(NGLLX,NGLLY,NGLLZ,nspec,0:NPROCTOT-1),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating stores')
+ else
+ ! dummy arrays
+ allocate(store_val_x_all(1,1), &
+ store_val_y_all(1,1), &
+ store_val_z_all(1,1), &
+ store_val_ux_all(1,1), &
+ store_val_uy_all(1,1), &
+ store_val_uz_all(1,1), &
+ idoubling_all(1,1), &
+ ibool_all(1,1,1,1,1),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy stores')
+ endif
+
+ ! gather info on master proc
+ call MPI_GATHER(xstore_dummy,nglob,CUSTOM_MPI_TYPE,store_val_x_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(ystore_dummy,nglob,CUSTOM_MPI_TYPE,store_val_y_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(zstore_dummy,nglob,CUSTOM_MPI_TYPE,store_val_z_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+ call MPI_GATHER(glob_data(1,:),nglob,CUSTOM_MPI_TYPE,store_val_ux_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(glob_data(2,:),nglob,CUSTOM_MPI_TYPE,store_val_uy_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(glob_data(3,:),nglob,CUSTOM_MPI_TYPE,store_val_uz_all,nglob,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+ call MPI_GATHER(ibool,NGLLX*NGLLY*NGLLZ*nspec,MPI_INTEGER,ibool_all, &
+ NGLLX*NGLLY*NGLLZ*nspec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(idoubling,nspec,MPI_INTEGER,idoubling_all,nspec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+
+
+ if( myrank == 0 ) then
+
+ ! write source and receiver VTK files for Paraview
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob*NPROCTOT, ' float'
+ do iproc=0, NPROCTOT-1
+ do i=1,nglob
+
+ !x,y,z store have been converted to r theta phi already, need to revert back for xyz output
+ rval = store_val_x_all(i,iproc)
+ thetaval = store_val_y_all(i,iproc)
+ phival = store_val_z_all(i,iproc)
+ call rthetaphi_2_xyz(xval,yval,zval,rval,thetaval,phival)
+
+ !write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOVTK,'(3e18.6)') xval,yval,zval
+ enddo
+ enddo
+ write(IOVTK,*) ""
+
+ ! defines cell on coarse corner points
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nspec*NPROCTOT,nspec*NPROCTOT*9
+ do iproc=0, NPROCTOT-1
+ do ispec=1,nspec
+
+ ! note: central cube elements are only shared and used in CHUNK_AB and CHUNK_AB_ANTIPODE
+ ! all other chunks ignore those elements
+
+ ! specific to inner core elements
+ ! exclude fictitious elements in central cube
+ if(idoubling_all(ispec,iproc) /= IFLAG_IN_FICTITIOUS_CUBE) then
+ ! valid cell
+ ! cell corner ids
+ write(IOVTK,'(9i12)') 8,ibool_all(1,1,1,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(NGLLX,1,1,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(NGLLX,NGLLY,1,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(1,NGLLY,1,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(1,1,NGLLZ,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(NGLLX,1,NGLLZ,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(NGLLX,NGLLY,NGLLZ,ispec,iproc)-1+iproc*nglob, &
+ ibool_all(1,NGLLY,NGLLZ,ispec,iproc)-1+iproc*nglob
+ else
+ ! fictitious elements in central cube
+ ! maps cell onto a randomly chosen point
+ write(IOVTK,'(9i12)') 8,ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1, &
+ ibool_all(1,1,1,1,iproc)-1
+ endif
+
+ enddo
+ enddo
+ write(IOVTK,*) ""
+
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec*NPROCTOT
+ write(IOVTK,*) (12,ispec=1,nspec*NPROCTOT)
+ write(IOVTK,*) ""
+
+ ! x components
+ write(IOVTK,'(a,i12)') "POINT_DATA ",nglob*NPROCTOT
+ write(IOVTK,'(a)') "SCALARS x_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do iproc=0, NPROCTOT-1
+ do i = 1,nglob
+ write(IOVTK,*) store_val_ux_all(i,iproc)
+ enddo
+ enddo
+ ! y components
+ write(IOVTK,'(a)') "SCALARS y_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do iproc=0, NPROCTOT-1
+ do i = 1,nglob
+ write(IOVTK,*) store_val_uy_all(i,iproc)
+ enddo
+ enddo
+ ! z components
+ write(IOVTK,'(a)') "SCALARS z_comp float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do iproc=0, NPROCTOT-1
+ do i = 1,nglob
+ write(IOVTK,*) store_val_uz_all(i,iproc)
+ enddo
+ enddo
+ ! norm
+ write(IOVTK,'(a)') "SCALARS norm float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do iproc=0, NPROCTOT-1
+ do i = 1,nglob
+ write(IOVTK,*) sqrt( store_val_ux_all(i,iproc)**2 &
+ + store_val_uy_all(i,iproc)**2 &
+ + store_val_uz_all(i,iproc)**2 )
+ enddo
+ enddo
+ write(IOVTK,*) ""
+
+ close(IOVTK)
+
+ endif
+
+ deallocate(store_val_x_all,store_val_y_all,store_val_z_all, &
+ store_val_ux_all,store_val_uy_all,store_val_uz_all, &
+ ibool_all)
+
+ end subroutine write_VTK_data_cr_all
More information about the CIG-COMMITS
mailing list