[cig-commits] r12668 - seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Aug 17 14:52:24 PDT 2008
Author: dkomati1
Date: 2008-08-17 14:52:24 -0700 (Sun, 17 Aug 2008)
New Revision: 12668
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/mantle_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90
Log:
fixed the spreading of the mesh in the crust at very high resolution, which used to create instabilities in the solver
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-17 21:00:46 UTC (rev 12667)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-17 21:52:24 UTC (rev 12668)
@@ -727,7 +727,7 @@
! creating mass matrix in this slice (will be fully assembled in the solver)
if(ipass == 2) rmass(:) = 0._CUSTOM_REAL
- if (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH) then
+ if (.not. PATCH_FOR_GORDON_BELL .and. (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH)) then
allocate(stretch_tab(2,ner(1)),STAT=ier )
if (ier /= 0) then
print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
@@ -781,7 +781,7 @@
! define topological coordinates of this mesh point
offset_x(ignod) = (ix_elem - 1) + iaddx(ignod) * ratio_sampling_array(ilayer)
offset_y(ignod) = (iy_elem - 1) + iaddy(ignod) * ratio_sampling_array(ilayer)
- if (ilayer == 1 .and. CASE_3D) then
+ if (.not. PATCH_FOR_GORDON_BELL .and. (ilayer == 1 .and. CASE_3D)) then
offset_z(ignod) = iaddz(ignod)
else
offset_z(ignod) = (iz_elem - 1) + iaddz(ignod)
@@ -790,10 +790,8 @@
call add_missing_nodes(offset_x,offset_y,offset_z)
! compute the actual position of all the grid points of that element
- if (ilayer == 1 .and. CASE_3D .and. .not. SUPPRESS_CRUSTAL_MESH) then
+ if (.not. PATCH_FOR_GORDON_BELL .and. (ilayer == 1 .and. CASE_3D .and. .not. SUPPRESS_CRUSTAL_MESH)) then
! crustal elements are stretched to be thinner in the upper crust than in lower crust in the 3D case
-! max ratio between size of upper crust elements and lower crust elements is given by the param MAX_RATIO_STRETCHING
-! to avoid stretching, set MAX_RATIO_STRETCHING = 1.d0 in constants.h
call compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &
ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD,iproc_xi,iproc_eta, &
NPROC_XI,NPROC_ETA,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
@@ -1060,7 +1058,7 @@
! end of loop on all the layers of the mesh
enddo
- if (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH) then
+ if (.not. PATCH_FOR_GORDON_BELL .and. (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH)) then
deallocate(stretch_tab,STAT=ier )
if (ier /= 0) then
print *,"ERROR can not deallocate stretch_tab in create_regions_mesh ier=",ier
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/mantle_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/mantle_model.f90 2008-08-17 21:00:46 UTC (rev 12667)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/mantle_model.f90 2008-08-17 21:52:24 UTC (rev 12668)
@@ -192,7 +192,6 @@
type (three_d_mantle_model_variables) D3MM_V
! three_d_mantle_model_variables
-
integer i,j
double precision qqwk(3,NK+1)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-08-17 21:00:46 UTC (rev 12667)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-08-17 21:52:24 UTC (rev 12668)
@@ -494,6 +494,12 @@
! right distribution is determined based upon maximum value of NEX
NEX_MAX = max(NEX_XI,NEX_ETA)
+!! never honor a 1D model in the crust for Gordon Bell runs
+ if(PATCH_FOR_GORDON_BELL) then
+ CASE_3D = .true.
+ HONOR_1D_SPHERICAL_MOHO = .false.
+ endif
+
!----
!---- case prem_onecrust by default
!----
@@ -1003,6 +1009,8 @@
RMOHO_FICTITIOUS_IN_MESHER = RMOHO
else
RMOHO_FICTITIOUS_IN_MESHER = (R80 + R_EARTH) / 2
+ if(PATCH_FOR_GORDON_BELL) &
+ RMOHO_FICTITIOUS_IN_MESHER = R80 + (R_EARTH - R80) * dble(NER_80_MOHO) / dble(NER_CRUST + NER_80_MOHO)
endif
call read_value_double_precision(RECORD_LENGTH_IN_MINUTES, 'solver.RECORD_LENGTH_IN_MINUTES')
@@ -1394,7 +1402,7 @@
rmaxs(14) = RICB / R_EARTH
rmins(14) = R_CENTRAL_CUBE / R_EARTH
- elseif (ONE_CRUST) then
+ else if (ONE_CRUST) then
NUMBER_OF_MESH_LAYERS = 13
layer_offset = 0
@@ -1525,15 +1533,27 @@
else
+!! DK DK this is the case we use for the Gordon Bell runs
+!! DK DK this is the case we use for the Gordon Bell runs
+!! DK DK this is the case we use for the Gordon Bell runs
+
NUMBER_OF_MESH_LAYERS = 14
layer_offset = 1
- if ((RMIDDLE_CRUST-RMOHO_FICTITIOUS_IN_MESHER)<(R_EARTH-RMIDDLE_CRUST)) then
+
+ if ((RMIDDLE_CRUST-RMOHO) < (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
+
+!! DK DK use only one layer in the upper part of the mesh if 3D model
+ if(PATCH_FOR_GORDON_BELL .and. CASE_3D) then
+ ner( 1) = NER_CRUST
+ ner( 2) = 0
+ endif
+
ner( 3) = NER_80_MOHO
ner( 4) = NER_220_80
ner( 5) = NER_400_220
@@ -1584,6 +1604,12 @@
r_top(2) = RMIDDLE_CRUST
r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
+!! DK DK use only one layer in the upper part of the mesh if 3D model
+ if(PATCH_FOR_GORDON_BELL .and. CASE_3D) then
+ r_top(1) = R_EARTH
+ r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
+ endif
+
r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
r_bottom(3) = R80
@@ -1627,6 +1653,12 @@
rmaxs(2) = RMIDDLE_CRUST / R_EARTH
rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+!! DK DK use only one layer in the upper part of the mesh if 3D model
+ if(PATCH_FOR_GORDON_BELL .and. CASE_3D) then
+ rmaxs(1) = ONE
+ rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
+ endif
+
rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
rmins(3) = R80 / R_EARTH
@@ -1658,7 +1690,9 @@
rmins(14) = R_CENTRAL_CUBE / R_EARTH
endif
- else
+
+ else !! DK DK if ADD_4TH_DOUBLING
+
if (SUPPRESS_CRUSTAL_MESH) then
ONE_CRUST = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90 2008-08-17 21:00:46 UTC (rev 12667)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90 2008-08-17 21:52:24 UTC (rev 12668)
@@ -921,7 +921,7 @@
numvar,ivarkern,varstr, &
refmdl,kerstr,hsplfl,dskker,ierror)
else
- write(6,"('the model ',a,' does not exits')") modeldef(1:len_trim(modeldef))
+ write(6,"('model ',a,' does not exits')") modeldef(1:len_trim(modeldef))
endif
! --- check arrays
@@ -995,7 +995,7 @@
end subroutine splcon
-! --- evaluate perturbations in per cent
+! --- evaluate perturbations in percent
subroutine subshsv(xcolat,xlon,xrad,dvsh,dvsv,dvph,dvpv, &
numker,numhpa,numcof,ihpa,lmax,nylm, &
More information about the cig-commits
mailing list