[cig-commits] r16304 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun Feb 21 21:07:47 PST 2010
Author: danielpeter
Date: 2010-02-21 21:07:47 -0800 (Sun, 21 Feb 2010)
New Revision: 16304
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90
Log:
updated stretching constant to have rmiddle_crust at 15km; updated auto_ner.f90 to use model specific radii (affects high-resolution and/or 1-chunks smaller than 90 degrees simulations)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -140,7 +140,8 @@
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)
+ R_CENTRAL_CUBE, CASE_3D, CRUSTAL, &
+ HONOR_1D_SPHERICAL_MOHO, REFERENCE_1D_MODEL)
implicit none
@@ -152,8 +153,10 @@
NER_670_600, NER_771_670, NER_TOPDDOUBLEPRIME_771, &
NER_CMB_TOPDDOUBLEPRIME, NER_OUTER_CORE, NER_TOP_CENTRAL_CUBE_ICB
double precision R_CENTRAL_CUBE
- logical CASE_3D
+ logical CASE_3D,CRUSTAL,HONOR_1D_SPHERICAL_MOHO
+ integer REFERENCE_1D_MODEL
+ ! local parameters
integer, parameter :: NUM_REGIONS = 14
integer, dimension(NUM_REGIONS) :: scaling
double precision, dimension(NUM_REGIONS) :: radius
@@ -161,23 +164,54 @@
double precision, dimension(NUM_REGIONS-1) :: ratio_bottom
integer, dimension(NUM_REGIONS-1) :: NER
integer NEX_ETA
+ double precision ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ RMOHO_FICTITIOUS_IN_MESHER,R80_FICTITIOUS_IN_MESHER
+ double precision RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS
! This is PREM in Kilometers, well ... kinda, not really ....
- radius(1) = 6371.00d0 ! Surface
- radius(2) = 6346.60d0 ! Moho - 1st Mesh Doubling Interface
- radius(3) = 6291.60d0 ! 80
- radius(4) = 6151.00d0 ! 220
- radius(5) = 5971.00d0 ! 400
- radius(6) = 5771.00d0 ! 600
- radius(7) = 5701.00d0 ! 670
- radius(8) = 5600.00d0 ! 771
- radius(9) = 4712.00d0 ! 1650 - 2nd Mesh Doubling: Geochemical Layering; Kellogg et al. 1999, Science
- radius(10) = 3630.00d0 ! D''
- radius(11) = 3480.00d0 ! CMB
- radius(12) = 2511.00d0 ! 3860 - 3rd Mesh Doubling Interface
- radius(13) = 1371.00d0 ! 5000 - 4th Mesh Doubling Interface
- radius(14) = 982.00d0 ! Top Central Cube
+ !radius(1) = 6371.00d0 ! Surface
+ !radius(2) = 6346.60d0 ! Moho - 1st Mesh Doubling Interface
+ !radius(3) = 6291.60d0 ! 80
+ !radius(4) = 6151.00d0 ! 220
+ !radius(5) = 5971.00d0 ! 400
+ !radius(6) = 5771.00d0 ! 600
+ !radius(7) = 5701.00d0 ! 670
+ !radius(8) = 5600.00d0 ! 771
+ !radius(9) = 4712.00d0 ! 1650 - 2nd Mesh Doubling: Geochemical Layering; Kellogg et al. 1999, Science
+ !radius(10) = 3630.00d0 ! D''
+ !radius(11) = 3480.00d0 ! CMB
+ !radius(12) = 2511.00d0 ! 3860 - 3rd Mesh Doubling Interface
+ !radius(13) = 1371.00d0 ! 5000 - 4th Mesh Doubling Interface
+ !radius(14) = 982.00d0 ! Top Central Cube
+ ! gets model specific radii used to determine number of elements in radial direction
+ call get_model_parameters_radii(REFERENCE_1D_MODEL,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771, &
+ RTOPDDOUBLEPRIME,RCMB,RICB, &
+ RMOHO_FICTITIOUS_IN_MESHER, &
+ R80_FICTITIOUS_IN_MESHER, &
+ RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS, &
+ HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL)
+
+ radius(1) = R_EARTH ! Surface
+ radius(2) = RMOHO_FICTITIOUS_IN_MESHER ! Moho - 1st Mesh Doubling Interface
+ radius(3) = R80 ! 80
+ radius(4) = R220 ! 220
+ radius(5) = R400 ! 400
+ radius(6) = R600 ! 600
+ radius(7) = R670 ! 670
+ radius(8) = R771 ! 771
+ radius(9) = 4712000.0d0 ! 1650 - 2nd Mesh Doubling: Geochemical Layering; Kellogg et al. 1999, Science
+ radius(10) = RTOPDDOUBLEPRIME ! D''
+ radius(11) = RCMB ! CMB
+ radius(12) = 2511000.0d0 ! 3860 - 3rd Mesh Doubling Interface
+ radius(13) = 1371000.0d0 ! 5000 - 4th Mesh Doubling Interface
+ radius(14) = 982000.0d0 ! Top Central Cube
+
+ ! radii in km
+ radius(:) = radius(:) / 1000.0d0
+
call find_r_central_cube(NEX_MAX, radius(14), NEX_ETA)
! Mesh Doubling
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-02-22 05:07:47 UTC (rev 16304)
@@ -510,7 +510,9 @@
!!!!!!!!!!!!!! end of parameters added for the thread-safe version of the code
! for the stretching of crustal elements in the case of 3D models
- double precision, parameter :: MAX_RATIO_CRUST_STRETCHING = 0.6d0
+! (values are chosen for 3D models to have RMOHO_FICTICIOUS at 35 km
+! and RMIDDLE_CRUST to become 15 km with stretching function stretch_tab)
+ double precision, parameter :: MAX_RATIO_CRUST_STRETCHING = 0.75d0
double precision, parameter :: RMOHO_STRETCH_ADJUSTEMENT = 5000.d0
double precision, parameter :: R80_STRETCH_ADJUSTEMENT = -40000.d0
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -418,10 +418,27 @@
enddo
endif
- ! crustal layer stretching: first layer's top and bottom radii will get stretched
+ ! crustal layer stretching: element layer's top and bottom radii will get stretched when in crust
+ ! (number of element layers in crust can vary for different resolutions and 1chunk simulations)
allocate(stretch_tab(2,ner(1)))
if (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH) then
+ ! stretching function determines top and bottom of each element layer in the
+ ! crust region (between r_top(1) and r_bottom(1)), where ner(1) is the
+ ! number of element layers in this crust region
call stretching_function(r_top(1),r_bottom(1),ner(1),stretch_tab)
+
+ ! RMIDDLE_CRUST so far is only used for 1D - models with two layers in the crust
+ ! (i.e. ONE_CRUST is set to .false.), those models do not use CASE_3D
+
+ ! all 3D models use this stretching function to honor a 3D crustal model
+ ! for those models, we set RMIDDLE_CRUST to the bottom of the first element layer
+ ! this value will be used in moho_stretching.f90 to decide whether or not elements
+ ! have to be stretched under oceanic crust.
+ !
+ ! note: stretch_tab uses (dimensionalized) radii from r_top and r_bottom
+ !(with stretch_tab( index_radius(1=top,2=bottom), index_layer( 1=first layer, 2=second layer, 3= ...) )
+ RMIDDLE_CRUST = stretch_tab(2,1)
+
endif
!----
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -960,7 +960,8 @@
implicit none
integer iregion_code
- double precision xmesh,ymesh,zmesh,r
+ ! note: r is the exact radius (and not r_prem with tolerance)
+ double precision xmesh,ymesh,zmesh,r
double precision vpv,vph,vsv,vsh,rho,eta_aniso,dvp
! the 21 coefficients for an anisotropic medium in reduced notation
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -149,10 +149,17 @@
vs = vss(6)
rho = rhos(6)
else if(x > x7 .or. elem_in_crust) then
+ ! takes lower crustal values only if x is slightly above moho depth or
+ ! if elem_in_crust is set
+ !
+ ! note: it looks like this does distinguish between GLL points at the exact moho boundary,
+ ! where the point is on the interface between both,
+ ! oceanic elements and mantle elements below
vp = vps(7)
vs = vss(7)
rho = rhos(7)
else
+ ! note: if x is exactly the moho depth this will return false
found_crust = .false.
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -219,7 +219,8 @@
NER_TOP_CENTRAL_CUBE_ICB,R_CENTRAL_CUBE, &
NEX_MAX,NCHUNKS,REFERENCE_1D_MODEL, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,&
- ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,ANISOTROPIC_INNER_CORE)
+ ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL, &
+ ANISOTROPIC_INNER_CORE)
! 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)
@@ -355,7 +356,8 @@
NER_TOP_CENTRAL_CUBE_ICB,R_CENTRAL_CUBE, &
NEX_MAX,NCHUNKS,REFERENCE_1D_MODEL, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,&
- ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,ANISOTROPIC_INNER_CORE)
+ ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL, &
+ ANISOTROPIC_INNER_CORE)
implicit none
@@ -375,7 +377,7 @@
double precision R_CENTRAL_CUBE
double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES
- logical ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,ANISOTROPIC_INNER_CORE
+ logical ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL,ANISOTROPIC_INNER_CORE
! local variables
integer multiplication_factor
@@ -653,7 +655,8 @@
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)
+ R_CENTRAL_CUBE, CASE_3D, CRUSTAL, &
+ HONOR_1D_SPHERICAL_MOHO, REFERENCE_1D_MODEL)
call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, &
MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD)
@@ -2086,7 +2089,7 @@
(NEX_PER_PROC_ETA / ratio_divide_central_cube) * &
(NEX_XI / ratio_divide_central_cube)
- if(minval(NSPEC) <= 0) stop 'negative NSPEC, there is a problem somewhere :) '
+ if(minval(NSPEC) <= 0) stop 'negative NSPEC, there is a problem somewhere, try to recompile :) '
end subroutine rcp_count_elements
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90 2010-02-22 03:00:20 UTC (rev 16303)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90 2010-02-22 05:07:47 UTC (rev 16304)
@@ -28,6 +28,9 @@
subroutine stretching_function(r_top,r_bottom,ner,stretch_tab)
! define stretch_tab which contains r_top and r_bottom for each element layer in the crust for 3D models.
+!
+! stretch_tab array uses indices index_radius & index_layer :
+! stretch_tab( index_radius (1=top,2=bottom) , index_layer (1=first layer, 2=second layer,..) )
implicit none
@@ -36,16 +39,18 @@
double precision :: r_top, r_bottom,value
integer :: ner,i
double precision, dimension (2,ner) :: stretch_tab
-! for increasing execution speed but have less precision in stretching, increase step
-! not very effective algorithm, but sufficient : used once per proc for meshing.
+ ! for increasing execution speed but have less precision in stretching, increase step
+ ! not very effective algorithm, but sufficient : used once per proc for meshing.
double precision, parameter :: step = 0.001
-! initialize array
+ ! initializes array
+ ! for example: 2 element layers (ner=2) for most probable resolutions (NEX < 1000) in the crust
+ ! then stretch_tab(2,1) = 0.5 = stretch_tab(2,2)
do i=1,ner
stretch_tab(2,i)=(1.d0/ner)
enddo
-! fill with ratio of the layer one thickness for each element
+ ! fill with ratio of the layer one thickness for each element
do while((stretch_tab(2,1) / stretch_tab(2,ner)) > MAX_RATIO_CRUST_STRETCHING)
if (modulo(ner,2) /= 0) then
value = -floor(ner/2.d0)*step
More information about the CIG-COMMITS
mailing list