[cig-commits] r11590 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . UTILS UTILS/doubling_brick UTILS/estimate_best_values_runs
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Mar 26 09:02:54 PDT 2008
Author: dkomati1
Date: 2008-03-26 09:02:54 -0700 (Wed, 26 Mar 2008)
New Revision: 11590
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/DATA
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/estimate_best_values_runs.f90
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/doubling_brick/generate_table_31_Manual.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/generate_table_NEX_manual.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/jp3d1994_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp1d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea1d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/sea99_s_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
- finished routine to compute best runs in terms of memory use for a given machine
- renamed "compteur" to "counter"
- renamed "compteur" to "counter"
- removed useless white spaces
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/doubling_brick/generate_table_31_Manual.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/doubling_brick/generate_table_31_Manual.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/doubling_brick/generate_table_31_Manual.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -1,35 +0,0 @@
-
-program generate
-
-! generate Table 3.1 of the SPECFEM3D_GLOBE manual
-
-implicit none
-
-integer NPROC_XI,c,NEX_XI,compteur
-
-integer, parameter :: NB_COLONNES = 7
-integer, dimension(NB_COLONNES) :: store_NEX_XI
-
-! base value depends if we implement three or four doublings
-!!!integer, parameter :: BASE_VALUE = 16
-integer, parameter :: BASE_VALUE = 32
-
-do NPROC_XI = 1,26
-
- compteur = 1
-
- do c = 1,20
- NEX_XI = BASE_VALUE * c * NPROC_XI
- if(NEX_XI >= 64 .and. compteur <= NB_COLONNES) then
- store_NEX_XI(compteur) = NEX_XI
- compteur = compteur + 1
- endif
- enddo
-
- write(*,"(i6,' | ',i6,' | ',i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
- NPROC_XI, 6*NPROC_XI**2,(store_NEX_XI(compteur), compteur = 1,NB_COLONNES)
-
-enddo
-
-end
-
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/DATA
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/DATA (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/DATA 2008-03-26 16:02:54 UTC (rev 11590)
@@ -0,0 +1 @@
+link ../../DATA/
\ No newline at end of file
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/DATA
___________________________________________________________________
Name: svn:special
+ *
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -1,273 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory, California Institute of Technology, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau / CNRS / INRIA
-! February 2008
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-! create file OUTPUT_FILES/values_from_mesher.h based upon DATA/Par_file
-! in order to compile the solver with the right array sizes
-
- subroutine create_header_file
-
- implicit none
-
- include "constants.h"
-
-! parameters read from parameter file
- integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
- NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
- NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
- NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NSOURCES,NTSTEP_BETWEEN_FRAMES, &
- NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
-
- double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
- CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
- RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
- R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
- MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH
-
- logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
- TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D, &
- RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
- SAVE_MESH_FILES,ATTENUATION,CASE_3D, &
- ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
- OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
- ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,first_time
-
- character(len=150) LOCAL_PATH,MODEL
-
-! parameters deduced from parameters read from file
- integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
- integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
- integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
- logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
- double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
- double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
-
-! this for all the regions
- integer, dimension(MAX_NUM_REGIONS) :: NSPEC, &
- NSPEC2D_XI, &
- NSPEC2D_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
- NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
- nglob
-
- double precision :: static_memory_size
- character(len=150) HEADER_FILE
-
- integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
- NSPEC_INNER_CORE_ATTENUATION, &
- NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
- NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
- NSPEC_CRUST_MANTLE_ADJOINT, &
- NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
- NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
- NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
- NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION
-
- integer :: iregion
- logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
- integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
- integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
- integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NGLOB1D_RADIAL_CORNER
- integer, dimension(MAX_NUM_REGIONS) :: NGLOB1D_RADIAL_TEMP
-
- integer :: c,compteur
-
- double precision :: mem_per_core,percent
-
- integer, parameter :: NB_COLONNES = 10
- integer, dimension(NB_COLONNES) :: store_NEX_XI
-
-! maximum total number of processors we want to see in the table
- integer, parameter :: MAX_NUMBER_OF_PROCS = 62976 !! Ranger !! 100000
-! integer, parameter :: MAX_NUMBER_OF_PROCS = 212992 ! current maximum on BlueGene at LLNL
-
-! amount of memory installed per processor core on the system (in gigabytes)
- double precision, parameter :: MAX_MEMORY_PER_CORE = 2.d0 !! Ranger
-
-! base value depends if we implement three or four doublings (default is three)
- integer, parameter :: NB_DOUBLING = 3
- integer :: BASE_VALUE
-
-! output in LaTeX format or in regular ASCII format
- logical, parameter :: OUTPUT_LATEX_FORMAT = .false.
-
-! ************** PROGRAM STARTS HERE **************
-
- call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
-
-! count the total number of sources in the CMTSOLUTION file
- NSOURCES = 1
-
- BASE_VALUE = 2**NB_DOUBLING
-
- first_time = .true.
-
- print *
- print *,'Number of GB of memory per core on the machine: ',MAX_MEMORY_PER_CORE
- print *
-
- print *,'total_proc, % of Ranger, NPROC_XI, NEX_XI, memory used per core, percentage:'
- print *
-
-! total number of processors is 6 * NPROC_XI^2
-!!!! do NPROC_XI = 1,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
-!! start after 1000 processors
- do NPROC_XI = 13,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
-
- compteur = 1
- c = 0
-
- do while (compteur <= NB_COLONNES)
- c = c + 1
- NEX_XI = BASE_VALUE * c * NPROC_XI
-
- if(NEX_XI >= 64 .and. mod(NEX_XI,2*BASE_VALUE) == 0) then
-
- store_NEX_XI(compteur) = NEX_XI
- compteur = compteur + 1
-
-! read the parameter file and compute additional parameters
- call read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
- NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
- NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
- NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
- ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
- CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
- RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
- R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
- MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
- TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
- ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
- ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
- MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
- PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
- ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
- NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NSPEC, &
- NSPEC2D_XI, &
- NSPEC2D_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB, &
- ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
- OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
- ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
- DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
-
- do iregion=1,MAX_NUM_REGIONS
- NGLOB1D_RADIAL_CORNER(iregion,:) = NGLOB1D_RADIAL(iregion)
- enddo
-
- if (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA) then
- NGLOB1D_RADIAL_CORNER(IREGION_OUTER_CORE,:) = NGLOB1D_RADIAL_CORNER(IREGION_OUTER_CORE,:) + &
- maxval(DIFF_NSPEC1D_RADIAL(:,:))*(NGLLZ-1)
- endif
-
- if(first_time) then
- first_time = .false.
- if(ATTENUATION) then
- print *,'ATTENUATION = .true.'
- else
- print *,'ATTENUATION = .false.'
- endif
- print *
- endif
-
-!! compute memory usage per processor core
-! evaluate the amount of static memory needed by the solver
- call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
- TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
- ONE_CRUST,doubling_index,this_region_has_a_doubling,&
- ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
- NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
- NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
- NSPEC_INNER_CORE_ATTENUATION, &
- NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
- NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
- NSPEC_CRUST_MANTLE_ADJOINT, &
- NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
- NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
- NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
- NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
-
- mem_per_core = static_memory_size/1073741824.d0
-
- percent = 100.d0 * mem_per_core / MAX_MEMORY_PER_CORE
-
- if(percent > 100.d0) goto 777
-
- if(percent < 0.d0) then
- goto 777
-! write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **mesher fails/temporary bug**')") &
-! 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
-
- else if(percent >= 93.d0) then
- write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **too high**')") &
- 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
- goto 777
-
- else if(percent >= 85.d0) then
- write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **excellent**')") &
- 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
-
- else if(percent >= 80.d0) then
- write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **very good**')") &
- 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
-
- else
- write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'%')") &
- 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
- endif
-
- endif
- enddo
-
- 777 continue
-
- print *
-
- enddo
-
- end subroutine create_header_file
-
Copied: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/estimate_best_values_runs.f90 (from rev 11585, seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/create_header_file.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/estimate_best_values_runs.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/estimate_best_values_runs.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -0,0 +1,278 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! output a table with the best parameters to use for large runs on a given machine
+! in order to use the total amount of memory available in the best possible way
+
+ program estimate_best_values_runs
+
+ implicit none
+
+ include "../../constants.h"
+
+ integer, parameter :: NB_COLUMNS_TABLE = 10
+
+! maximum total number of processors we want to see in the table
+! integer, parameter :: MAX_NUMBER_OF_PROCS = 100000
+! integer, parameter :: MAX_NUMBER_OF_PROCS = 4000 ! current maximum on pangu at Caltech
+! integer, parameter :: MAX_NUMBER_OF_PROCS = 10240 ! current maximum on MareNostrum in Barcelona
+ integer, parameter :: MAX_NUMBER_OF_PROCS = 62976 ! current maximum on Ranger in Texas
+! integer, parameter :: MAX_NUMBER_OF_PROCS = 212992 ! current maximum on BlueGene at LLNL
+
+! minimum total number of processors we want to see in the table
+! integer, parameter :: MIN_NUMBER_OF_PROCS = 100
+ integer, parameter :: MIN_NUMBER_OF_PROCS = 1000
+
+! amount of memory installed per processor core on the system (in gigabytes)
+! double precision, parameter :: MAX_MEMORY_PER_CORE = 1.5d0 ! pangu at Caltech
+ double precision, parameter :: MAX_MEMORY_PER_CORE = 2.0d0 ! Ranger in Texas, MareNostrum in Barcelona
+
+! base value depends if we implement three or four doublings (default is three)
+ integer, parameter :: NB_DOUBLING = 3
+ integer :: BASE_VALUE
+
+! parameters read from parameter file
+ integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NSOURCES,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+
+ double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH
+
+ logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
+ TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D, &
+ RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
+ SAVE_MESH_FILES,ATTENUATION,CASE_3D, &
+ ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,first_time
+
+ character(len=150) LOCAL_PATH,MODEL
+
+! parameters deduced from parameters read from file
+ integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+
+! this for all the regions
+ integer, dimension(MAX_NUM_REGIONS) :: NSPEC,NSPEC2D_XI,NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,nglob
+
+ double precision :: static_memory_size
+ integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION
+
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,EMULATE_ONLY
+ 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 :: c,counter,offset
+
+ double precision :: mem_per_core,percent
+
+! ************** PROGRAM STARTS HERE **************
+
+! count the total number of sources in the CMTSOLUTION file
+ NSOURCES = 1
+
+ BASE_VALUE = 2**NB_DOUBLING
+
+ first_time = .true.
+
+ print *
+ print *,'Number of GB of memory per core installed in the machine: ',MAX_MEMORY_PER_CORE
+ print *
+
+ print *,'NPROC, % of the machine, NPROC_XI, NEX_XI, memory used per core, percentage:'
+ print *
+
+! make sure we do not start below the lower limit
+ if(6 * int(sqrt(MIN_NUMBER_OF_PROCS / 6.d0))**2 == MIN_NUMBER_OF_PROCS) then
+ offset = 0
+ else
+ offset = 1
+ endif
+
+! the total number of processors is 6 * NPROC_XI^2
+ do NPROC_XI = int(sqrt(MIN_NUMBER_OF_PROCS / 6.d0)) + offset,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
+
+ counter = 1
+ c = 0
+
+ do while (counter <= NB_COLUMNS_TABLE)
+ c = c + 1
+ NEX_XI = BASE_VALUE * c * NPROC_XI
+
+ if(NEX_XI >= 64 .and. mod(NEX_XI,2*BASE_VALUE) == 0) then
+
+ counter = counter + 1
+
+! read the parameter file and compute additional parameters
+ EMULATE_ONLY = .true.
+ call read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
+ ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
+ MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
+ PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
+ ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB, &
+ ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,EMULATE_ONLY)
+
+ if(first_time) then
+ first_time = .false.
+ if(ATTENUATION) then
+ print *,'ATTENUATION = .true.'
+ else
+ print *,'ATTENUATION = .false.'
+ endif
+ print *
+ endif
+
+!! compute memory usage per processor core
+! evaluate the amount of static memory needed by the solver
+ call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
+ ONE_CRUST,doubling_index,this_region_has_a_doubling,&
+ ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
+ NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
+ NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+
+ mem_per_core = static_memory_size/1073741824.d0
+
+ percent = 100.d0 * mem_per_core / MAX_MEMORY_PER_CORE
+
+ if(percent > 100.d0) goto 777
+
+ if(percent < 0.d0) then
+ goto 777
+! write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **mesher fails/temporary bug**')") &
+! 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+
+ else if(percent >= 93.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **too high**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+ goto 777
+
+ else if(percent >= 85.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **excellent**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+ goto 777
+
+ else if(percent >= 80.d0) then
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'% **very good**')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+
+ else
+ write(*,"(' ',i5,' ',f6.2,'% ',i3,' ',i4,' ',f6.2,' ',f6.2,'%')") &
+ 6*NPROC_XI**2,100.d0*6*NPROC_XI**2/dble(MAX_NUMBER_OF_PROCS),NPROC_XI,NEX_XI,mem_per_core,percent
+ endif
+
+ endif
+ enddo
+
+ 777 continue
+
+ print *
+
+ enddo
+
+ end program estimate_best_values_runs
+
+!
+!---- include subroutines from the main code
+!
+
+ include "../../read_compute_parameters.f90"
+
+ include "../../memory_eval.f90"
+
+ include "../../read_value_parameters.f90"
+
+ include "../../get_value_parameters.f90"
+
+ include "../../auto_ner.f90"
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/estimate_best_values_runs/read_compute_parameters.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -1,2338 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory, California Institute of Technology, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau / CNRS / INRIA
-! February 2008
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
- NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
- NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
- NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
- ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
- CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
- RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
- R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
- MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
- TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
- ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
- ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
- MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
- PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
- ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
- NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NSPEC, &
- NSPEC2D_XI, &
- NSPEC2D_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
- NGLOB, &
- ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
- OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
- ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
- DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
-
-
- implicit none
-
- include "constants.h"
-
- integer NEX_ETA_dummy , NEX_XI_dummy, NPROC_ETA_dummy , NPROC_XI_dummy
-
-! parameters read from parameter file
- integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
- NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
- NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
- NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
-
- double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
- CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
- RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
- R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
- MOVIE_TOP_KM,MOVIE_TOP,MOVIE_BOTTOM_KM,MOVIE_BOTTOM, &
- MOVIE_EAST_DEG,MOVIE_EAST,MOVIE_WEST_DEG,MOVIE_WEST,MOVIE_NORTH_DEG,MOVIE_NORTH,MOVIE_SOUTH_DEG,MOVIE_SOUTH
-
- logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
- TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_VOLUME_COARSE,ATTENUATION_3D, &
- RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
- SAVE_MESH_FILES,ATTENUATION, &
- ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
- OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
- ROTATE_SEISMOGRAMS_RT,WRITE_SEISMOGRAMS_BY_MASTER,&
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE
-
- character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
-
-! local variables
- integer NEX_MAX
-
- double precision RECORD_LENGTH_IN_MINUTES,ELEMENT_WIDTH
-
- integer, external :: err_occurred
-
-! parameters to be computed based upon parameters above read from file
- integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
-
- integer, dimension(MAX_NUM_REGIONS) :: NSPEC, &
- NSPEC2D_XI, &
- NSPEC2D_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
- NGLOB
-
- integer nblocks_xi,nblocks_eta
-
- integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
- double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
- logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
-
- integer :: ielem,elem_doubling_mantle,elem_doubling_middle_outer_core,elem_doubling_bottom_outer_core
- double precision :: DEPTH_SECOND_DOUBLING_REAL,DEPTH_THIRD_DOUBLING_REAL, &
- DEPTH_FOURTH_DOUBLING_REAL,distance,distance_min,zval
-
-! honor PREM Moho or not
-! doing so drastically reduces the stability condition and therefore the time step
- logical :: HONOR_1D_SPHERICAL_MOHO,CASE_3D
-
- integer :: ifirst_region, ilast_region, iter_region, iter_layer, doubling, padding, tmp_sum, tmp_sum_xi, tmp_sum_eta
- integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
- double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
- integer :: NUMBER_OF_MESH_LAYERS,layer_offset,nspec2D_xi_sb,nspec2D_eta_sb, &
- nb_lay_sb, nspec_sb, nglob_vol, nglob_surf, nglob_edge
-
- integer :: multiplication_factor
-
-! for the cut doublingbrick improvement
- logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
- integer :: lastdoubling_layer, cut_doubling, nglob_int_surf_xi, nglob_int_surf_eta,nglob_ext_surf,&
- normal_doubling, nglob_center_edge, nglob_corner_edge, nglob_border_edge
- integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
- integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
-
-! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
- call open_parameter_file
-
- call read_value_integer(SIMULATION_TYPE, 'solver.SIMULATION_TYPE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(SAVE_FORWARD, 'solver.SAVE_FORWARD')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
- call read_value_integer(NCHUNKS, 'mesher.NCHUNKS')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
- stop 'NCHUNKS must be either 1, 2, 3 or 6'
-
- call read_value_double_precision(ANGULAR_WIDTH_XI_IN_DEGREES, 'mesher.ANGULAR_WIDTH_XI_IN_DEGREES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(ANGULAR_WIDTH_ETA_IN_DEGREES, 'mesher.ANGULAR_WIDTH_ETA_IN_DEGREES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(CENTER_LATITUDE_IN_DEGREES, 'mesher.CENTER_LATITUDE_IN_DEGREES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(CENTER_LONGITUDE_IN_DEGREES, 'mesher.CENTER_LONGITUDE_IN_DEGREES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(GAMMA_ROTATION_AZIMUTH, 'mesher.GAMMA_ROTATION_AZIMUTH')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
-! this MUST be 90 degrees for two chunks or more to match geometrically
- if(NCHUNKS > 1 .and. abs(ANGULAR_WIDTH_XI_IN_DEGREES - 90.d0) > 0.00000001d0) &
- stop 'ANGULAR_WIDTH_XI_IN_DEGREES must be 90 for more than one chunk'
-
-! this can be any value in the case of two chunks
- if(NCHUNKS > 2 .and. abs(ANGULAR_WIDTH_ETA_IN_DEGREES - 90.d0) > 0.00000001d0) &
- stop 'ANGULAR_WIDTH_ETA_IN_DEGREES must be 90 for more than two chunks'
-
-! include central cube or not
-! use regular cubed sphere instead of cube for large distances
- if(NCHUNKS == 6) then
- INCLUDE_CENTRAL_CUBE = .true.
- INFLATE_CENTRAL_CUBE = .false.
- else
- INCLUDE_CENTRAL_CUBE = .false.
- INFLATE_CENTRAL_CUBE = .true.
- endif
-
-! number of elements at the surface along the two sides of the first chunk
- call read_value_integer(NEX_XI_dummy, 'mesher.NEX_XI')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NEX_ETA_dummy, 'mesher.NEX_ETA')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NPROC_XI_dummy, 'mesher.NPROC_XI')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NPROC_ETA_dummy, 'mesher.NPROC_ETA')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-!! DK DK
- NEX_ETA = NEX_XI
- NPROC_ETA = NPROC_XI
-
-! define the velocity model
- call read_value_string(MODEL, 'model.name')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
-! use PREM as the 1D reference model by default
- REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
-
-! HONOR_1D_SPHERICAL_MOHO: honor PREM Moho or not: doing so drastically reduces
-! the stability condition and therefore the time step, resulting in expensive
-! calculations. If not, honor a fictitious Moho at the depth of 40 km
-! in order to have even radial sampling from the d220 to the Earth surface.
-
-! ONE_CRUST: in order to increase stability and therefore to allow cheaper
-! simulations (larger time step), 1D models can be run with just one average crustal
-! layer instead of two.
-
-! CASE_3D : this flag allows the stretching of the elements in the crustal
-! layers in the case of 3D models. The purpose of this stretching is to squeeze more
-! GLL points per km in the upper part of the crust than in the lower part.
- HONOR_1D_SPHERICAL_MOHO = .false.
- ONE_CRUST = .false.
- CASE_3D = .false.
-
-! default is no 3D model
- THREE_D_MODEL = 0
-
- if(MODEL == '1D_isotropic_prem') then
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
-
- else if(MODEL == '1D_transversely_isotropic_prem') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
-
- else if(MODEL == '1D_iasp91' .or. MODEL == '1D_1066a' .or. MODEL == '1D_ak135') then
- if(MODEL == '1D_iasp91') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_IASP91
- else if(MODEL == '1D_1066a') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_1066A
- else if(MODEL == '1D_ak135') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
- else
- stop 'reference 1D Earth model unknown'
- endif
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
-
- else if(MODEL == '1D_ref') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
-
- else if(MODEL == '1D_ref_iso') then
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
-
- else if(MODEL == '1D_isotropic_prem_onecrust') then
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
- ONE_CRUST = .true.
-
- else if(MODEL == '1D_transversely_isotropic_prem_onecrust') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
- ONE_CRUST = .true.
-
- else if(MODEL == '1D_iasp91_onecrust' .or. MODEL == '1D_1066a_onecrust' .or. MODEL == '1D_ak135_onecrust') then
- if(MODEL == '1D_iasp91_onecrust') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_IASP91
- else if(MODEL == '1D_1066a_onecrust') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_1066A
- else if(MODEL == '1D_ak135_onecrust') then
- REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
- else
- stop 'reference 1D Earth model unknown'
- endif
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- HONOR_1D_SPHERICAL_MOHO = .true.
- ONE_CRUST = .true.
-
- else if(MODEL == 'transversely_isotropic_prem_plus_3D_crust_2.0') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
-
- else if(MODEL == 's20rts') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
- THREE_D_MODEL = THREE_D_MODEL_S20RTS
-
- else if(MODEL == 's362ani') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
- THREE_D_MODEL = THREE_D_MODEL_S362ANI
-
- else if(MODEL == 's362iso') then
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
- THREE_D_MODEL = THREE_D_MODEL_S362ANI
-
- else if(MODEL == 's362wmani') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
- THREE_D_MODEL = THREE_D_MODEL_S362WMANI
-
- else if(MODEL == 's362ani_prem') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
- THREE_D_MODEL = THREE_D_MODEL_S362ANI_PREM
-
- else if(MODEL == 's29ea') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .true.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
- REFERENCE_1D_MODEL = REFERENCE_MODEL_REF
- THREE_D_MODEL = THREE_D_MODEL_S29EA
-
- else if(MODEL == '3D_attenuation') then
- TRANSVERSE_ISOTROPY = .false.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .true.
- ONE_CRUST = .true.
- CASE_3D = .true.
-
- else if(MODEL == '3D_anisotropic') then
- TRANSVERSE_ISOTROPY = .true.
- ISOTROPIC_3D_MANTLE = .false.
- ANISOTROPIC_3D_MANTLE = .true.
- ANISOTROPIC_INNER_CORE = .false.
- CRUSTAL = .false.
- ATTENUATION_3D = .false.
- ONE_CRUST = .true.
- CASE_3D = .true.
-
- else
- stop 'model not implemented, edit read_compute_parameters.f90 and recompile'
- endif
-
-! set time step, radial distribution of elements, and attenuation period range
-! right distribution is determined based upon maximum value of NEX
- NEX_MAX = max(NEX_XI,NEX_ETA)
-
-!----
-!---- case prem_onecrust by default
-!----
- if (SUPPRESS_CRUSTAL_MESH) then
- multiplication_factor=2
- else
- multiplication_factor=1
- endif
-
- ! element width = 0.5625000 degrees = 62.54715 km
- if(NEX_MAX*multiplication_factor <= 160) then
- DT = 0.252d0
-
- MIN_ATTENUATION_PERIOD = 30
- MAX_ATTENUATION_PERIOD = 1500
-
- NER_CRUST = 1
- NER_80_MOHO = 1
- NER_220_80 = 2
- NER_400_220 = 2
- NER_600_400 = 2
- NER_670_600 = 1
- NER_771_670 = 1
- NER_TOPDDOUBLEPRIME_771 = 15
- NER_CMB_TOPDDOUBLEPRIME = 1
- NER_OUTER_CORE = 16
- NER_TOP_CENTRAL_CUBE_ICB = 2
- R_CENTRAL_CUBE = 950000.d0
-
- ! element width = 0.3515625 degrees = 39.09196 km
- else if(NEX_MAX*multiplication_factor <= 256) then
- DT = 0.225d0
-
- MIN_ATTENUATION_PERIOD = 20
- MAX_ATTENUATION_PERIOD = 1000
-
- NER_CRUST = 1
- NER_80_MOHO = 1
- NER_220_80 = 2
- NER_400_220 = 3
- NER_600_400 = 3
- NER_670_600 = 1
- NER_771_670 = 1
- NER_TOPDDOUBLEPRIME_771 = 22
- NER_CMB_TOPDDOUBLEPRIME = 2
- NER_OUTER_CORE = 24
- NER_TOP_CENTRAL_CUBE_ICB = 3
- R_CENTRAL_CUBE = 965000.d0
-
- ! element width = 0.2812500 degrees = 31.27357 km
- else if(NEX_MAX*multiplication_factor <= 320) then
- DT = 0.16d0
-
- MIN_ATTENUATION_PERIOD = 15
- MAX_ATTENUATION_PERIOD = 750
-
- NER_CRUST = 1
- NER_80_MOHO = 1
- NER_220_80 = 3
- NER_400_220 = 4
- NER_600_400 = 4
- NER_670_600 = 1
- NER_771_670 = 2
- NER_TOPDDOUBLEPRIME_771 = 29
- NER_CMB_TOPDDOUBLEPRIME = 2
- NER_OUTER_CORE = 32
- NER_TOP_CENTRAL_CUBE_ICB = 4
- R_CENTRAL_CUBE = 940000.d0
-
- ! element width = 0.1875000 degrees = 20.84905 km
- else if(NEX_MAX*multiplication_factor <= 480) then
- DT = 0.11d0
-
- MIN_ATTENUATION_PERIOD = 10
- MAX_ATTENUATION_PERIOD = 500
-
- NER_CRUST = 1
- NER_80_MOHO = 2
- NER_220_80 = 4
- NER_400_220 = 5
- NER_600_400 = 6
- NER_670_600 = 2
- NER_771_670 = 2
- NER_TOPDDOUBLEPRIME_771 = 44
- NER_CMB_TOPDDOUBLEPRIME = 3
- NER_OUTER_CORE = 48
- NER_TOP_CENTRAL_CUBE_ICB = 5
- R_CENTRAL_CUBE = 988000.d0
-
- ! element width = 0.1757812 degrees = 19.54598 km
- else if(NEX_MAX*multiplication_factor <= 512) then
- DT = 0.1125d0
-
- MIN_ATTENUATION_PERIOD = 9
- MAX_ATTENUATION_PERIOD = 500
-
- NER_CRUST = 1
- NER_80_MOHO = 2
- NER_220_80 = 4
- NER_400_220 = 6
- NER_600_400 = 6
- NER_670_600 = 2
- NER_771_670 = 3
- NER_TOPDDOUBLEPRIME_771 = 47
- NER_CMB_TOPDDOUBLEPRIME = 3
- NER_OUTER_CORE = 51
- NER_TOP_CENTRAL_CUBE_ICB = 5
- R_CENTRAL_CUBE = 1010000.d0
-
- ! element width = 0.1406250 degrees = 15.63679 km
- else if(NEX_MAX*multiplication_factor <= 640) then
- DT = 0.09d0
-
- MIN_ATTENUATION_PERIOD = 8
- MAX_ATTENUATION_PERIOD = 400
-
- NER_CRUST = 2
- NER_80_MOHO = 3
- NER_220_80 = 5
- NER_400_220 = 7
- NER_600_400 = 8
- NER_670_600 = 3
- NER_771_670 = 3
- NER_TOPDDOUBLEPRIME_771 = 59
- NER_CMB_TOPDDOUBLEPRIME = 4
- NER_OUTER_CORE = 64
- NER_TOP_CENTRAL_CUBE_ICB = 6
- R_CENTRAL_CUBE = 1020000.d0
-
- ! element width = 0.1041667 degrees = 11.58280 km
- else if(NEX_MAX*multiplication_factor <= 864) then
- DT = 0.0667d0
-
- MIN_ATTENUATION_PERIOD = 6
- MAX_ATTENUATION_PERIOD = 300
-
- NER_CRUST = 2
- NER_80_MOHO = 4
- NER_220_80 = 6
- NER_400_220 = 10
- NER_600_400 = 10
- NER_670_600 = 3
- NER_771_670 = 4
- NER_TOPDDOUBLEPRIME_771 = 79
- NER_CMB_TOPDDOUBLEPRIME = 5
- NER_OUTER_CORE = 86
- NER_TOP_CENTRAL_CUBE_ICB = 9
- R_CENTRAL_CUBE = 990000.d0
-
- ! element width = 7.8125000E-02 degrees = 8.687103 km
- else if(NEX_MAX*multiplication_factor <= 1152) then
- DT = 0.05d0
-
- MIN_ATTENUATION_PERIOD = 4
- MAX_ATTENUATION_PERIOD = 200
-
- NER_CRUST = 3
- NER_80_MOHO = 6
- NER_220_80 = 8
- NER_400_220 = 13
- NER_600_400 = 13
- NER_670_600 = 4
- NER_771_670 = 6
- NER_TOPDDOUBLEPRIME_771 = 106
- NER_CMB_TOPDDOUBLEPRIME = 7
- NER_OUTER_CORE = 116
- NER_TOP_CENTRAL_CUBE_ICB = 12
- R_CENTRAL_CUBE = 985000.d0
-
- ! element width = 7.2115384E-02 degrees = 8.018865 km
- else if(NEX_MAX*multiplication_factor <= 1248) then
- DT = 0.0462d0
-
- MIN_ATTENUATION_PERIOD = 4
- MAX_ATTENUATION_PERIOD = 200
-
- NER_CRUST = 3
- NER_80_MOHO = 6
- NER_220_80 = 9
- NER_400_220 = 14
- NER_600_400 = 14
- NER_670_600 = 5
- NER_771_670 = 6
- NER_TOPDDOUBLEPRIME_771 = 114
- NER_CMB_TOPDDOUBLEPRIME = 8
- NER_OUTER_CORE = 124
- NER_TOP_CENTRAL_CUBE_ICB = 13
- R_CENTRAL_CUBE = 985000.d0
-
- else
-
-! scale with respect to 1248 if above that limit
- DT = 0.0462d0 * 1248.d0 / (2.d0*NEX_MAX)
-
- MIN_ATTENUATION_PERIOD = 4
- MAX_ATTENUATION_PERIOD = 200
-
- NER_CRUST = nint(3 * 2.d0*NEX_MAX / 1248.d0)
- NER_80_MOHO = nint(6 * 2.d0*NEX_MAX / 1248.d0)
- NER_220_80 = nint(9 * 2.d0*NEX_MAX / 1248.d0)
- NER_400_220 = nint(14 * 2.d0*NEX_MAX / 1248.d0)
- NER_600_400 = nint(14 * 2.d0*NEX_MAX / 1248.d0)
- NER_670_600 = nint(5 * 2.d0*NEX_MAX / 1248.d0)
- NER_771_670 = nint(6 * 2.d0*NEX_MAX / 1248.d0)
- NER_TOPDDOUBLEPRIME_771 = nint(114 * 2.d0*NEX_MAX / 1248.d0)
- NER_CMB_TOPDDOUBLEPRIME = nint(8 * 2.d0*NEX_MAX / 1248.d0)
- NER_OUTER_CORE = nint(124 * 2.d0*NEX_MAX / 1248.d0)
- NER_TOP_CENTRAL_CUBE_ICB = nint(13 * 2.d0*NEX_MAX / 1248.d0)
- R_CENTRAL_CUBE = 985000.d0
-
-!! removed this limit else
-!! removed this limit stop 'problem with this value of NEX_MAX'
- endif
-
-!----
-!---- change some values in the case of regular PREM with two crustal layers or of 3D models
-!----
-
-! case of regular PREM with two crustal layers: change the time step for small meshes
-! because of a different size of elements in the radial direction in the crust
- if (HONOR_1D_SPHERICAL_MOHO) then
- if (.not. ONE_CRUST) then
- ! case 1D + two crustal layers
- if (NER_CRUST<2) NER_CRUST=2
- if(NEX_MAX*multiplication_factor <= 160) then
- DT = 0.20d0
- else if(NEX_MAX*multiplication_factor <= 256) then
- DT = 0.20d0
- endif
- endif
- else
- ! case 3D
- if (NER_CRUST<2) NER_CRUST=2
- if(NEX_MAX*multiplication_factor <= 160) then
- DT = 0.15d0
- else if(NEX_MAX*multiplication_factor <= 256) then
- DT = 0.17d0
- else if(NEX_MAX*multiplication_factor <= 320) then
- DT = 0.155d0
- endif
- endif
-
- if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
- DT = DT*0.20d0
- endif
-
-
- if( .not. ATTENUATION_RANGE_PREDEFINED ) then
- call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
- MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD)
- endif
-
- if(ANGULAR_WIDTH_XI_IN_DEGREES < 90.0d0 .or. &
- ANGULAR_WIDTH_ETA_IN_DEGREES < 90.0d0 .or. &
- NEX_MAX > 1248) then
-
- call auto_ner(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
- NER_CRUST, NER_80_MOHO, NER_220_80, NER_400_220, NER_600_400, &
- NER_670_600, NER_771_670, NER_TOPDDOUBLEPRIME_771, &
- NER_CMB_TOPDDOUBLEPRIME, NER_OUTER_CORE, NER_TOP_CENTRAL_CUBE_ICB, &
- R_CENTRAL_CUBE, CASE_3D)
-
- call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
- MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD)
-
- call auto_time_stepping(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, DT)
-
-goto 777
- write(*,*)'##############################################################'
- write(*,*)
- write(*,*)' Auto Radial Meshing Code '
- write(*,*)' Consult read_compute_parameters.f90 and auto_ner.f90 '
- write(*,*)' This should only be invoked for chunks less than 90 degrees'
- write(*,*)' and for chunks greater than 1248 elements wide'
- write(*,*)
- write(*,*)'CHUNK WIDTH: ', ANGULAR_WIDTH_XI_IN_DEGREES
- write(*,*)'NEX: ', NEX_MAX
- write(*,*)'NER_CRUST: ', NER_CRUST
- write(*,*)'NER_80_MOHO: ', NER_80_MOHO
- write(*,*)'NER_220_80: ', NER_220_80
- write(*,*)'NER_400_220: ', NER_400_220
- write(*,*)'NER_600_400: ', NER_600_400
- write(*,*)'NER_670_600: ', NER_670_600
- write(*,*)'NER_771_670: ', NER_771_670
- write(*,*)'NER_TOPDDOUBLEPRIME_771: ', NER_TOPDDOUBLEPRIME_771
- write(*,*)'NER_CMB_TOPDDOUBLEPRIME: ', NER_CMB_TOPDDOUBLEPRIME
- write(*,*)'NER_OUTER_CORE: ', NER_OUTER_CORE
- write(*,*)'NER_TOP_CENTRAL_CUBE_ICB: ', NER_TOP_CENTRAL_CUBE_ICB
- write(*,*)'R_CENTRAL_CUBE: ', R_CENTRAL_CUBE
- write(*,*)'multiplication factor: ', multiplication_factor
- write(*,*)
- write(*,*)'DT: ',DT
- write(*,*)'MIN_ATTENUATION_PERIOD ',MIN_ATTENUATION_PERIOD
- write(*,*)'MAX_ATTENUATION_PERIOD ',MAX_ATTENUATION_PERIOD
- write(*,*)
- write(*,*)'##############################################################'
-777 continue
-
- if (HONOR_1D_SPHERICAL_MOHO) then
- if (.not. ONE_CRUST) then
- ! case 1D + two crustal layers
- if (NER_CRUST<2) NER_CRUST=2
- endif
- else
- ! case 3D
- if (NER_CRUST<2) NER_CRUST=2
- endif
- endif
-
-
-! take a 5% safety margin on the maximum stable time step
-! which was obtained by trial and error
- DT = DT * (1.d0 - 0.05d0)
-
- call read_value_logical(OCEANS, 'model.OCEANS')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(ELLIPTICITY, 'model.ELLIPTICITY')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(TOPOGRAPHY, 'model.TOPOGRAPHY')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(GRAVITY, 'model.GRAVITY')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(ROTATION, 'model.ROTATION')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(ATTENUATION, 'model.ATTENUATION')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
- call read_value_logical(ABSORBING_CONDITIONS, 'solver.ABSORBING_CONDITIONS')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
- if(ABSORBING_CONDITIONS .and. NCHUNKS == 6) stop 'cannot have absorbing conditions in the full Earth'
-
- if(ABSORBING_CONDITIONS .and. NCHUNKS == 3) stop 'absorbing conditions not supported for three chunks yet'
-
- if(ATTENUATION_3D .and. .not. ATTENUATION) stop 'need ATTENUATION to use ATTENUATION_3D'
-
-! radii in PREM or IASP91
-! and normalized density at fluid-solid interface on fluid size for coupling
-! ROCEAN: radius of the ocean (m)
-! RMIDDLE_CRUST: radius of the middle crust (m)
-! RMOHO: radius of the Moho (m)
-! R80: radius of 80 km discontinuity (m)
-! R120: radius of 120 km discontinuity (m) in IASP91
-! R220: radius of 220 km discontinuity (m)
-! R400: radius of 400 km discontinuity (m)
-! R600: radius of 600 km 2nd order discontinuity (m)
-! R670: radius of 670 km discontinuity (m)
-! R771: radius of 771 km 2nd order discontinuity (m)
-! RTOPDDOUBLEPRIME: radius of top of D" 2nd order discontinuity (m)
-! RCMB: radius of CMB (m)
-! RICB: radius of ICB (m)
-
-! by default there is no d120 discontinuity, except in IASP91, therefore set to fictitious value
- R120 = -1.d0
-
-! value common to all models
- RHO_OCEANS = 1020.0 / RHOAV
-
- if(REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91) then
-
-! IASP91
- ROCEAN = 6371000.d0
- RMIDDLE_CRUST = 6351000.d0
- RMOHO = 6336000.d0
- R80 = 6291000.d0
- R120 = 6251000.d0
- R220 = 6161000.d0
- R400 = 5961000.d0
-! there is no d600 discontinuity in IASP91 therefore this value is useless
-! but it needs to be there for compatibility with other subroutines
- R600 = R_EARTH - 600000.d0
- R670 = 5711000.d0
- R771 = 5611000.d0
- RTOPDDOUBLEPRIME = 3631000.d0
- RCMB = 3482000.d0
- RICB = 1217000.d0
-
- RHO_TOP_OC = 9900.2379 / RHOAV
- RHO_BOTTOM_OC = 12168.6383 / RHOAV
-
- else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
-
-! our implementation of AK135 has not been checked carefully yet
-! therefore let us doublecheck it carefully one day
-
-! values below corrected by Ying Zhou <yingz at gps.caltech.edu>
-
-! AK135 without the 300 meters of mud layer
- ROCEAN = 6368000.d0
- RMIDDLE_CRUST = 6361000.d0
- RMOHO = 6353000.d0
- R80 = 6291000.d0
- R220 = 6161000.d0
- R400 = 5961000.d0
- R670 = 5711000.d0
- RTOPDDOUBLEPRIME = 3631000.d0
- RCMB = 3479500.d0
- RICB = 1217500.d0
-
-! values for AK135 that are not discontinuities
- R600 = 5771000.d0
- R771 = 5611000.d0
-
- RHO_TOP_OC = 9914.5000 / RHOAV
- RHO_BOTTOM_OC = 12139.1000 / RHOAV
-
- else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
-
-! values below corrected by Ying Zhou <yingz at gps.caltech.edu>
-
-! 1066A
- RMOHO = 6360000.d0
- R400 = 5950000.d0
- R600 = 5781000.d0
- R670 = 5700000.d0
- RCMB = 3484300.d0
- RICB = 1229480.d0
-
-! values for 1066A that are not discontinuities
- RTOPDDOUBLEPRIME = 3631000.d0
- R220 = 6161000.d0
- R771 = 5611000.d0
-! RMIDDLE_CRUST used only for high resolution FFSW1C model, with 3 elements crust simulations
-! mid_crust = 10 km
- RMIDDLE_CRUST = 6361000.d0
- R80 = 6291000.d0
-
-! model 1066A has no oceans, therefore we use the radius of the Earth instead
- ROCEAN = R_EARTH
-
- RHO_TOP_OC = 9917.4500 / RHOAV
- RHO_BOTTOM_OC = 12160.6500 / RHOAV
-
- else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
-
-! REF
- ROCEAN = 6368000.d0
- RMIDDLE_CRUST = 6356000.d0
- RMOHO = 6346600.d0
- R80 = 6291000.d0
- R220 = 6151000.d0
- R400 = 5961000.d0
- R600 = 5771000.d0
- R670 = 5721000.d0
- R771 = 5600000.d0
- RTOPDDOUBLEPRIME = 3630000.d0
- RCMB = 3479958.d0
- RICB = 1221491.d0
-
- RHO_TOP_OC = 9903.48 / RHOAV
- RHO_BOTTOM_OC = 12166.35 / RHOAV
-
- else
-
-! PREM
- ROCEAN = 6368000.d0
- RMIDDLE_CRUST = 6356000.d0
- RMOHO = 6346600.d0
- R80 = 6291000.d0
- R220 = 6151000.d0
- R400 = 5971000.d0
- R600 = 5771000.d0
- R670 = 5701000.d0
- R771 = 5600000.d0
- RTOPDDOUBLEPRIME = 3630000.d0
- RCMB = 3480000.d0
- RICB = 1221000.d0
-
- RHO_TOP_OC = 9903.4384 / RHOAV
- RHO_BOTTOM_OC = 12166.5885 / RHOAV
-
- endif
-
-! honor the PREM Moho or define a fictitious Moho in order to have even radial sampling
-! from the d220 to the Earth surface
- if(HONOR_1D_SPHERICAL_MOHO) then
- RMOHO_FICTITIOUS_IN_MESHER = RMOHO
- else
- RMOHO_FICTITIOUS_IN_MESHER = (R80 + R_EARTH) / 2
- endif
-
- call read_value_double_precision(RECORD_LENGTH_IN_MINUTES, 'solver.RECORD_LENGTH_IN_MINUTES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
-! compute total number of time steps, rounded to next multiple of 100
- NSTEP = 100 * (int(RECORD_LENGTH_IN_MINUTES * 60.d0 / (100.d0*DT)) + 1)
-
- call read_value_logical(MOVIE_SURFACE, 'solver.MOVIE_SURFACE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(MOVIE_VOLUME, 'solver.MOVIE_VOLUME')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NTSTEP_BETWEEN_FRAMES, 'solver.NTSTEP_BETWEEN_FRAMES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(HDUR_MOVIE, 'solver.HDUR_MOVIE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
-! computes a default hdur_movie that creates nice looking movies.
-! Sets HDUR_MOVIE as the minimum period the mesh can resolve
- if(HDUR_MOVIE <= TINYVAL) &
- HDUR_MOVIE = 1.1d0*max(240.d0/NEX_XI*18.d0*ANGULAR_WIDTH_XI_IN_DEGREES/90.d0, &
- 240.d0/NEX_ETA*18.d0*ANGULAR_WIDTH_ETA_IN_DEGREES/90.d0)
-
- call read_value_integer(MOVIE_VOLUME_TYPE, 'solver.MOVIE_VOLUME_TYPE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(MOVIE_VOLUME_COARSE,'solver.MOVIE_VOLUME_COARSE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_TOP_KM, 'solver.MOVIE_TOP_KM')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_BOTTOM_KM, 'solver.MOVIE_BOTTOM_KM')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_WEST_DEG, 'solver.MOVIE_WEST_DEG')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_EAST_DEG, 'solver.MOVIE_EAST_DEG')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_NORTH_DEG, 'solver.MOVIE_NORTH_DEG')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_double_precision(MOVIE_SOUTH_DEG, 'solver.MOVIE_SOUTH_DEG')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(MOVIE_START, 'solver.MOVIE_START')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(MOVIE_STOP, 'solver.MOVIE_STOP')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- MOVIE_TOP = (R_EARTH_KM-MOVIE_TOP_KM)/R_EARTH_KM
- MOVIE_BOTTOM = (R_EARTH_KM-MOVIE_BOTTOM_KM)/R_EARTH_KM
- MOVIE_EAST = MOVIE_EAST_DEG * DEGREES_TO_RADIANS
- MOVIE_WEST = MOVIE_WEST_DEG * DEGREES_TO_RADIANS
- MOVIE_NORTH = (90.0d0 - MOVIE_NORTH_DEG) * DEGREES_TO_RADIANS ! converting from latitude to colatitude
- MOVIE_SOUTH = (90.0d0 - MOVIE_SOUTH_DEG) * DEGREES_TO_RADIANS
-
- call read_value_logical(SAVE_MESH_FILES, 'mesher.SAVE_MESH_FILES')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NUMBER_OF_RUNS, 'solver.NUMBER_OF_RUNS')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NUMBER_OF_THIS_RUN, 'solver.NUMBER_OF_THIS_RUN')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_string(LOCAL_PATH, 'LOCAL_PATH')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NTSTEP_BETWEEN_OUTPUT_INFO, 'solver.NTSTEP_BETWEEN_OUTPUT_INFO')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NTSTEP_BETWEEN_OUTPUT_SEISMOS, 'solver.NTSTEP_BETWEEN_OUTPUT_SEISMOS')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NTSTEP_BETWEEN_READ_ADJSRC, 'solver.NTSTEP_BETWEEN_READ_ADJSRC')
- if(err_occurred() /= 0) return
-
- call read_value_logical(OUTPUT_SEISMOS_ASCII_TEXT, 'solver.OUTPUT_SEISMOS_ASCII_TEXT')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(OUTPUT_SEISMOS_SAC_ALPHANUM, 'solver.OUTPUT_SEISMOS_SAC_ALPHANUM')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(OUTPUT_SEISMOS_SAC_BINARY, 'solver.OUTPUT_SEISMOS_SAC_BINARY')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(ROTATE_SEISMOGRAMS_RT, 'solver.ROTATE_SEISMOGRAMS_RT')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(WRITE_SEISMOGRAMS_BY_MASTER, 'solver.WRITE_SEISMOGRAMS_BY_MASTER')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(SAVE_ALL_SEISMOS_IN_ONE_FILE, 'solver.SAVE_ALL_SEISMOS_IN_ONE_FILE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(USE_BINARY_FOR_LARGE_FILE, 'solver.USE_BINARY_FOR_LARGE_FILE')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
- call read_value_logical(RECEIVERS_CAN_BE_BURIED, 'solver.RECEIVERS_CAN_BE_BURIED')
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_logical(PRINT_SOURCE_TIME_FUNCTION, 'solver.PRINT_SOURCE_TIME_FUNCTION')
-
- if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
-
-! close parameter file
- call close_parameter_file
-!--- check that parameters make sense
-
- if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
- stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5, check the Par_file'
-
-! subsets used to save seismograms must not be larger than the whole time series,
-! otherwise we waste memory
- if(NTSTEP_BETWEEN_OUTPUT_SEISMOS > NSTEP) then
- NTSTEP_BETWEEN_OUTPUT_SEISMOS = NSTEP
- if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
- stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then modified NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5'
- endif
-
-! check that reals are either 4 or 8 bytes
- if(CUSTOM_REAL /= SIZE_REAL .and. CUSTOM_REAL /= SIZE_DOUBLE) stop 'wrong size of CUSTOM_REAL for reals'
-
-! check that the parameter file is correct
- if(NGNOD /= 27) stop 'number of control nodes must be 27'
- if(NGNOD == 27 .and. NGNOD2D /= 9) stop 'elements with 27 points should have NGNOD2D = 9'
-
-! for the number of standard linear solids for attenuation
- if(N_SLS /= 3) stop 'number of SLS must be 3'
-
-! check number of slices in each direction
- if(NCHUNKS < 1) stop 'must have at least one chunk'
- if(NPROC_XI < 1) stop 'NPROC_XI must be at least 1'
- if(NPROC_ETA < 1) stop 'NPROC_ETA must be at least 1'
-
-! check number of chunks
- if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
- stop 'only one, two, three or six chunks can be meshed'
-
-! check that the central cube can be included
- if(INCLUDE_CENTRAL_CUBE .and. NCHUNKS /= 6) stop 'need six chunks to include central cube'
-
-! check that sphere can be cut into slices without getting negative Jacobian
- if(NEX_XI < 48) stop 'NEX_XI must be greater than 48 to cut the sphere into slices with positive Jacobian'
- if(NEX_ETA < 48) stop 'NEX_ETA must be greater than 48 to cut the sphere into slices with positive Jacobian'
-
-! check that mesh can be coarsened in depth three or four times
- CUT_SUPERBRICK_XI=.false.
- CUT_SUPERBRICK_ETA=.false.
-
- if (SUPPRESS_CRUSTAL_MESH .and. .not. ADD_4TH_DOUBLING) then
- if(mod(NEX_XI,8) /= 0) stop 'NEX_XI must be a multiple of 8'
- if(mod(NEX_ETA,8) /= 0) stop 'NEX_ETA must be a multiple of 8'
- if(mod(NEX_XI/4,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 4*NPROC_XI'
- if(mod(NEX_ETA/4,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 4*NPROC_ETA'
- if(mod(NEX_XI/8,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
- if(mod(NEX_ETA/8,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
- elseif (SUPPRESS_CRUSTAL_MESH .or. .not. ADD_4TH_DOUBLING) then
- if(mod(NEX_XI,16) /= 0) stop 'NEX_XI must be a multiple of 16'
- if(mod(NEX_ETA,16) /= 0) stop 'NEX_ETA must be a multiple of 16'
- if(mod(NEX_XI/8,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 8*NPROC_XI'
- if(mod(NEX_ETA/8,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 8*NPROC_ETA'
- if(mod(NEX_XI/16,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
- if(mod(NEX_ETA/16,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
- else
- if(mod(NEX_XI,32) /= 0) stop 'NEX_XI must be a multiple of 32'
- if(mod(NEX_ETA,32) /= 0) stop 'NEX_ETA must be a multiple of 32'
- if(mod(NEX_XI/16,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 16*NPROC_XI'
- if(mod(NEX_ETA/16,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 16*NPROC_ETA'
- if(mod(NEX_XI/32,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
- if(mod(NEX_ETA/32,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
- endif
-
-! check that topology is correct if more than two chunks
- if(NCHUNKS > 2 .and. NEX_XI /= NEX_ETA) stop 'must have NEX_XI = NEX_ETA for more than two chunks'
- if(NCHUNKS > 2 .and. NPROC_XI /= NPROC_ETA) stop 'must have NPROC_XI = NPROC_ETA for more than two chunks'
-
-! check that IASP91, AK135, or 1066A is isotropic
- if((REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91 .or. &
- REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135 .or. &
- REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) .and. TRANSVERSE_ISOTROPY) &
- stop 'models IASP91, AK135 and 1066A are currently isotropic'
-
- ELEMENT_WIDTH = ANGULAR_WIDTH_XI_IN_DEGREES/dble(NEX_MAX) * DEGREES_TO_RADIANS
-
-!
-!--- compute additional parameters
-!
-
-! number of elements horizontally in each slice (i.e. per processor)
-! these two values MUST be equal in all cases
- NEX_PER_PROC_XI = NEX_XI / NPROC_XI
- NEX_PER_PROC_ETA = NEX_ETA / NPROC_ETA
-
-! total number of processors in each of the six chunks
- NPROC = NPROC_XI * NPROC_ETA
-
-! total number of processors in the full Earth composed of the six chunks
- NPROCTOT = NCHUNKS * NPROC
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! definition of general mesh parameters below
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! find element below top of which we should implement the second doubling in the mantle
-! locate element closest to optimal value
- distance_min = HUGEVAL
- do ielem = 2,NER_TOPDDOUBLEPRIME_771
- zval = RTOPDDOUBLEPRIME + ielem * (R771 - RTOPDDOUBLEPRIME) / dble(NER_TOPDDOUBLEPRIME_771)
- distance = abs(zval - (R_EARTH - DEPTH_SECOND_DOUBLING_OPTIMAL))
- if(distance < distance_min) then
- elem_doubling_mantle = ielem
- distance_min = distance
- DEPTH_SECOND_DOUBLING_REAL = R_EARTH - zval
- endif
- enddo
-
-! find element below top of which we should implement the third doubling in the middle of the outer core
-! locate element closest to optimal value
- distance_min = HUGEVAL
-! start at element number 4 because we need at least two elements below for the fourth doubling
-! implemented at the bottom of the outer core
- do ielem = 4,NER_OUTER_CORE
- zval = RICB + ielem * (RCMB - RICB) / dble(NER_OUTER_CORE)
- distance = abs(zval - (R_EARTH - DEPTH_THIRD_DOUBLING_OPTIMAL))
- if(distance < distance_min) then
- elem_doubling_middle_outer_core = ielem
- distance_min = distance
- DEPTH_THIRD_DOUBLING_REAL = R_EARTH - zval
- endif
- enddo
-
- if (ADD_4TH_DOUBLING) then
-! find element below top of which we should implement the fourth doubling in the middle of the outer core
-! locate element closest to optimal value
- distance_min = HUGEVAL
-! end two elements before the top because we need at least two elements above for the third doubling
-! implemented in the middle of the outer core
- do ielem = 2,NER_OUTER_CORE-2
- zval = RICB + ielem * (RCMB - RICB) / dble(NER_OUTER_CORE)
- distance = abs(zval - (R_EARTH - DEPTH_FOURTH_DOUBLING_OPTIMAL))
- if(distance < distance_min) then
- elem_doubling_bottom_outer_core = ielem
- distance_min = distance
- DEPTH_FOURTH_DOUBLING_REAL = R_EARTH - zval
- endif
- enddo
-! make sure that the two doublings in the outer core are found in the right order
- if(elem_doubling_bottom_outer_core >= elem_doubling_middle_outer_core) &
- stop 'error in location of the two doublings in the outer core'
- endif
-
- ratio_sampling_array(15) = 0
-
-! define all the layers of the mesh
- if (.not. ADD_4TH_DOUBLING) then
-
- if (SUPPRESS_CRUSTAL_MESH) then
-
- ONE_CRUST = .false.
- OCEANS= .false.
- TOPOGRAPHY = .false.
- CRUSTAL = .false.
-
- NUMBER_OF_MESH_LAYERS = 14
- layer_offset = 1
-
- ! now only one region
- ner( 1) = NER_CRUST + NER_80_MOHO
- ner( 2) = 0
- ner( 3) = 0
-
- ner( 4) = NER_220_80
- ner( 5) = NER_400_220
- ner( 6) = NER_600_400
- ner( 7) = NER_670_600
- ner( 8) = NER_771_670
- ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner(10) = elem_doubling_mantle
- ner(11) = NER_CMB_TOPDDOUBLEPRIME
- ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(13) = elem_doubling_middle_outer_core
- ner(14) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1:9) = 1
- ratio_sampling_array(10:12) = 2
- ratio_sampling_array(13:14) = 4
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1:3) = IFLAG_CRUST !!!!! IFLAG_80_MOHO
- doubling_index(4) = IFLAG_220_80
- doubling_index(5:7) = IFLAG_670_220
- doubling_index(8:11) = IFLAG_MANTLE_NORMAL
- doubling_index(12:13) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(14) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(10) = .true.
- this_region_has_a_doubling(13) = .true.
- lastdoubling_layer = 13
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- r_top(1) = R_EARTH
- r_bottom(1) = R80
-
- r_top(2) = RMIDDLE_CRUST !!!! now fictitious
- r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
-
- r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_bottom(3) = R80 !!!! now fictitious
-
- r_top(4) = R80
- r_bottom(4) = R220
-
- r_top(5) = R220
- r_bottom(5) = R400
-
- r_top(6) = R400
- r_bottom(6) = R600
-
- r_top(7) = R600
- r_bottom(7) = R670
-
- r_top(8) = R670
- r_bottom(8) = R771
-
- r_top(9) = R771
- r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(10) = RTOPDDOUBLEPRIME
-
- r_top(11) = RTOPDDOUBLEPRIME
- r_bottom(11) = RCMB
-
- r_top(12) = RCMB
- r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(13) = RICB
-
- r_top(15) = RICB
- r_bottom(15) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = R80 / R_EARTH
-
- rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
- rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
-
- rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmins(3) = R80 / R_EARTH !!!! now fictitious
-
- rmaxs(4) = R80 / R_EARTH
- rmins(4) = R220 / R_EARTH
-
- rmaxs(5) = R220 / R_EARTH
- rmins(5) = R400 / R_EARTH
-
- rmaxs(6) = R400 / R_EARTH
- rmins(6) = R600 / R_EARTH
-
- rmaxs(7) = R600 / R_EARTH
- rmins(7) = R670 / R_EARTH
-
- rmaxs(8) = R670 / R_EARTH
- rmins(8) = R771 / R_EARTH
-
- rmaxs(9:10) = R771 / R_EARTH
- rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(11) = RCMB / R_EARTH
-
- rmaxs(12:13) = RCMB / R_EARTH
- rmins(12:13) = RICB / R_EARTH
-
- rmaxs(14) = RICB / R_EARTH
- rmins(14) = R_CENTRAL_CUBE / R_EARTH
-
- elseif (ONE_CRUST) then
-
- NUMBER_OF_MESH_LAYERS = 13
- layer_offset = 0
-
- ner( 1) = NER_CRUST
- ner( 2) = NER_80_MOHO
- ner( 3) = NER_220_80
- ner( 4) = NER_400_220
- ner( 5) = NER_600_400
- ner( 6) = NER_670_600
- ner( 7) = NER_771_670
- ner( 8) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner( 9) = elem_doubling_mantle
- ner(10) = NER_CMB_TOPDDOUBLEPRIME
- ner(11) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(12) = elem_doubling_middle_outer_core
- ner(13) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1) = 1
- ratio_sampling_array(2:8) = 2
- ratio_sampling_array(9:11) = 4
- ratio_sampling_array(12:13) = 8
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1) = IFLAG_CRUST
- doubling_index(2) = IFLAG_80_MOHO
- doubling_index(3) = IFLAG_220_80
- doubling_index(4:6) = IFLAG_670_220
- doubling_index(7:10) = IFLAG_MANTLE_NORMAL
- doubling_index(11:12) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(13) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(2) = .true.
- this_region_has_a_doubling(9) = .true.
- this_region_has_a_doubling(12) = .true.
- lastdoubling_layer = 12
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- !!!!!!!!!!! DK DK UGLY: beware, is there a bug when 3D crust crosses anisotropy in the mantle?
- !!!!!!!!!!! DK DK UGLY: i.e. if there is no thick crust there, some elements above the Moho
- !!!!!!!!!!! DK DK UGLY: should be anisotropic but anisotropy is currently only
- !!!!!!!!!!! DK DK UGLY: stored between d220 and MOHO to save memory? Clarify this one day.
- !!!!!!!!!!! DK DK UGLY: The Moho stretching and squishing that Jeroen added to V4.0
- !!!!!!!!!!! DK DK UGLY: should partly deal with this problem.
-
- r_top(1) = R_EARTH
- r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
-
- r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(2) = R80
-
- r_top(3) = R80
- r_bottom(3) = R220
-
- r_top(4) = R220
- r_bottom(4) = R400
-
- r_top(5) = R400
- r_bottom(5) = R600
-
- r_top(6) = R600
- r_bottom(6) = R670
-
- r_top(7) = R670
- r_bottom(7) = R771
-
- r_top(8) = R771
- r_bottom(8) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(9) = RTOPDDOUBLEPRIME
-
- r_top(10) = RTOPDDOUBLEPRIME
- r_bottom(10) = RCMB
-
- r_top(11) = RCMB
- r_bottom(11) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(12) = RICB
-
- r_top(13) = RICB
- r_bottom(13) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
-
- rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(2) = R80 / R_EARTH
-
- rmaxs(3) = R80 / R_EARTH
- rmins(3) = R220 / R_EARTH
-
- rmaxs(4) = R220 / R_EARTH
- rmins(4) = R400 / R_EARTH
-
- rmaxs(5) = R400 / R_EARTH
- rmins(5) = R600 / R_EARTH
-
- rmaxs(6) = R600 / R_EARTH
- rmins(6) = R670 / R_EARTH
-
- rmaxs(7) = R670 / R_EARTH
- rmins(7) = R771 / R_EARTH
-
- rmaxs(8:9) = R771 / R_EARTH
- rmins(8:9) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(10) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(10) = RCMB / R_EARTH
-
- rmaxs(11:12) = RCMB / R_EARTH
- rmins(11:12) = RICB / R_EARTH
-
- rmaxs(13) = RICB / R_EARTH
- rmins(13) = R_CENTRAL_CUBE / R_EARTH
-
- else
-
- NUMBER_OF_MESH_LAYERS = 14
- layer_offset = 1
- if ((RMIDDLE_CRUST-RMOHO_FICTITIOUS_IN_MESHER)<(R_EARTH-RMIDDLE_CRUST)) then
- ner( 1) = ceiling (NER_CRUST / 2.d0)
- ner( 2) = floor (NER_CRUST / 2.d0)
- else
- ner( 1) = floor (NER_CRUST / 2.d0)
- ner( 2) = ceiling (NER_CRUST / 2.d0)
- endif
- ner( 3) = NER_80_MOHO
- ner( 4) = NER_220_80
- ner( 5) = NER_400_220
- ner( 6) = NER_600_400
- ner( 7) = NER_670_600
- ner( 8) = NER_771_670
- ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner(10) = elem_doubling_mantle
- ner(11) = NER_CMB_TOPDDOUBLEPRIME
- ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(13) = elem_doubling_middle_outer_core
- ner(14) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1:2) = 1
- ratio_sampling_array(3:9) = 2
- ratio_sampling_array(10:12) = 4
- ratio_sampling_array(13:14) = 8
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1:2) = IFLAG_CRUST
- doubling_index(3) = IFLAG_80_MOHO
- doubling_index(4) = IFLAG_220_80
- doubling_index(5:7) = IFLAG_670_220
- doubling_index(8:11) = IFLAG_MANTLE_NORMAL
- doubling_index(12:13) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(14) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(3) = .true.
- this_region_has_a_doubling(10) = .true.
- this_region_has_a_doubling(13) = .true.
- this_region_has_a_doubling(14) = .false.
- lastdoubling_layer = 13
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- r_top(1) = R_EARTH
- r_bottom(1) = RMIDDLE_CRUST
-
- r_top(2) = RMIDDLE_CRUST
- r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
-
- r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(3) = R80
-
- r_top(4) = R80
- r_bottom(4) = R220
-
- r_top(5) = R220
- r_bottom(5) = R400
-
- r_top(6) = R400
- r_bottom(6) = R600
-
- r_top(7) = R600
- r_bottom(7) = R670
-
- r_top(8) = R670
- r_bottom(8) = R771
-
- r_top(9) = R771
- r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(10) = RTOPDDOUBLEPRIME
-
- r_top(11) = RTOPDDOUBLEPRIME
- r_bottom(11) = RCMB
-
- r_top(12) = RCMB
- r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(13) = RICB
-
- r_top(14) = RICB
- r_bottom(14) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = RMIDDLE_CRUST / R_EARTH
-
- rmaxs(2) = RMIDDLE_CRUST / R_EARTH
- rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
-
- rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(3) = R80 / R_EARTH
-
- rmaxs(4) = R80 / R_EARTH
- rmins(4) = R220 / R_EARTH
-
- rmaxs(5) = R220 / R_EARTH
- rmins(5) = R400 / R_EARTH
-
- rmaxs(6) = R400 / R_EARTH
- rmins(6) = R600 / R_EARTH
-
- rmaxs(7) = R600 / R_EARTH
- rmins(7) = R670 / R_EARTH
-
- rmaxs(8) = R670 / R_EARTH
- rmins(8) = R771 / R_EARTH
-
- rmaxs(9:10) = R771 / R_EARTH
- rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(11) = RCMB / R_EARTH
-
- rmaxs(12:13) = RCMB / R_EARTH
- rmins(12:13) = RICB / R_EARTH
-
- rmaxs(14) = RICB / R_EARTH
- rmins(14) = R_CENTRAL_CUBE / R_EARTH
-
- endif
- else
- if (SUPPRESS_CRUSTAL_MESH) then
-
- ONE_CRUST = .false.
- OCEANS= .false.
- TOPOGRAPHY = .false.
- CRUSTAL = .false.
-
- NUMBER_OF_MESH_LAYERS = 15
- layer_offset = 1
-
- ! now only one region
- ner( 1) = NER_CRUST + NER_80_MOHO
- ner( 2) = 0
- ner( 3) = 0
-
- ner( 4) = NER_220_80
- ner( 5) = NER_400_220
- ner( 6) = NER_600_400
- ner( 7) = NER_670_600
- ner( 8) = NER_771_670
- ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner(10) = elem_doubling_mantle
- ner(11) = NER_CMB_TOPDDOUBLEPRIME
- ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(13) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
- ner(14) = elem_doubling_bottom_outer_core
- ner(15) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1:9) = 1
- ratio_sampling_array(10:12) = 2
- ratio_sampling_array(13) = 4
- ratio_sampling_array(14:15) = 8
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1:3) = IFLAG_CRUST !!!!! IFLAG_80_MOHO
- doubling_index(4) = IFLAG_220_80
- doubling_index(5:7) = IFLAG_670_220
- doubling_index(8:11) = IFLAG_MANTLE_NORMAL
- doubling_index(12:14) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(15) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(10) = .true.
- this_region_has_a_doubling(13) = .true.
- this_region_has_a_doubling(14) = .true.
- lastdoubling_layer = 14
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- r_top(1) = R_EARTH
- r_bottom(1) = R80
-
- r_top(2) = RMIDDLE_CRUST !!!! now fictitious
- r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
-
- r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_bottom(3) = R80 !!!! now fictitious
-
- r_top(4) = R80
- r_bottom(4) = R220
-
- r_top(5) = R220
- r_bottom(5) = R400
-
- r_top(6) = R400
- r_bottom(6) = R600
-
- r_top(7) = R600
- r_bottom(7) = R670
-
- r_top(8) = R670
- r_bottom(8) = R771
-
- r_top(9) = R771
- r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(10) = RTOPDDOUBLEPRIME
-
- r_top(11) = RTOPDDOUBLEPRIME
- r_bottom(11) = RCMB
-
- r_top(12) = RCMB
- r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
-
- r_top(14) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
- r_bottom(14) = RICB
-
- r_top(15) = RICB
- r_bottom(15) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = R80 / R_EARTH
-
- rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
- rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
-
- rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmins(3) = R80 / R_EARTH !!!! now fictitious
-
- rmaxs(4) = R80 / R_EARTH
- rmins(4) = R220 / R_EARTH
-
- rmaxs(5) = R220 / R_EARTH
- rmins(5) = R400 / R_EARTH
-
- rmaxs(6) = R400 / R_EARTH
- rmins(6) = R600 / R_EARTH
-
- rmaxs(7) = R600 / R_EARTH
- rmins(7) = R670 / R_EARTH
-
- rmaxs(8) = R670 / R_EARTH
- rmins(8) = R771 / R_EARTH
-
- rmaxs(9:10) = R771 / R_EARTH
- rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(11) = RCMB / R_EARTH
-
- rmaxs(12:14) = RCMB / R_EARTH
- rmins(12:14) = RICB / R_EARTH
-
- rmaxs(15) = RICB / R_EARTH
- rmins(15) = R_CENTRAL_CUBE / R_EARTH
-
- elseif (ONE_CRUST) then
-
- NUMBER_OF_MESH_LAYERS = 14
- layer_offset = 0
-
- ner( 1) = NER_CRUST
- ner( 2) = NER_80_MOHO
- ner( 3) = NER_220_80
- ner( 4) = NER_400_220
- ner( 5) = NER_600_400
- ner( 6) = NER_670_600
- ner( 7) = NER_771_670
- ner( 8) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner( 9) = elem_doubling_mantle
- ner(10) = NER_CMB_TOPDDOUBLEPRIME
- ner(11) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(12) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
- ner(13) = elem_doubling_bottom_outer_core
- ner(14) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1) = 1
- ratio_sampling_array(2:8) = 2
- ratio_sampling_array(9:11) = 4
- ratio_sampling_array(12) = 8
- ratio_sampling_array(13:14) = 16
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1) = IFLAG_CRUST
- doubling_index(2) = IFLAG_80_MOHO
- doubling_index(3) = IFLAG_220_80
- doubling_index(4:6) = IFLAG_670_220
- doubling_index(7:10) = IFLAG_MANTLE_NORMAL
- doubling_index(11:13) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(14) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(2) = .true.
- this_region_has_a_doubling(9) = .true.
- this_region_has_a_doubling(12) = .true.
- this_region_has_a_doubling(13) = .true.
- lastdoubling_layer = 13
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- !!!!!!!!!!! DK DK UGLY: beware, is there a bug when 3D crust crosses anisotropy in the mantle?
- !!!!!!!!!!! DK DK UGLY: i.e. if there is no thick crust there, some elements above the Moho
- !!!!!!!!!!! DK DK UGLY: should be anisotropic but anisotropy is currently only
- !!!!!!!!!!! DK DK UGLY: stored between d220 and MOHO to save memory? Clarify this one day.
- !!!!!!!!!!! DK DK UGLY: The Moho stretching and squishing that Jeroen added to V4.0
- !!!!!!!!!!! DK DK UGLY: should partly deal with this problem.
-
- r_top(1) = R_EARTH
- r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
-
- r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(2) = R80
-
- r_top(3) = R80
- r_bottom(3) = R220
-
- r_top(4) = R220
- r_bottom(4) = R400
-
- r_top(5) = R400
- r_bottom(5) = R600
-
- r_top(6) = R600
- r_bottom(6) = R670
-
- r_top(7) = R670
- r_bottom(7) = R771
-
- r_top(8) = R771
- r_bottom(8) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(9) = RTOPDDOUBLEPRIME
-
- r_top(10) = RTOPDDOUBLEPRIME
- r_bottom(10) = RCMB
-
- r_top(11) = RCMB
- r_bottom(11) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(12) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
-
- r_top(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
- r_bottom(13) = RICB
-
- r_top(14) = RICB
- r_bottom(14) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
-
- rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(2) = R80 / R_EARTH
-
- rmaxs(3) = R80 / R_EARTH
- rmins(3) = R220 / R_EARTH
-
- rmaxs(4) = R220 / R_EARTH
- rmins(4) = R400 / R_EARTH
-
- rmaxs(5) = R400 / R_EARTH
- rmins(5) = R600 / R_EARTH
-
- rmaxs(6) = R600 / R_EARTH
- rmins(6) = R670 / R_EARTH
-
- rmaxs(7) = R670 / R_EARTH
- rmins(7) = R771 / R_EARTH
-
- rmaxs(8:9) = R771 / R_EARTH
- rmins(8:9) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(10) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(10) = RCMB / R_EARTH
-
- rmaxs(11:13) = RCMB / R_EARTH
- rmins(11:13) = RICB / R_EARTH
-
- rmaxs(14) = RICB / R_EARTH
- rmins(14) = R_CENTRAL_CUBE / R_EARTH
-
- else
-
- NUMBER_OF_MESH_LAYERS = 15
- layer_offset = 1
- if ((RMIDDLE_CRUST-RMOHO_FICTITIOUS_IN_MESHER)<(R_EARTH-RMIDDLE_CRUST)) then
- ner( 1) = ceiling (NER_CRUST / 2.d0)
- ner( 2) = floor (NER_CRUST / 2.d0)
- else
- ner( 1) = floor (NER_CRUST / 2.d0)
- ner( 2) = ceiling (NER_CRUST / 2.d0)
- endif
- ner( 3) = NER_80_MOHO
- ner( 4) = NER_220_80
- ner( 5) = NER_400_220
- ner( 6) = NER_600_400
- ner( 7) = NER_670_600
- ner( 8) = NER_771_670
- ner( 9) = NER_TOPDDOUBLEPRIME_771 - elem_doubling_mantle
- ner(10) = elem_doubling_mantle
- ner(11) = NER_CMB_TOPDDOUBLEPRIME
- ner(12) = NER_OUTER_CORE - elem_doubling_middle_outer_core
- ner(13) = elem_doubling_middle_outer_core - elem_doubling_bottom_outer_core
- ner(14) = elem_doubling_bottom_outer_core
- ner(15) = NER_TOP_CENTRAL_CUBE_ICB
-
- ! value of the doubling ratio in each radial region of the mesh
- ratio_sampling_array(1:2) = 1
- ratio_sampling_array(3:9) = 2
- ratio_sampling_array(10:12) = 4
- ratio_sampling_array(13) = 8
- ratio_sampling_array(14:15) = 16
-
- ! value of the doubling index flag in each radial region of the mesh
- doubling_index(1:2) = IFLAG_CRUST
- doubling_index(3) = IFLAG_80_MOHO
- doubling_index(4) = IFLAG_220_80
- doubling_index(5:7) = IFLAG_670_220
- doubling_index(8:11) = IFLAG_MANTLE_NORMAL
- doubling_index(12:14) = IFLAG_OUTER_CORE_NORMAL
- doubling_index(15) = IFLAG_INNER_CORE_NORMAL
-
- ! define the three regions in which we implement a mesh doubling at the top of that region
- this_region_has_a_doubling(:) = .false.
- this_region_has_a_doubling(3) = .true.
- this_region_has_a_doubling(10) = .true.
- this_region_has_a_doubling(13) = .true.
- this_region_has_a_doubling(14) = .true.
- lastdoubling_layer = 14
-
- ! define the top and bottom radii of all the regions of the mesh in the radial direction
- ! the first region is the crust at the surface of the Earth
- ! the last region is in the inner core near the center of the Earth
-
- r_top(1) = R_EARTH
- r_bottom(1) = RMIDDLE_CRUST
-
- r_top(2) = RMIDDLE_CRUST
- r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
-
- r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(3) = R80
-
- r_top(4) = R80
- r_bottom(4) = R220
-
- r_top(5) = R220
- r_bottom(5) = R400
-
- r_top(6) = R400
- r_bottom(6) = R600
-
- r_top(7) = R600
- r_bottom(7) = R670
-
- r_top(8) = R670
- r_bottom(8) = R771
-
- r_top(9) = R771
- r_bottom(9) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
-
- r_top(10) = R_EARTH - DEPTH_SECOND_DOUBLING_REAL
- r_bottom(10) = RTOPDDOUBLEPRIME
-
- r_top(11) = RTOPDDOUBLEPRIME
- r_bottom(11) = RCMB
-
- r_top(12) = RCMB
- r_bottom(12) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
-
- r_top(13) = R_EARTH - DEPTH_THIRD_DOUBLING_REAL
- r_bottom(13) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
-
- r_top(14) = R_EARTH - DEPTH_FOURTH_DOUBLING_REAL
- r_bottom(14) = RICB
-
- r_top(15) = RICB
- r_bottom(15) = R_CENTRAL_CUBE
-
- ! new definition of rmins & rmaxs
- rmaxs(1) = ONE
- rmins(1) = RMIDDLE_CRUST / R_EARTH
-
- rmaxs(2) = RMIDDLE_CRUST / R_EARTH
- rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
-
- rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(3) = R80 / R_EARTH
-
- rmaxs(4) = R80 / R_EARTH
- rmins(4) = R220 / R_EARTH
-
- rmaxs(5) = R220 / R_EARTH
- rmins(5) = R400 / R_EARTH
-
- rmaxs(6) = R400 / R_EARTH
- rmins(6) = R600 / R_EARTH
-
- rmaxs(7) = R600 / R_EARTH
- rmins(7) = R670 / R_EARTH
-
- rmaxs(8) = R670 / R_EARTH
- rmins(8) = R771 / R_EARTH
-
- rmaxs(9:10) = R771 / R_EARTH
- rmins(9:10) = RTOPDDOUBLEPRIME / R_EARTH
-
- rmaxs(11) = RTOPDDOUBLEPRIME / R_EARTH
- rmins(11) = RCMB / R_EARTH
-
- rmaxs(12:14) = RCMB / R_EARTH
- rmins(12:14) = RICB / R_EARTH
-
- rmaxs(15) = RICB / R_EARTH
- rmins(15) = R_CENTRAL_CUBE / R_EARTH
- endif
- endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! calculation of number of elements (NSPEC) below
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ratio_divide_central_cube = maxval(ratio_sampling_array)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 1D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-! theoretical number of spectral elements in radial direction
-do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
- if(iter_region == IREGION_CRUST_MANTLE) then
- ifirst_region = 1
- ilast_region = 10 + layer_offset
- else if(iter_region == IREGION_OUTER_CORE) then
- ifirst_region = 11 + layer_offset
- ilast_region = NUMBER_OF_MESH_LAYERS - 1
- else if(iter_region == IREGION_INNER_CORE) then
- ifirst_region = NUMBER_OF_MESH_LAYERS
- ilast_region = NUMBER_OF_MESH_LAYERS
- else
- stop 'incorrect region code detected'
- endif
- NSPEC1D_RADIAL(iter_region) = sum(ner(ifirst_region:ilast_region))
-enddo
-
-! difference of radial number of element for outer core if the superbrick is cut
- DIFF_NSPEC1D_RADIAL(:,:) = 0
- if (CUT_SUPERBRICK_XI) then
- if (CUT_SUPERBRICK_ETA) then
- DIFF_NSPEC1D_RADIAL(2,1) = 1
- DIFF_NSPEC1D_RADIAL(3,1) = 2
- DIFF_NSPEC1D_RADIAL(4,1) = 1
-
- DIFF_NSPEC1D_RADIAL(1,2) = 1
- DIFF_NSPEC1D_RADIAL(2,2) = 2
- DIFF_NSPEC1D_RADIAL(3,2) = 1
-
- DIFF_NSPEC1D_RADIAL(1,3) = 1
- DIFF_NSPEC1D_RADIAL(3,3) = 1
- DIFF_NSPEC1D_RADIAL(4,3) = 2
-
- DIFF_NSPEC1D_RADIAL(1,4) = 2
- DIFF_NSPEC1D_RADIAL(2,4) = 1
- DIFF_NSPEC1D_RADIAL(4,4) = 1
- else
- DIFF_NSPEC1D_RADIAL(2,1) = 1
- DIFF_NSPEC1D_RADIAL(3,1) = 1
-
- DIFF_NSPEC1D_RADIAL(1,2) = 1
- DIFF_NSPEC1D_RADIAL(4,2) = 1
- endif
- else
- if (CUT_SUPERBRICK_ETA) then
- DIFF_NSPEC1D_RADIAL(3,1) = 1
- DIFF_NSPEC1D_RADIAL(4,1) = 1
-
- DIFF_NSPEC1D_RADIAL(1,2) = 1
- DIFF_NSPEC1D_RADIAL(2,2) = 1
- endif
- endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 2D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! exact number of surface elements for faces along XI and ETA
-
-do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
- if(iter_region == IREGION_CRUST_MANTLE) then
- ifirst_region = 1
- ilast_region = 10 + layer_offset
- else if(iter_region == IREGION_OUTER_CORE) then
- ifirst_region = 11 + layer_offset
- ilast_region = NUMBER_OF_MESH_LAYERS - 1
- else if(iter_region == IREGION_INNER_CORE) then
- ifirst_region = NUMBER_OF_MESH_LAYERS
- ilast_region = NUMBER_OF_MESH_LAYERS
- else
- stop 'incorrect region code detected'
- endif
- tmp_sum_xi = 0
- tmp_sum_eta = 0
- do iter_layer = ifirst_region, ilast_region
- if (this_region_has_a_doubling(iter_layer)) then
- if (ner(iter_layer) == 1) then
- nb_lay_sb = 1
- nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK_1L
- nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK_1L
- else
- nb_lay_sb = 2
- nspec2D_xi_sb = NSPEC2D_XI_SUPERBRICK
- nspec2D_eta_sb = NSPEC2D_ETA_SUPERBRICK
- endif
- doubling = 1
- else
- doubling = 0
- nb_lay_sb = 0
- nspec2D_xi_sb = 0
- nspec2D_eta_sb = 0
- endif
-
- tmp_sum_xi = tmp_sum_xi + ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)) * &
- (ner(iter_layer) - doubling*nb_lay_sb)) + &
- doubling * ((NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)) * (nspec2D_xi_sb/2))
-
- tmp_sum_eta = tmp_sum_eta + ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)) * &
- (ner(iter_layer) - doubling*nb_lay_sb)) + &
- doubling * ((NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)) * (nspec2D_eta_sb/2))
- enddo
- NSPEC2D_XI(iter_region) = tmp_sum_xi
- NSPEC2D_ETA(iter_region) = tmp_sum_eta
- if (iter_region == IREGION_INNER_CORE .and. INCLUDE_CENTRAL_CUBE) then
- NSPEC2D_XI(iter_region) = NSPEC2D_XI(iter_region) + &
- ((NEX_PER_PROC_XI / ratio_divide_central_cube)*(NEX_XI / ratio_divide_central_cube))
- NSPEC2D_ETA(iter_region) = NSPEC2D_ETA(iter_region) + &
- ((NEX_PER_PROC_ETA / ratio_divide_central_cube)*(NEX_XI / ratio_divide_central_cube))
- endif
-enddo
-
-! difference of number of surface elements along xi or eta for outer core if the superbrick is cut
- DIFF_NSPEC2D_XI(:,:) = 0
- DIFF_NSPEC2D_ETA(:,:) = 0
- if (CUT_SUPERBRICK_XI) then
- if (CUT_SUPERBRICK_ETA) then
- DIFF_NSPEC2D_XI(2,1) = 2
- DIFF_NSPEC2D_XI(1,2) = 2
- DIFF_NSPEC2D_XI(2,3) = 2
- DIFF_NSPEC2D_XI(1,4) = 2
-
- DIFF_NSPEC2D_ETA(2,1) = 1
- DIFF_NSPEC2D_ETA(2,2) = 1
- DIFF_NSPEC2D_ETA(1,3) = 1
- DIFF_NSPEC2D_ETA(1,4) = 1
- else
- DIFF_NSPEC2D_ETA(2,1) = 1
- DIFF_NSPEC2D_ETA(1,2) = 1
- endif
- else
- if (CUT_SUPERBRICK_ETA) then
- DIFF_NSPEC2D_XI(2,1) = 2
- DIFF_NSPEC2D_XI(1,2) = 2
- endif
- endif
- DIFF_NSPEC2D_XI(:,:) = DIFF_NSPEC2D_XI(:,:) * (NEX_PER_PROC_XI / ratio_divide_central_cube)
- DIFF_NSPEC2D_ETA(:,:) = DIFF_NSPEC2D_ETA(:,:) * (NEX_PER_PROC_ETA / ratio_divide_central_cube)
-
-! exact number of surface elements on the bottom and top boundaries
-
-! in the crust and mantle
- NSPEC2D_TOP(IREGION_CRUST_MANTLE) = (NEX_XI/ratio_sampling_array(1))*(NEX_ETA/ratio_sampling_array(1))/NPROC
- NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE) = (NEX_XI/ratio_sampling_array(10+layer_offset))*&
- (NEX_ETA/ratio_sampling_array(10+layer_offset))/NPROC
-
-! in the outer core with mesh doubling
- if (ADD_4TH_DOUBLING) then
- NSPEC2D_TOP(IREGION_OUTER_CORE) = (NEX_XI/(ratio_divide_central_cube/4))*(NEX_ETA/(ratio_divide_central_cube/4))/NPROC
- NSPEC2D_BOTTOM(IREGION_OUTER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
- else
- NSPEC2D_TOP(IREGION_OUTER_CORE) = (NEX_XI/(ratio_divide_central_cube/2))*(NEX_ETA/(ratio_divide_central_cube/2))/NPROC
- NSPEC2D_BOTTOM(IREGION_OUTER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
- endif
-
-! in the top of the inner core
- NSPEC2D_TOP(IREGION_INNER_CORE) = (NEX_XI/ratio_divide_central_cube)*(NEX_ETA/ratio_divide_central_cube)/NPROC
- NSPEC2D_BOTTOM(IREGION_INNER_CORE) = NSPEC2D_TOP(IREGION_INNER_CORE)
-
-! maximum number of surface elements on vertical boundaries of the slices
- NSPEC2DMAX_XMIN_XMAX(:) = NSPEC2D_ETA(:)
- NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE) = NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE) + maxval(DIFF_NSPEC2D_ETA(:,:))
- NSPEC2DMAX_YMIN_YMAX(:) = NSPEC2D_XI(:)
- NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE) = NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE) + maxval(DIFF_NSPEC2D_XI(:,:))
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 3D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! exact number of spectral elements in each region
-
-do iter_region = IREGION_CRUST_MANTLE,IREGION_INNER_CORE
- if(iter_region == IREGION_CRUST_MANTLE) then
- ifirst_region = 1
- ilast_region = 10 + layer_offset
- else if(iter_region == IREGION_OUTER_CORE) then
- ifirst_region = 11 + layer_offset
- ilast_region = NUMBER_OF_MESH_LAYERS - 1
- else if(iter_region == IREGION_INNER_CORE) then
- ifirst_region = NUMBER_OF_MESH_LAYERS
- ilast_region = NUMBER_OF_MESH_LAYERS
- else
- stop 'incorrect region code detected'
- endif
- tmp_sum = 0;
- do iter_layer = ifirst_region, ilast_region
- if (this_region_has_a_doubling(iter_layer)) then
- if (ner(iter_layer) == 1) then
- nb_lay_sb = 1
- nspec_sb = NSPEC_SUPERBRICK_1L
- else
- nb_lay_sb = 2
- nspec_sb = NSPEC_DOUBLING_SUPERBRICK
- endif
- doubling = 1
- else
- doubling = 0
- nb_lay_sb = 0
- nspec_sb = 0
- endif
- tmp_sum = tmp_sum + ((NEX_XI / ratio_sampling_array(iter_layer)) * (NEX_ETA / ratio_sampling_array(iter_layer)) * &
- (ner(iter_layer) - doubling*nb_lay_sb)) + &
- doubling * ((NEX_XI / ratio_sampling_array(iter_layer)) * (NEX_ETA / ratio_sampling_array(iter_layer)) * &
- (nspec_sb/4))
- enddo
- NSPEC(iter_region) = tmp_sum / NPROC
-enddo
-
- if(INCLUDE_CENTRAL_CUBE) NSPEC(IREGION_INNER_CORE) = NSPEC(IREGION_INNER_CORE) + &
- (NEX_PER_PROC_XI / ratio_divide_central_cube) * &
- (NEX_PER_PROC_ETA / ratio_divide_central_cube) * &
- (NEX_XI / ratio_divide_central_cube)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! calculation of number of points (NGLOB) below
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 1D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! theoretical number of Gauss-Lobatto points in radial direction
- NGLOB1D_RADIAL(:) = NSPEC1D_RADIAL(:)*(NGLLZ-1)+1
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 2D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! 2-D addressing and buffers for summation between slices
-! we add one to number of points because of the flag after the last point
- NGLOB2DMAX_XMIN_XMAX(:) = NSPEC2DMAX_XMIN_XMAX(:)*NGLLY*NGLLZ + 1
- NGLOB2DMAX_YMIN_YMAX(:) = NSPEC2DMAX_YMIN_YMAX(:)*NGLLX*NGLLZ + 1
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!
-!!!!!! 3D case
-!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! exact number of global points in each region
-
-! initialize array
- NGLOB(:) = 0
-
-! in the inner core (no doubling region + eventually central cube)
- if(INCLUDE_CENTRAL_CUBE) then
- NGLOB(IREGION_INNER_CORE) = ((NEX_PER_PROC_XI/ratio_divide_central_cube) &
- *(NGLLX-1)+1)*((NEX_PER_PROC_ETA/ratio_divide_central_cube) &
- *(NGLLY-1)+1)*((NER_TOP_CENTRAL_CUBE_ICB + NEX_XI / ratio_divide_central_cube)*(NGLLZ-1)+1)
- else
- NGLOB(IREGION_INNER_CORE) = ((NEX_PER_PROC_XI/ratio_divide_central_cube) &
- *(NGLLX-1)+1)*((NEX_PER_PROC_ETA/ratio_divide_central_cube) &
- *(NGLLY-1)+1)*((NER_TOP_CENTRAL_CUBE_ICB)*(NGLLZ-1)+1)
- endif
-
-! in the crust-mantle and outercore
- do iter_region = IREGION_CRUST_MANTLE,IREGION_OUTER_CORE
- if(iter_region == IREGION_CRUST_MANTLE) then
- ifirst_region = 1
- ilast_region = 10 + layer_offset
- else if(iter_region == IREGION_OUTER_CORE) then
- ifirst_region = 11 + layer_offset
- ilast_region = NUMBER_OF_MESH_LAYERS - 1
- else
- stop 'incorrect region code detected'
- endif
- tmp_sum = 0;
- do iter_layer = ifirst_region, ilast_region
- nglob_int_surf_eta=0
- nglob_int_surf_xi=0
- nglob_ext_surf = 0
- nglob_center_edge = 0
- nglob_corner_edge = 0
- nglob_border_edge = 0
- if (this_region_has_a_doubling(iter_layer)) then
- if (iter_region == IREGION_OUTER_CORE .and. iter_layer == lastdoubling_layer .and. &
- (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA)) then
- doubling = 1
- normal_doubling = 0
- cut_doubling = 1
- nb_lay_sb = 2
- nglob_edge = 0
- nglob_surf = 0
- nglob_vol = 8*NGLLX**3 - 12*NGLLX**2 + 6*NGLLX - 1
- nglob_int_surf_eta = 6*NGLLX**2 - 7*NGLLX + 2
- nglob_int_surf_xi = 5*NGLLX**2 - 5*NGLLX + 1
- nglob_ext_surf = 4*NGLLX**2-4*NGLLX+1
- nglob_center_edge = 4*(NGLLX-1)+1
- nglob_corner_edge = 2*(NGLLX-1)+1
- nglob_border_edge = 3*(NGLLX-1)+1
- else
- if (ner(iter_layer) == 1) then
- nb_lay_sb = 1
- nglob_vol = 28*NGLLX**3 - 62*NGLLX**2 + 47*NGLLX - 12
- nglob_surf = 6*NGLLX**2-8*NGLLX+3
- nglob_edge = NGLLX
- else
- nb_lay_sb = 2
- nglob_vol = 32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13
- nglob_surf = 8*NGLLX**2-11*NGLLX+4
- nglob_edge = 2*NGLLX-1
- endif
- doubling = 1
- normal_doubling = 1
- cut_doubling = 0
- endif
- padding = -1
- else
- doubling = 0
- normal_doubling = 0
- cut_doubling = 0
- padding = 0
- nb_lay_sb = 0
- nglob_vol = 0
- nglob_surf = 0
- nglob_edge = 0
- endif
- if (iter_layer == ilast_region) padding = padding +1
- nblocks_xi = NEX_PER_PROC_XI / ratio_sampling_array(iter_layer)
- nblocks_eta = NEX_PER_PROC_ETA / ratio_sampling_array(iter_layer)
-
- tmp_sum = tmp_sum + &
- ((nblocks_xi)*(NGLLX-1)+1) * ((nblocks_eta)*(NGLLX-1)+1) * ((ner(iter_layer) - doubling*nb_lay_sb)*(NGLLX-1)+padding)+&
- normal_doubling * ((((nblocks_xi*nblocks_eta)/4)*nglob_vol) - &
- (((nblocks_eta/2-1)*nblocks_xi/2+(nblocks_xi/2-1)*nblocks_eta/2)*nglob_surf) + &
- ((nblocks_eta/2-1)*(nblocks_xi/2-1)*nglob_edge)) + &
- cut_doubling*(nglob_vol*(nblocks_xi*nblocks_eta) - &
- ( nblocks_eta*(int(nblocks_xi/2)*nglob_int_surf_xi + int((nblocks_xi-1)/2)*nglob_ext_surf) + &
- nblocks_xi*(int(nblocks_eta/2)*nglob_int_surf_eta + int((nblocks_eta-1)/2)*nglob_ext_surf)&
- ) + &
- ( int(nblocks_xi/2)*int(nblocks_eta/2)*nglob_center_edge + &
- int((nblocks_xi-1)/2)*int((nblocks_eta-1)/2)*nglob_corner_edge + &
- ((int(nblocks_eta/2)*int((nblocks_xi-1)/2))+(int((nblocks_eta-1)/2)*int(nblocks_xi/2)))*nglob_border_edge&
- ))
- enddo
- NGLOB(iter_region) = tmp_sum
- enddo
-
-!!! example :
-!!! nblocks_xi/2=5
-!!! ____________________________________
-!!! I I I I I I
-!!! I I I I I I
-!!! I I I I I I
-!!! nblocks_eta/2=3 I______+______+______+______+______I
-!!! I I I I I I
-!!! I I I I I I
-!!! I I I I I I
-!!! I______+______+______+______+______I
-!!! I I I I I I
-!!! I I I I I I
-!!! I I I I I I
-!!! I______I______I______I______I______I
-!!!
-!!! NGLOB for this doubling layer = 3*5*Volume - ((3-1)*5+(5-1)*3)*Surface + (3-1)*(5-1)*Edge
-!!!
-!!! 32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13 -> nb GLL points in a superbrick (Volume)
-!!! 8*NGLLX**2-11*NGLLX+4 -> nb GLL points on a superbrick side (Surface)
-!!! 2*NGLLX-1 -> nb GLL points on a corner edge of a superbrick (Edge)
-
-!!! for the one layer superbrick :
-!!! NGLOB = 28.NGLL^3 - 62.NGLL^2 + 47.NGLL - 12 (Volume)
-!!! NGLOB = 6.NGLL^2 - 8.NGLL + 3 (Surface)
-!!! NGLOB = NGLL (Edge)
-!!!
-!!! those results were obtained by using the script UTILS/doubling_brick/count_nglob_analytical.pl
-!!! with an opendx file of the superbrick's geometry
-
-!!! for the basic doubling bricks (two layers)
-!!! NGLOB = 8.NGLL^3 - 12.NGLL^2 + 6.NGLL - 1 (VOLUME)
-!!! NGLOB = 5.NGLL^2 - 5.NGLL + 1 (SURFACE 1)
-!!! NGLOB = 6.NGLL^2 - 7.NGLL + 2 (SURFACE 2)
-
- end subroutine read_compute_parameters
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/generate_table_NEX_manual.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/generate_table_NEX_manual.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/generate_table_NEX_manual.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -7,10 +7,10 @@
implicit none
- integer NPROC_XI,c,NEX_XI,compteur
+ integer NPROC_XI,c,NEX_XI,counter
- integer, parameter :: NB_COLONNES = 10
- integer, dimension(NB_COLONNES) :: store_NEX_XI
+ integer, parameter :: NB_COLUMNS_TABLE = 10
+ integer, dimension(NB_COLUMNS_TABLE) :: store_NEX_XI
! maximum total number of processors we want to see in the table
integer, parameter :: MAX_NUMBER_OF_PROCS = 100000
@@ -28,28 +28,28 @@
! total number of processors is 6 * NPROC_XI^2
do NPROC_XI = 1,int(sqrt(MAX_NUMBER_OF_PROCS / 6.d0))
- compteur = 1
+ counter = 1
c = 0
- do while (compteur <= NB_COLONNES)
+ do while (counter <= NB_COLUMNS_TABLE)
c = c + 1
NEX_XI = BASE_VALUE * c * NPROC_XI
if(NEX_XI >= 64 .and. mod(NEX_XI,2*BASE_VALUE) == 0) then
- store_NEX_XI(compteur) = NEX_XI
- compteur = compteur + 1
+ store_NEX_XI(counter) = NEX_XI
+ counter = counter + 1
endif
enddo
if(OUTPUT_LATEX_FORMAT) then
write(*,"(i6,' &')") NPROC_XI
write(*,"(i6,' &')") 6*NPROC_XI**2
- do compteur = 1,NB_COLONNES-1
- write(*,"(i6,' &')") store_NEX_XI(compteur)
+ do counter = 1,NB_COLUMNS_TABLE-1
+ write(*,"(i6,' &')") store_NEX_XI(counter)
enddo
- write(*,*) store_NEX_XI(NB_COLONNES),' ',achar(92),'tabularnewline'
+ write(*,*) store_NEX_XI(NB_COLUMNS_TABLE),' ',achar(92),'tabularnewline'
else
write(*,"(i6,' | ',i6,' | ',i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
- NPROC_XI, 6*NPROC_XI**2,(store_NEX_XI(compteur), compteur = 1,NB_COLONNES)
+ NPROC_XI, 6*NPROC_XI**2,(store_NEX_XI(counter), counter = 1,NB_COLUMNS_TABLE)
endif
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -242,7 +242,7 @@
AM_V%Qmu(:) = Mref_V%Qmu_ref(:)
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_JP1D) then
AM_V%Qr(:) = (/ 0.0d0, RICB, RICB, RCMB, RCMB, R670, R670, R220, R220, R120, R120, R_EARTH /)
- AM_V%Qmu(:) = (/ 84.6d0, 84.6d0, 0.0d0, 0.0d0, 312.0d0, 312.0d0, 143.0d0, 143.0d0, 80.0d0, 80.0d0, 600.0d0, 600.0d0 /)
+ AM_V%Qmu(:) = (/ 84.6d0, 84.6d0, 0.0d0, 0.0d0, 312.0d0, 312.0d0, 143.0d0, 143.0d0, 80.0d0, 80.0d0, 600.0d0, 600.0d0 /)
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
AM_V%Qr(:) = SEA1DM_V%radius_sea1d(:)
AM_V%Qmu(:) = SEA1DM_V%Qmu_sea1d(:)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -148,7 +148,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -152,7 +152,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
! get the base pathname for output files
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -141,7 +141,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
print *
print *,'There are ',NPROCTOT,' slices numbered from 0 to ',NPROCTOT-1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -146,7 +146,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -221,7 +221,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
if(.not. SAVE_MESH_FILES) stop 'AVS or DX files were not saved by the mesher'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -229,7 +229,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_header_file.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_header_file.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -140,7 +140,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
! count the total number of sources in the CMTSOLUTION file
call count_number_of_sources(NSOURCES)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_AVS_DX.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_AVS_DX.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -769,7 +769,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
end subroutine read_AVS_DX_parameters
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_GMT_global.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_movie_GMT_global.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -161,7 +161,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
if(.not. MOVIE_SURFACE) stop 'movie frames were not saved by the solver'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -230,7 +230,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
@@ -508,7 +508,7 @@
endif
endif
elseif(THREE_D_MODEL == THREE_D_MODEL_SEA99) then
-! sea99
+! sea99
dvs = ZERO
dvp = ZERO
drho = ZERO
@@ -525,7 +525,7 @@
drho = ZERO
if(theta>=(PI/2.d0 - LAT_MAX*DEGREES_TO_RADIANS) .and. theta<=(PI/2.d0 - LAT_MIN*DEGREES_TO_RADIANS) &
.and. phi>=LON_MIN*DEGREES_TO_RADIANS .and. phi<=LON_MAX*DEGREES_TO_RADIANS) then
- if(r_prem > (R_EARTH - DEP_MAX*1000.d0)/R_EARTH) then
+ if(r_prem > (R_EARTH - DEP_MAX*1000.d0)/R_EARTH) then
call iso3d_dpzhao_model(r,theta,phi,vp,vs,dvp,dvs,rho,found_crust,JP3DM_V)
vpv=vpv*(1.0d0+dvp)
vph=vph*(1.0d0+dvp)
@@ -606,7 +606,7 @@
vsh=vsh*(1.0d0+dvs)
endif
elseif(THREE_D_MODEL == THREE_D_MODEL_SEA99) then
-! sea99
+! sea99
dvs = ZERO
dvp = ZERO
drho = ZERO
@@ -731,7 +731,7 @@
if(THREE_D_MODEL == THREE_D_MODEL_SEA99_JP3D .or. THREE_D_MODEL == THREE_D_MODEL_JP3D) then
if(theta>=(PI/2.d0 - LAT_MAX*DEGREES_TO_RADIANS) .and. theta<=(PI/2.d0 - LAT_MIN*DEGREES_TO_RADIANS) &
.and. phi>=LON_MIN*DEGREES_TO_RADIANS .and. phi<=LON_MAX*DEGREES_TO_RADIANS) then
- if(r > (R_EARTH - DEP_MAX*1000.d0)/R_EARTH) then
+ if(r > (R_EARTH - DEP_MAX*1000.d0)/R_EARTH) then
call iso3d_dpzhao_model(r,theta,phi,vpc,vsc,dvp,dvs,rhoc,found_crust,JP3DM_V)
if(found_crust) then
vpv=vpc
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/jp3d1994_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/jp3d1994_model.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/jp3d1994_model.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -45,13 +45,13 @@
! longitude: 130-145 E
! depth : 0 - 500 km
!
-! Dapeng Zhao
+! Dapeng Zhao
! Dept. of Earth & Planet. Sci
! Washington University
! St. Louis, MO 63130
! U.S.A.
! dapeng at izu.wustl.edu
-!=========================================================================
+!=========================================================================
subroutine read_iso3d_dpzhao_model(JP3DM_V)
implicit none
@@ -136,7 +136,7 @@
!==========================================================================
subroutine iso3d_dpzhao_model(radius,theta,phi,vp,vs,dvp,dvs,rho,found_crust,JP3DM_V)
implicit none
-
+
include "constants.h"
! jp3d_model_variables
@@ -212,9 +212,9 @@
double precision :: PE,RE,HE,H1,H2,H3,scaleval
integer :: LAY
-
+
found_crust = .false.
-
+
PE = theta
RE = phi
HE = (ONE - radius)*R_EARTH_KM
@@ -245,7 +245,7 @@
dvp = 0.01d0*dvp
dvs = 1.5d0*dvp
vp = vp*(1.0d0+dvp)
- vs = vs*(1.0d0+dvs)
+ vs = vs*(1.0d0+dvs)
! determine rho
if(LAY .eq. 1) then
@@ -575,7 +575,7 @@
SUBROUTINE VEL3(PE,RE,HE,V,LAY,JP3DM_V)
implicit none
-
+
include "constants.h"
! jp3d_model_variables
type jp3d_model_variables
@@ -644,11 +644,11 @@
type (jp3d_model_variables) JP3DM_V
! jp3d_model_variables
-
+
double precision :: PE,RE,HE,V
-
+
integer :: LAY
-
+
JP3DM_V%P = 90.0-PE/DEGREES_TO_RADIANS
JP3DM_V%R = RE/DEGREES_TO_RADIANS
JP3DM_V%H = HE
@@ -675,7 +675,7 @@
CALL VABPS(MPB,MRB,MHB,JP3DM_V%VELBP,V,JP3DM_V)
ELSE
END IF
-
+
RETURN
END SUBROUTINE VEL3
@@ -1071,7 +1071,7 @@
IF(HE.LT.HM) THEN
CALL JPMODEL(IPS,HM,VM,JP3DM_V)
V = VM-(HM-HE)*0.003
- ELSE
+ ELSE
CALL JPMODEL(IPS,HE,V,JP3DM_V)
END IF
ELSE
@@ -1160,7 +1160,7 @@
5.227,5.463,5.670,5.850,6.125,6.295,6.395, &
6.483,6.564,6.637,6.706,6.770,6.833,6.893, &
6.953,7.012,7.074,7.137,7.199,7.258,7.314,7.304/
- DATA RA1/1.00,0.99,0.98,0.97,0.96,0.95,0.94,0.93, &
+ DATA RA1/1.00,0.99,0.98,0.97,0.96,0.95,0.94,0.93, &
0.92,0.91,0.90,0.88,0.86,0.84,0.82,0.80, &
0.78,0.76,0.74,0.72,0.70,0.68,0.66,0.64, &
0.62,0.60,0.58,0.56,0.55/
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -666,7 +666,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
if(err_occurred() /= 0) then
call exit_MPI(myrank,'an error occurred while reading the parameter file')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp1d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp1d.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp1d.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -30,9 +30,9 @@
R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST)
implicit none
-
+
include "constants.h"
-
+
! given a normalized radius x, gives the non-dimesionalized density rho,
! speeds vp and vs, and the quality factors Qkappa and Qmu
@@ -103,7 +103,7 @@
!
!--- inner core
!
- if (r >= 0.d0 .and. r <= RICB) then
+ if (r >= 0.d0 .and. r <= RICB) then
rho=13.0885d0-8.8381d0*x*x
vp=11.24094-4.09689*x**2
vs=3.56454-3.45241*x**2
@@ -112,7 +112,7 @@
!
!--- outer core
!
- else if (r > RICB .and. r <= RCMB) then
+ else if (r > RICB .and. r <= RCMB) then
rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
vp=10.03904+3.75665*x-13.67046*x**2
vs=0.0d0
@@ -121,7 +121,7 @@
!
!--- D" at the base of the mantle
!
- else if (r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
+ else if (r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
vp=14.49470-1.47089*x
vs=8.16616-1.58206*x
@@ -130,13 +130,13 @@
!
!--- mantle: from top of D" to d670
!
- else if(r > RTOPDDOUBLEPRIME .and. r <= R771) then
+ else if(r > RTOPDDOUBLEPRIME .and. r <= R771) then
rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
vp=-355.58324*x**4 + 1002.03178*x**3 - 1057.3873425*x**2 + 487.0891011*x - 68.520645
vs=-243.33862*x**4 + 668.06411*x**3 - 685.20113*x**2 + 308.04893*x - 43.737642
Qmu=312.0d0
Qkappa=57827.0d0
- else if(r > R771 .and. r <= R670) then
+ else if(r > R771 .and. r <= R670) then
rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
vp=-174.468866*x**2 + 286.37769*x - 106.034798
vs=-81.0865*x*x + 129.67095*x - 45.268933
@@ -145,42 +145,42 @@
!
!--- mantle: above d670
!
- else if(r > R670 .and. r <= 5871000.d0) then
+ else if(r > R670 .and. r <= 5871000.d0) then
vp=-300.510146*x*x + 511.17372648*x - 206.265832
vs=-139.78275*x*x + 233.3097462*x - 91.0129372
rho=3.3d0 + (vs-4.4d0)*0.7d0
Qmu=143.0d0
Qkappa=57827.0d0
-
- else if(r > 5871000.d0 .and. r <= R400) then
+
+ else if(r > 5871000.d0 .and. r <= R400) then
vp=-601.0202917*x*x + 1063.3823*x - 459.9388738
vs=-145.2465705*x*x + 243.2807524*x - 95.561877
rho=3.3d0 + (vs - 4.4d0)*0.7d0
Qmu=143.0d0
Qkappa=57827.0d0
-
+
else if(r > R400 .and. r <= R220) then
vp=25.042512155*x*x - 68.8367583*x + 51.4120272
vs=15.540158021*x*x - 40.2087657*x + 28.9578929
rho=3.3d0 + (vs - 4.4d0)*0.7d0
Qmu=143.0d0
Qkappa=57827.0d0
-
+
else if(r > R220 .and. r <= R80) then
vp=27.0989608 - 19.473338*x
vs=13.920596 - 9.6309917*x
rho=3.3d0 + (vs - 4.4d0)*0.7d0
Qmu=80.0d0
Qkappa=57827.0d0
-
- else if(r > R80 .and. r <= RMOHO) then
+
+ else if(r > R80 .and. r <= RMOHO) then
vp=26.7663028 - 19.13645*x
vs=13.4601434 - 9.164683*x
rho=3.3d0 + (vs - 4.4d0)*0.7d0
Qmu=600.0d0
Qkappa=57827.0d0
-
- else if(r > RMOHO .and. r <= RMIDDLE_CRUST) then
+
+ else if(r > RMOHO .and. r <= RMIDDLE_CRUST) then
rho=2.9d0
vp = 6.7d0
vs = 3.8d0
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea1d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea1d.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea1d.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -144,7 +144,7 @@
integer i
! define all the values in the model
-
+
SEA1DM_V%radius_sea1d(1)= 0.0000000000
SEA1DM_V%radius_sea1d(2)= 101425.0000000000
SEA1DM_V%radius_sea1d(3)= 202850.0000000000
@@ -152,7 +152,7 @@
SEA1DM_V%radius_sea1d(5)= 405700.0000000000
SEA1DM_V%radius_sea1d(6)= 507125.0000000000
SEA1DM_V%radius_sea1d(7)= 608550.0000000000
- SEA1DM_V%radius_sea1d(8)= 709975.0000000000
+ SEA1DM_V%radius_sea1d(8)= 709975.0000000000
SEA1DM_V%radius_sea1d(9)= 811400.0000000000
SEA1DM_V%radius_sea1d(10)= 912825.0000000000
SEA1DM_V%radius_sea1d(11)= 1014250.0000000000
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -56,7 +56,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,EMULATE_ONLY)
implicit none
@@ -71,7 +71,8 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP, &
+ NEX_XI_read,NEX_ETA_read,NPROC_XI_read,NPROC_ETA_read
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
@@ -88,7 +89,7 @@
ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,WRITE_SEISMOGRAMS_BY_MASTER,&
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,EMULATE_ONLY
character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
@@ -184,15 +185,26 @@
endif
! number of elements at the surface along the two sides of the first chunk
- call read_value_integer(NEX_XI, 'mesher.NEX_XI')
+ call read_value_integer(NEX_XI_read, 'mesher.NEX_XI')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NEX_ETA, 'mesher.NEX_ETA')
+ call read_value_integer(NEX_ETA_read, 'mesher.NEX_ETA')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NPROC_XI, 'mesher.NPROC_XI')
+ call read_value_integer(NPROC_XI_read, 'mesher.NPROC_XI')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
- call read_value_integer(NPROC_ETA, 'mesher.NPROC_ETA')
+ call read_value_integer(NPROC_ETA_read, 'mesher.NPROC_ETA')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ if(.not. EMULATE_ONLY) then
+ NEX_XI = NEX_XI_read
+ NEX_ETA = NEX_ETA_read
+ NPROC_XI = NPROC_XI_read
+ NPROC_ETA = NPROC_ETA_read
+ else
+! this is used in UTILS/estimate_best_values_runs.f90 only, to estimate memory use
+ NEX_ETA = NEX_XI
+ NPROC_ETA = NPROC_XI
+ endif
+
! define the velocity model
call read_value_string(MODEL, 'model.name')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
@@ -734,34 +746,35 @@
call auto_time_stepping(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, DT)
- 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(*,*)'##############################################################'
+!! 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
@@ -939,7 +952,7 @@
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
-! SEA1D without the 2 km of mud layer or the 3km water layer
+! SEA1D without the 2 km of mud layer or the 3km water layer
ROCEAN = 6371000.d0
RMIDDLE_CRUST = 6361000.d0
RMOHO = 6346000.d0
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/sea99_s_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/sea99_s_model.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/sea99_s_model.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -36,7 +36,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
@@ -54,7 +54,7 @@
!----------------------- choose input file: ------------------
! relative anomaly
-
+
open(1,file='DATA/Lebedev_sea99/sea99_dvsvs')
!----------------------- read input file: ------------------
@@ -75,18 +75,18 @@
stop 'alonmin,alonmax,sea99_nlon'
endif
read(1,*)
- do j = 1, SEA99M_V%sea99_ndep
- do ia = 1, SEA99M_V%sea99_nlat
+ do j = 1, SEA99M_V%sea99_ndep
+ do ia = 1, SEA99M_V%sea99_nlat
read (1,*) (SEA99M_V%sea99_vs(ia,io,j), io = 1, SEA99M_V%sea99_nlon)
enddo
enddo
-
+
end subroutine read_sea99_s_model
subroutine sea99_s_model(radius,theta,phi,dvs,SEA99M_V)
-
+
implicit none
-
+
include "constants.h"
! sea99_s_model_variables
@@ -94,7 +94,7 @@
sequence
integer :: sea99_ndep
integer :: sea99_nlat
- integer :: sea99_nlon
+ integer :: sea99_nlon
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
@@ -108,7 +108,7 @@
! sea99_s_model_variables
integer :: id1,i,ilat,ilon
- double precision :: alat1,alon1,radius,theta,phi,dvs
+ double precision :: alat1,alon1,radius,theta,phi,dvs
double precision :: xxx,yyy,dep,pla,plo,xd1,dd1,dd2,ddd(2)
!----------------------- depth in the model ------------------
dep=R_EARTH_KM*(R_UNIT_SPHERE - radius)
@@ -138,7 +138,7 @@
ilon = int((plo - SEA99M_V%alonmin)/SEA99M_V%sea99_ddeg) + 1
alat1 = SEA99M_V%alatmin + (ilat-1)*SEA99M_V%sea99_ddeg
alon1 = SEA99M_V%alonmin + (ilon-1)*SEA99M_V%sea99_ddeg
-
+
do i = 1, 2
xxx = (pla-alat1)/SEA99M_V%sea99_ddeg
yyy = SEA99M_V%sea99_vs(ilat+1,ilon,id1+i-1)-SEA99M_V%sea99_vs(ilat,ilon,id1+i-1)
@@ -151,7 +151,7 @@
enddo
dvs = ddd(1) + (ddd(2)-ddd(1)) * xd1
if(dvs>1.d0) dvs=0.0d0
-
+
end subroutine sea99_s_model
-
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-03-26 15:59:15 UTC (rev 11589)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-03-26 16:02:54 UTC (rev 11590)
@@ -851,7 +851,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
if(err_occurred() /= 0) then
call exit_MPI(myrank,'an error occurred while reading the parameter file')
More information about the cig-commits
mailing list