[cig-commits] r16391 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Mar 8 08:36:55 PST 2010
Author: danielpeter
Date: 2010-03-08 08:36:54 -0800 (Mon, 08 Mar 2010)
New Revision: 16391
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL
seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_crustmaps.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_surface.f90
Log:
Sedimentary wavespeeds are superimposed on the mesh if sediment thickness exceeds 2 km; adds regional moho meshing routines for europe and asia; fixes model type declarations
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/INSTALL 2010-03-08 16:36:54 UTC (rev 16391)
@@ -189,3 +189,34 @@
you can either try to use e.g. -mcmodel=medium as additional compiler flag,
or better, change your resolution settings for NPROC_XI,NPROC_ETA and NEX_XI,NEX_ETA.
+4. compilation fails with:
+ ...
+ /opt/ibmcmp/xlf/bg/11.1/bglib/libxl.a(mf2x2.o): In function `mf2x2':
+ mf2x2.c:(.text+0x10): relocation truncated to fit: R_PPC_LOCAL24PC against symbol `_GLOBAL_OFFSET_TABLE_'
+ ...
+
+ this could happen when using IBM xlf90 compiler and gcc to link binaries.
+
+ fix: add the flag -Wl,-relax in the Makefile, e.g.:
+ CFLAGS = -O2 -Wl,-relax
+ should solve the problem.
+
+
+5. compilation fails with:
+ ...
+ obj/libspecfem_mesher.a(read_value_parameters.o): In function `read_value_integer':
+ (.text+0x78): undefined reference to `param_read'
+ ...
+
+ this might happen when the default settings add/suppress trailing underscores
+ to function routine names.
+
+ fix: use corresponding flags for your fortran & C compilers, or remove the underscores at the
+ end of the function names in param_reader.c, like e.g.
+
+ void param_read(char * ..
+
+
+
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -53,8 +53,8 @@
integer numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
integer iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
- !daniel: instead of taking the first value which appears for a global point, average the values
- ! if there are more than one gll points for a global point (points on element corners, edges, faces)
+ ! instead of taking the first value which appears for a global point, average the values
+ ! if there are more than one gll points for a global point (points on element corners, edges, faces)
logical,parameter:: AVERAGE_GLOBALPOINTS = .false.
integer:: ibool_count(NGLOB_CRUST_MANTLE)
real(kind=CUSTOM_REAL):: ibool_dat(NGLOB_CRUST_MANTLE)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -120,8 +120,17 @@
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO ) then
! Stretch mesh to honor smoothed moho thickness from crust2.0
- call moho_stretching_honor_crust(myrank,xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
- R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
+
+ ! differentiate between regional and global meshing
+ if( REGIONAL_MOHO_MESH ) then
+ call moho_stretching_honor_crust_reg(myrank, &
+ xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
+ R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
+ else
+ call moho_stretching_honor_crust(myrank, &
+ xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
+ R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
+ endif
endif
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-03-08 16:36:54 UTC (rev 16391)
@@ -506,9 +506,21 @@
! (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
+ double precision, parameter :: RMOHO_STRETCH_ADJUSTEMENT = 5000.d0 ! moho up to 35km
+ double precision, parameter :: R80_STRETCH_ADJUSTEMENT = -40000.d0 ! r80 down to 120km
+! adapted regional moho stretching
+! 1 chunk simulations, 3-layer crust
+ logical, parameter :: REGIONAL_MOHO_MESH = .false.
+ logical, parameter :: REGIONAL_MOHO_MESH_EUROPE = .false. ! used only for fixing time step
+ logical, parameter :: REGIONAL_MOHO_MESH_ASIA = .false. ! used only for fixing time step
+ logical, parameter :: HONOR_DEEP_MOHO = .false.
+! uncomment for e.g. Europe case, where deep moho is rare
+!! double precision, parameter :: RMOHO_STRETCH_ADJUSTEMENT = -15000.d0 ! moho mesh boundary down to 55km
+! uncomment for deep moho cases, e.g. Asia case (Himalayan moho)
+!! double precision, parameter :: RMOHO_STRETCH_ADJUSTEMENT = -20000.d0 ! moho mesh boundary down to 60km
+
+
! to suppress the crustal layers
! (replaced by an extension of the mantle: R_EARTH is not modified, but no more crustal doubling)
logical, parameter :: SUPPRESS_CRUSTAL_MESH = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -425,8 +425,14 @@
! 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)
-
+
+ ! differentiate between regional meshes or global meshes
+ if( REGIONAL_MOHO_MESH ) then
+ call stretching_function_regional(r_top(1),r_bottom(1),ner(1),stretch_tab)
+ else
+ call stretching_function(r_top(1),r_bottom(1),ner(1),stretch_tab)
+ endif
+
! 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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2010-03-08 16:36:54 UTC (rev 16391)
@@ -165,10 +165,15 @@
# or
#
# CC = gcc
- # CFLAGS = -O3 -m64
+ # CFLAGS = -O3 -m64
#
# for the C compiler when using -q64 for the Fortran compiler
#
+ # on IBM xlf90 compiler:
+ # when encountering errors: ...relocation truncated to fit: R_PPC_LOCAL24PC...
+ # one should also use additional flags:
+ # CFLAGS = -Wl,-relax
+ #
if test x"$FLAGS_NO_CHECK" = x; then
FLAGS_NO_CHECK="-O3 -qsave -qstrict -q64 -qtune=auto -qarch=auto -qcache=auto -qfree=f90 -qsuffix=f=f90 -qhalt=w -qlanglvl=2003pure -qflttrap=overflow:zerodivide:invalid:enable -qsigtrap -qinitauto=7FBFFFFF -Q -Q+rank,swap_all,mxm_m1_m2_5points,mxm_m1_m1_5points,mxm_m2_m1_5points"
# on MareNostrum at the Barcelona SuperComputing Center (Spain) use
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -260,10 +260,12 @@
! EUcrust
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -316,19 +318,23 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
! GLL model_variables
type model_gll_variables
+ sequence
! tomographic iteration model on GLL points
+ double precision :: scale_velocity,scale_density
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
- double precision :: scale_velocity,scale_density
logical :: MODEL_GLL
+ logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
type (model_gll_variables) MGLL_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_crust.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -131,12 +131,12 @@
found_crust = .true.
if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS/R_EARTH_KM) then
+ .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
vp = vps(3)
vs = vss(3)
rho = rhos(3)
else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS/R_EARTH_KM) then
+ .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
vp = vps(4)
vs = vss(4)
rho = rhos(4)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_crustmaps.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_crustmaps.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_crustmaps.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -394,12 +394,12 @@
found_crust = .true.
if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS/R_EARTH_KM) then
+ .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
vp = vps(1)
vs = vss(1)
rho = rhos(1)
else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS/R_EARTH_KM) then
+ .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
vp = vps(2)
vs = vss(2)
rho = rhos(2)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -45,10 +45,12 @@
! EUcrust
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -87,10 +89,12 @@
include "constants.h"
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -154,10 +158,12 @@
implicit none
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -201,10 +207,12 @@
include "constants.h"
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -247,7 +255,13 @@
scaleval = dsqrt(PI*GRAV*RHOAV)
- if( x > x3 ) then
+ if( x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
+ .and. h_basement > MINIMUM_SEDIMENT_THICKNESS) then
+ ! above sediment basement, returns average upper crust value
+ ! since no special sediment values are given
+ found_crust = .true.
+ vp = EUCM_V%eucrust_vp_uppercrust(i+j*ilons) *1000.0d0/(R_EARTH*scaleval)
+ crust_eu = vp
return
else if( x > x4 ) then
found_crust = .true.
@@ -287,10 +301,12 @@
logical :: found
type model_eucrust_variables
+ sequence
double precision, dimension(:),pointer :: eucrust_lat,eucrust_lon,&
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -46,10 +46,12 @@
! GLL model_variables
type model_gll_variables
+ sequence
! tomographic iteration model on GLL points
+ double precision :: scale_velocity,scale_density
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
- double precision :: scale_velocity,scale_density
logical :: MODEL_GLL
+ logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
type (model_gll_variables) MGLL_V
@@ -95,10 +97,12 @@
! GLL model_variables
type model_gll_variables
+ sequence
! tomographic iteration model on GLL points
+ double precision :: scale_velocity,scale_density
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
- double precision :: scale_velocity,scale_density
logical :: MODEL_GLL
+ logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
type (model_gll_variables) MGLL_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -102,10 +102,12 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -151,10 +153,12 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -332,10 +336,12 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -444,10 +450,12 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -545,10 +553,12 @@
! point profile model_variables
type model_ppm_variables
+ sequence
double precision,dimension(:),pointer :: dvs,lat,lon,depth
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
+ integer :: dummy ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -117,17 +117,8 @@
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
- ! offset will be gamma * elevation
- ! scaling cartesian coordinates xyz rather than spherical r/theta/phi involves division of offset by r
- stretch_factor = ONE + gamma * elevation/r
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- xelm(ia) = x * stretch_factor
- yelm(ia) = y * stretch_factor
- zelm(ia) = z * stretch_factor
-
- ! recalculate radius to decide whether this element is in the crust
- r = dsqrt(xelm(ia)*xelm(ia) + yelm(ia)*yelm(ia) + zelm(ia)*zelm(ia))
-
else if ( moho > R_middlecrust ) then
! moho above middle crust
! elements in first layer will squeeze into crust above moho
@@ -153,17 +144,8 @@
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
- ! offset will be gamma * elevation
- ! scaling cartesian coordinates xyz rather than spherical r/theta/phi involves division of offset by r
- stretch_factor = ONE + gamma * elevation/r
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- xelm(ia) = x * stretch_factor
- yelm(ia) = y * stretch_factor
- zelm(ia) = z * stretch_factor
-
- ! recalculate radius to decide whether this element is in the crust
- r = dsqrt(xelm(ia)*xelm(ia) + yelm(ia)*yelm(ia) + zelm(ia)*zelm(ia))
-
end if
! counts corners in above moho
@@ -196,10 +178,423 @@
end subroutine moho_stretching_honor_crust
+
!
+!------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine moho_stretching_honor_crust_reg(myrank, &
+ xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
+ R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
+
+! regional routine: for REGIONAL_MOHO_MESH adaptations
+!
+! uses a 3-layer crust region
+!
+! stretching the moho according to the crust 2.0
+! input: myrank, xelm, yelm, zelm, RMOHO_FICTITIOUS_IN_MESHER R220,RMIDDLE_CRUST, CM_V
+! Dec, 30, 2009
+
+ implicit none
+
+ include "constants.h"
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+ double precision R220,RMIDDLE_CRUST
+ double precision RMOHO_FICTITIOUS_IN_MESHER
+ integer :: myrank
+ logical :: elem_in_crust,elem_in_mantle
+
+ ! local parameters
+ integer:: ia,count_crust,count_mantle
+ double precision:: r,theta,phi,lat,lon
+ double precision:: vpc,vsc,rhoc,moho
+ logical:: found_crust
+
+ double precision, parameter :: RADIANS_TO_DEGREES = 180.d0 / PI
+ double precision, parameter :: PI_OVER_TWO = PI / 2.0d0
+ double precision :: x,y,z
+
+ ! loops over element's anchor points
+ count_crust = 0
+ count_mantle = 0
+ do ia = 1,NGNOD
+
+ ! anchor point location
+ x = xelm(ia)
+ y = yelm(ia)
+ z = zelm(ia)
+
+ call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
+ call reduce(theta,phi)
+
+ lat = 90.d0 - theta * RADIANS_TO_DEGREES
+ lon = phi * RADIANS_TO_DEGREES
+ if( lon > 180.d0 ) lon = lon - 360.0d0
+
+ ! initializes
+ moho = 0.d0
+
+ ! gets smoothed moho depth
+ call meshfem3D_model_crust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
+
+ ! checks moho depth
+ if( abs(moho) < TINYVAL ) call exit_mpi(myrank,'error moho depth to honor')
+
+ moho = ONE - moho
+
+ ! checks if moho will be honored by elements
+ !
+ ! note: we will honor the moho, if the moho depth is
+ ! - above 15km
+ ! - between 25km and 45km
+ ! - below 60 km (in HONOR_DEEP_MOHO case)
+ ! otherwise, the moho will be "interpolated" within the element
+ if( HONOR_DEEP_MOHO) then
+ call stretch_deep_moho(ia,xelm,yelm,zelm,x,y,z,r,moho,R220)
+ else
+ call stretch_moho(ia,xelm,yelm,zelm,x,y,z,r,moho,R220)
+ endif
+
+ ! counts corners in above moho
+ ! note: uses a small tolerance
+ if ( r >= 0.9999d0*moho ) then
+ count_crust = count_crust + 1
+ endif
+ ! counts corners below moho
+ ! again within a small tolerance
+ if ( r <= 1.0001d0*moho ) then
+ count_mantle = count_mantle + 1
+ endif
+
+ end do
+
+ ! sets flag when all corners are above moho
+ if( count_crust == NGNOD) then
+ elem_in_crust = .true.
+ end if
+ ! sets flag when all corners are below moho
+ if( count_mantle == NGNOD) then
+ elem_in_mantle = .true.
+ end if
+
+ ! small stretch check: stretching should affect only points above R220
+ if( r*R_EARTH < R220 ) then
+ print*,'error moho stretching: ',r*R_EARTH,R220,moho*R_EARTH
+ call exit_mpi(myrank,'incorrect moho stretching')
+ endif
+
+ end subroutine moho_stretching_honor_crust_reg
+
+
+
+!
!-------------------------------------------------------------------------------------------------
!
+ subroutine stretch_deep_moho(ia,xelm,yelm,zelm,x,y,z,r,moho,R220)
+
+! honors deep moho (below 60 km), otherwise keeps the mesh boundary at r60 fixed
+
+ implicit none
+
+ include "constants.h"
+
+ integer ia
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision :: x,y,z
+
+ double precision :: r,moho,R220
+
+ ! local parameters
+ double precision :: elevation,gamma
+ ! radii for stretching criteria
+ double precision,parameter :: R15=6356000.d0/R_EARTH
+ double precision,parameter :: R25=6346000.d0/R_EARTH
+ double precision,parameter :: R30=6341000.d0/R_EARTH
+ double precision,parameter :: R35=6336000.d0/R_EARTH
+ double precision,parameter :: R40=6331000.d0/R_EARTH
+ double precision,parameter :: R45=6326000.d0/R_EARTH
+ double precision,parameter :: R50=6321000.d0/R_EARTH
+ double precision,parameter :: R55=6316000.d0/R_EARTH
+ double precision,parameter :: R60=6311000.d0/R_EARTH
+
+ if( RMOHO_STRETCH_ADJUSTEMENT /= -20000.d0 ) &
+ stop 'wrong moho stretch adjustement for stretch_deep_moho'
+
+ if ( moho < R25 .and. moho > R45 ) then
+ ! moho between r25 and r45
+
+ ! stretches mesh at r35 to moho depth
+ elevation = moho - R35
+ if ( r >=R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35))
+ else if ( r < R35 .and. r > R60 ) then
+ gamma = (( r - R60)/( R35 - R60)) ! keeps r60 fixed
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ else if ( moho < R45 ) then
+ ! moho below r45
+
+ ! moves mesh at r35 down to r45
+ elevation = R45 - R35
+ if ( r>= R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35)) ! moves r35 down to r45
+ else if ( r<R35 .and. r>R60 ) then
+ gamma=((r-R60)/(R35-R60)) ! keeps r60 fixed
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ ! add deep moho here
+ if ( moho < R60) then
+ ! moho below r60
+
+ ! stretches mesh at r60 to moho
+ elevation = moho - R60
+ if ( r <R45.and. r >= R60) then
+ gamma=(R45-r)/(R45-R60)
+ else if (r<R60) then
+ gamma=(r-R220/R_EARTH)/(R60-R220/R_EARTH)
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+ end if
+
+ else if (moho > R25) then
+ ! moho above r25
+
+ ! moves mesh at r35 up to r25
+ elevation = R25-R35
+ if (r>=R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35)) ! stretches r35 up to r25
+ else if (r<R35 .and. r>R60 ) then
+ gamma=(r-R60)/(R35-R60) ! keeps r60 fixed
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ ! add shallow moho here
+ if ( moho > R15 ) then
+ ! moho above r15
+
+ ! stretches mesh at r15 to moho depth
+ elevation = moho-R15
+ if (r>=R15) then
+ gamma=(R_UNIT_SPHERE-r)/(R_UNIT_SPHERE-R15)
+ else if (r<R15.and.R>R25) then
+ gamma=(r-R25)/(R15-R25) ! keeps mesh at r25 fixed
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+ end if
+ end if
+
+ end subroutine stretch_deep_moho
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine stretch_moho(ia,xelm,yelm,zelm,x,y,z,r,moho,R220)
+
+! honors shallow and middle depth moho, deep moho will be interpolated within elements
+! mesh will get stretched down to r220
+
+ implicit none
+
+ include "constants.h"
+
+ integer ia
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision :: r,moho,R220
+ double precision :: x,y,z
+
+ ! local parameters
+ double precision :: elevation,gamma
+ ! radii for stretching criteria
+ double precision,parameter :: R15=6356000.d0/R_EARTH
+ double precision,parameter :: R25=6346000.d0/R_EARTH
+ double precision,parameter :: R30=6341000.d0/R_EARTH
+ double precision,parameter :: R35=6336000.d0/R_EARTH
+ double precision,parameter :: R40=6331000.d0/R_EARTH
+ double precision,parameter :: R45=6326000.d0/R_EARTH
+ double precision,parameter :: R50=6321000.d0/R_EARTH
+ double precision,parameter :: R55=6316000.d0/R_EARTH
+ double precision,parameter :: R60=6311000.d0/R_EARTH
+
+ ! moho between 25km and 45 km
+ if ( moho < R25 .and. moho > R45 ) then
+
+ elevation = moho - R35
+ if ( r >=R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35))
+ else if ( r<R35.and.r>R220/R_EARTH) then
+ gamma = ((r-R220/R_EARTH)/(R35-R220/R_EARTH))
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ else if ( moho < R45 ) then
+ ! moho below 45 km
+
+ ! moves mesh at r35 down to r45
+ elevation = R45 - R35
+ if ( r>= R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35))
+ else if ( r<R35.and.r>R220/R_EARTH) then
+ gamma=((r-R220/R_EARTH)/(R35-R220/R_EARTH))
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ else if (moho > R25) then
+ ! moho above 25km
+
+ ! moves mesh at r35 up to r25
+ elevation = R25-R35
+ if (r>=R35.and.r<R15) then
+ gamma=((R15-r)/(R15-R35))
+ else if (r<R35.and.r>R220/R_EARTH) then
+ gamma=(r-R220/R_EARTH)/(R35-R220/R_EARTH)
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+ if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ stop 'incorrect value of gamma for moho from crust 2.0'
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+ ! add shallow moho here
+ if ( moho >R15) then
+ elevation = moho-R15
+ if (r>=R15) then
+ gamma=(R_UNIT_SPHERE-r)/(R_UNIT_SPHERE-R15)
+ else if (r<R15.and.R>R25) then
+ gamma=(r-R25)/(R15-R25)
+ if (abs(gamma)<SMALLVAL) then
+ gamma=0.0d0
+ end if
+ else
+ gamma=0.0d0
+ end if
+
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+ end if
+ endif
+
+ end subroutine stretch_moho
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+
+! moves a point to a new location defined by gamma,elevation and r
+ implicit none
+
+ include "constants.h"
+
+ integer ia
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision :: x,y,z
+
+ double precision :: r,elevation,gamma
+
+ ! local parameters
+ double precision :: stretch_factor
+
+ ! stretch factor
+ ! offset will be gamma * elevation
+ ! scaling cartesian coordinates xyz rather than spherical r/theta/phi involves division of offset by r
+ stretch_factor = ONE + gamma * elevation/r
+
+ ! new point location
+ x = x * stretch_factor
+ y = y * stretch_factor
+ z = z * stretch_factor
+
+ ! stores new point location
+ xelm(ia) = x
+ yelm(ia) = y
+ zelm(ia) = z
+
+ ! new radius
+ r = dsqrt(xelm(ia)*xelm(ia) + yelm(ia)*yelm(ia) + zelm(ia)*zelm(ia))
+
+ end subroutine move_point
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
! obsolete...
!
! subroutine moho_stretching(myrank,xelm,yelm,zelm,RMOHO,R220)
@@ -212,7 +607,8 @@
! integer, parameter :: NL_OCEAN_CONTINENT = 12
!
!! spherical harmonic coefficients of the ocean-continent function (km)
-! double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT),B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
+! double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT), &
+! B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
!
! common /smooth_moho/ A_lm,B_lm
!
@@ -260,7 +656,8 @@
!! stretching between R220 and RMOHO
! gamma = (r - R220/R_EARTH) / (RMOHO/R_EARTH - R220/R_EARTH)
! endif
-! if(gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for Moho topography')
+! if(gamma < -0.0001 .or. gamma > 1.0001) &
+! call exit_MPI(myrank,'incorrect value of gamma for Moho topography')
!
! xelm(ia) = xelm(ia)*(ONE + gamma * elevation / r)
! yelm(ia) = yelm(ia)*(ONE + gamma * elevation / r)
@@ -282,14 +679,16 @@
! integer, parameter :: NL_OCEAN_CONTINENT = 12
!
!! spherical harmonic coefficients of the ocean-continent function (km)
-! double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT),B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
+! double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT), &
+! B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
!
! common /smooth_moho/ A_lm,B_lm
!
!! integer l,m
!!
!! ocean-continent function (km)
-!! open(unit=10,file='DATA/ocean_continent_function/ocean_continent_function.txt',status='old',action='read')
+!! open(unit=10,file='DATA/ocean_continent_function/ocean_continent_function.txt', &
+!! status='old',action='read')
!! do l=0,NL_OCEAN_CONTINENT
!! read(10,*) A_lm(l,0),(A_lm(l,m),B_lm(l,m),m=1,l)
!! enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -742,12 +742,30 @@
! if( THREE_D_MODEL == THREE_D_MODEL_PPM ) &
! DT = DT * (1.d0 - 0.2d0)
-
-
! takes a 5% safety margin on the maximum stable time step
! which was obtained by trial and error
DT = DT * (1.d0 - 0.05d0)
+ ! adapts number of element layers in crust and time step for regional simulations
+ if( REGIONAL_MOHO_MESH ) then
+ ! hard coded number of crustal element layers and time step
+
+ ! checks
+ if( NCHUNKS > 1 ) stop 'regional moho mesh: NCHUNKS error in rcp_set_timestep_and_layers'
+ if( HONOR_1D_SPHERICAL_MOHO ) return
+
+ ! enforce 3 element layers
+ NER_CRUST = 3
+
+ ! increased stability, empirical
+ DT = DT*(1.d0 + 0.5d0)
+
+ if( REGIONAL_MOHO_MESH_EUROPE ) DT = 0.14 ! europe
+ if( REGIONAL_MOHO_MESH_ASIA ) DT = 0.10 ! asia & middle east
+
+ endif
+
+
end subroutine rcp_set_timestep_and_layers
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/stretching_function.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -78,3 +78,72 @@
end subroutine stretching_function
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+subroutine stretching_function_regional(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
+
+ include "constants.h"
+
+ 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.
+! double precision, parameter :: step = 0.001
+!
+! ! 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
+! 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
+! else
+! value = (0.5d0-floor(ner/2.d0))*step
+! endif
+! do i=1,ner
+! stretch_tab(2,i) = stretch_tab(2,i) + value
+! value = value + step
+! enddo
+! enddo
+!
+! deduce r_top and r_bottom
+! ! r_top
+! stretch_tab(1,1) = r_top
+! do i=2,ner
+! stretch_tab(1,i) = sum(stretch_tab(2,i:ner))*(r_top-r_bottom) + r_bottom
+! enddo
+!
+! ! r_bottom
+! stretch_tab(2,ner) = r_bottom
+! do i=1,ner-1
+! stretch_tab(2,i) = stretch_tab(1,i+1)
+! enddo
+
+ if( ner /= 3 ) stop 'error regional stretching function: ner value'
+
+ stretch_tab(1,1) = r_top
+ stretch_tab(1,2) = 6356000.d0 ! 15km second layer top
+ stretch_tab(1,3) = 6336000.d0 ! 35km third layer top
+
+ stretch_tab(2,1) = 6356000.d0 ! bottom first layer
+ stretch_tab(2,2) = 6336000.d0 ! bottom second layer
+ stretch_tab(2,3) = r_bottom ! bottom third layer
+
+end subroutine stretching_function_regional
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_surface.f90 2010-03-05 16:10:39 UTC (rev 16390)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_surface.f90 2010-03-08 16:36:54 UTC (rev 16391)
@@ -76,10 +76,10 @@
do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
ispec = ibelm_top_crust_mantle(ispec2D)
- !daniel: in case of global, NCHUNKS_VAL == 6 simulations, be aware that for
- ! the cubed sphere, the mapping changes for different chunks,
- ! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates.
- ! for future consideration, like in create_movie_GMT_global.f90 ...
+ ! in case of global, NCHUNKS_VAL == 6 simulations, be aware that for
+ ! the cubed sphere, the mapping changes for different chunks,
+ ! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates.
+ ! for future consideration, like in create_movie_GMT_global.f90 ...
k = NGLLZ
! loop on all the points inside the element
More information about the CIG-COMMITS
mailing list