[cig-commits] [commit] devel, master: Added initial comprehensive earth model functionality to mesher (0a011ec)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:29:42 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 0a011ecccf823e34325aaa711bc5e824bdfb6652
Author: Michael Afanasiev <michael.afanasiev at erdw.ethz.ch>
Date: Wed Sep 3 14:43:17 2014 +0200
Added initial comprehensive earth model functionality to mesher
>---------------------------------------------------------------
0a011ecccf823e34325aaa711bc5e824bdfb6652
DATA/Par_file | 2 +-
DATA/cemRequest/readme.txt | 10 +++++++++
src/meshfem3D/create_meshes.f90 | 4 ++++
src/meshfem3D/create_regions_mesh.F90 | 21 ++++++++++++++++++-
src/meshfem3D/get_model.f90 | 3 ++-
.../{meshfem3D_models.f90 => meshfem3D_models.F90} | 21 ++++++++++++++++++-
src/meshfem3D/meshfem3D_par.f90 | 3 ++-
src/meshfem3D/model_cem.f90 | 2 --
src/meshfem3D/rules.mk | 2 --
src/shared/broadcast_computed_parameters.f90 | 7 +++++--
src/shared/get_model_parameters.f90 | 24 +++++++++++++++++-----
src/shared/read_compute_parameters.f90 | 3 ++-
src/shared/shared_par.f90 | 3 ++-
13 files changed, 87 insertions(+), 18 deletions(-)
diff --git a/DATA/Par_file b/DATA/Par_file
index 789b043..b413b75 100644
--- a/DATA/Par_file
+++ b/DATA/Par_file
@@ -38,7 +38,7 @@ NPROC_ETA = 2
# to take the 1D crustal model from the
# associated reference model rather than the default 3D crustal model
# e.g. s20rts_1Dcrust, s362ani_1Dcrust, etc.
-MODEL = 1D_transversely_isotropic_prem
+MODEL = CEM_REQUEST
# parameters describing the Earth model
OCEANS = .true.
diff --git a/DATA/cemRequest/readme.txt b/DATA/cemRequest/readme.txt
new file mode 100644
index 0000000..ede09eb
--- /dev/null
+++ b/DATA/cemRequest/readme.txt
@@ -0,0 +1,10 @@
+--- README ---
+
+COMPREHENSIVE EARTH MODEL [CEM]
+This is an empty directory that holds CEM request files, when the CEM_ACCEPT or CEM_REQUEST flag is specified in the Par_file.
+
+See: Towards a comprehensive earth model across the scales
+ Michael Afanasiev, Andreas Fichtner, and Daniel Peter
+ Geophysical Research Abstracts
+ Vol. 16, EGU2014-11244, 2014
+ EGU General Assembly 2014
diff --git a/src/meshfem3D/create_meshes.f90 b/src/meshfem3D/create_meshes.f90
index f10b02d..31641f2 100644
--- a/src/meshfem3D/create_meshes.f90
+++ b/src/meshfem3D/create_meshes.f90
@@ -115,6 +115,10 @@
NSPEC2D_BOTTOM(iregion_code),NSPEC2D_TOP(iregion_code), &
mod(iproc_xi_slice(myrank),2),mod(iproc_eta_slice(myrank),2), &
ipass)
+
+ ! If we're in the request stage of CEM, exit.
+ if (CEM_REQUEST) exit
+
enddo
! deallocate arrays used for that region
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index b648172..d6f7f26 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -60,7 +60,7 @@
use meshfem3D_models_par,only: &
SAVE_BOUNDARY_MESH,SUPPRESS_CRUSTAL_MESH,REGIONAL_MOHO_MESH, &
- OCEANS
+ OCEANS,CEM_REQUEST
use create_MPI_interfaces_par, only: &
NGLOB1D_RADIAL_MAX,iboolcorner,iboolfaces, &
@@ -185,6 +185,25 @@
njmin,njmax, nkmin_xi,nkmin_eta, NSPEC2DMAX_XMIN_XMAX, &
NSPEC2DMAX_YMIN_YMAX, NSPEC2D_BOTTOM)
+ ! Only go into here if we're requesting xyz files for CEM
+#if defined (CEM)
+ if (CEM_REQUEST) then
+
+ call build_global_coordinates (nspec, nglob_theor, iregion_code)
+ call write_cem_request (iregion_code)
+ call synchronize_all ( )
+
+ deallocate(ibool1D_leftxi_lefteta,ibool1D_rightxi_lefteta, &
+ ibool1D_leftxi_righteta,ibool1D_rightxi_righteta, &
+ xyz1D_leftxi_lefteta,xyz1D_rightxi_lefteta, &
+ xyz1D_leftxi_righteta,xyz1D_rightxi_righteta,iboolleft_xi, &
+ iboolright_xi,iboolleft_eta,iboolright_eta,nimin,nimax, &
+ njmin, njmax,nkmin_xi,nkmin_eta,iboolfaces,iboolcorner)
+
+ end if
+#endif
+
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
case( 2 ) !!!!!!!!!!! second pass of the mesher
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
diff --git a/src/meshfem3D/get_model.f90 b/src/meshfem3D/get_model.f90
index fc8ff77..db391b6 100644
--- a/src/meshfem3D/get_model.f90
+++ b/src/meshfem3D/get_model.f90
@@ -158,7 +158,8 @@
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
! gets the 3-D crustal model
- if (CRUSTAL) then
+ ! M.A. don't overwrite crust if using CEM.
+ if (CRUSTAL .and. .not. CEM_ACCEPT) then
if (.not. elem_in_mantle) &
call meshfem3D_models_get3Dcrust_val(iregion_code,xmesh,ymesh,zmesh,r, &
vpv,vph,vsv,vsh,rho,eta_aniso,dvp, &
diff --git a/src/meshfem3D/meshfem3D_models.f90 b/src/meshfem3D/meshfem3D_models.F90
similarity index 98%
rename from src/meshfem3D/meshfem3D_models.f90
rename to src/meshfem3D/meshfem3D_models.F90
index 2b6e1c0..14fa012 100644
--- a/src/meshfem3D/meshfem3D_models.f90
+++ b/src/meshfem3D/meshfem3D_models.F90
@@ -142,6 +142,13 @@
if (ANISOTROPIC_3D_MANTLE) &
call model_aniso_mantle_broadcast(myrank)
+ ! Enclose this in an ifdef so we don't link to netcdf
+ ! if we don't need it.
+#if defined (CEM)
+ if (CEM_REQUEST .or. CEM_ACCEPT) &
+ call model_cem_broadcast(myrank)
+#endif
+
! crustal model
if (CRUSTAL) &
call meshfem3D_crust_broadcast(myrank)
@@ -343,12 +350,15 @@
RCMB,R670,RMOHO, &
xmesh,ymesh,zmesh,r, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
- c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
+ ispec,i,j,k)
use meshfem3D_models_par
+
implicit none
+ intent (in) :: ispec, i, j, k
integer iregion_code
double precision r_prem
double precision rho,dvp
@@ -367,6 +377,9 @@
real(kind=4) :: xcolat,xlon,xrad,dvpv,dvph,dvsv,dvsh
logical :: found_crust,suppress_mantle_extension
+ ! CEM needs these to determine iglob
+ integer :: ispec, i, j, k
+
! initializes perturbation values
dvs = ZERO
dvp = ZERO
@@ -595,6 +608,12 @@
endif
endif ! ANISOTROPIC_3D_MANTLE
+#if defined (CEM)
+ if (CEM_ACCEPT) then
+ call request_cem (vsh, vsv, vph, vpv, rho, iregion_code, ispec, i, j, k)
+ end if
+#endif
+
!> Hejun
! Assign Attenuation after get 3-D crustal model
! This is here to identify how and where to include 3D attenuation
diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90
index 3b33eeb..2653640 100644
--- a/src/meshfem3D/meshfem3D_par.f90
+++ b/src/meshfem3D/meshfem3D_par.f90
@@ -50,7 +50,8 @@
HONOR_1D_SPHERICAL_MOHO,CRUSTAL,ONE_CRUST,CASE_3D,TRANSVERSE_ISOTROPY, &
ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE, &
ATTENUATION_3D, &
- ANISOTROPIC_INNER_CORE
+ ANISOTROPIC_INNER_CORE, &
+ CEM_REQUEST, CEM_ACCEPT
implicit none
diff --git a/src/meshfem3D/model_cem.f90 b/src/meshfem3D/model_cem.f90
index c4be778..7d6c363 100644
--- a/src/meshfem3D/model_cem.f90
+++ b/src/meshfem3D/model_cem.f90
@@ -8,7 +8,6 @@ module cem_par
double precision, parameter :: scale_GPa = (RHOAV / 1000.d0) * &
((R_EARTH * scaleval / 1000.d0) ** 2)
-
integer, dimension (:), allocatable :: regCode
integer, parameter :: shuOn=1
integer, parameter :: comLvl=9
@@ -67,7 +66,6 @@ subroutine model_cem_broadcast ( myrank )
call return_populated_arrays (reg3Bc, "vpp", 3)
call synchronize_all ()
- if (myrank == 0 ) print *, "Finished reading in netcdf."
end if
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index c0e4a59..11163a7 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -230,8 +230,6 @@ endif
# conditional CEM model
ifeq ($(CEM),yes)
meshfem3D_OBJECTS += $O/model_cem.checknetcdf.o
-NETCDF_INCLUDE = -I$(netcdf_specfem_dir)/include -I$(hdf5_specfem_dir)/include -I$(zlib_specfem_dir)/include
-LDFLAGS += -L$(netcdf_specfem_dir)/lib -lnetcdff -lnetcdf -L$(hdf5_specfem_dir)/lib -lhdf5_hl -lhdf5 -L$(zlib_specfem_dir)/lib -lz
endif
diff --git a/src/shared/broadcast_computed_parameters.f90 b/src/shared/broadcast_computed_parameters.f90
index 45f8c0d..b57ad41 100644
--- a/src/shared/broadcast_computed_parameters.f90
+++ b/src/shared/broadcast_computed_parameters.f90
@@ -39,7 +39,7 @@
integer, parameter :: nparam_i = 44
integer, dimension(nparam_i) :: bcast_integer
- integer, parameter :: nparam_l = 57
+ integer, parameter :: nparam_l = 59
logical, dimension(nparam_l) :: bcast_logical
integer, parameter :: nparam_dp = 34
@@ -99,7 +99,8 @@
ADIOS_ENABLED,ADIOS_FOR_FORWARD_ARRAYS, &
ADIOS_FOR_MPI_ARRAYS,ADIOS_FOR_ARRAYS_SOLVER, &
ADIOS_FOR_SOLVER_MESHFILES,ADIOS_FOR_AVS_DX,&
- ADIOS_FOR_KERNELS,ADIOS_FOR_MODELS, ADIOS_FOR_UNDO_ATTENUATION &
+ ADIOS_FOR_KERNELS,ADIOS_FOR_MODELS, ADIOS_FOR_UNDO_ATTENUATION, &
+ CEM_REQUEST, CEM_ACCEPT &
/)
bcast_double_precision = (/ &
@@ -265,6 +266,8 @@
ADIOS_FOR_KERNELS = bcast_logical(55)
ADIOS_FOR_MODELS = bcast_logical(56)
ADIOS_FOR_UNDO_ATTENUATION = bcast_logical(57)
+ CEM_REQUEST = bcast_logical(58)
+ CEM_ACCEPT = bcast_logical(59)
! double precisions
DT = bcast_double_precision(1)
diff --git a/src/shared/get_model_parameters.f90 b/src/shared/get_model_parameters.f90
index 5a9d255..987e0e9 100644
--- a/src/shared/get_model_parameters.f90
+++ b/src/shared/get_model_parameters.f90
@@ -32,7 +32,8 @@
OCEANS,TOPOGRAPHY, &
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)
+ R80_FICTITIOUS_IN_MESHER,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS, &
+ CEM_REQUEST, CEM_ACCEPT)
use constants
@@ -44,7 +45,7 @@
logical ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ATTENUATION_3D, &
CASE_3D,CRUSTAL,HETEROGEN_3D_MANTLE,HONOR_1D_SPHERICAL_MOHO,&
- ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY
+ ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY,CEM_REQUEST,CEM_ACCEPT
logical OCEANS,TOPOGRAPHY
@@ -59,7 +60,7 @@
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ATTENUATION_3D, &
CASE_3D,CRUSTAL,HETEROGEN_3D_MANTLE,HONOR_1D_SPHERICAL_MOHO, &
ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY, &
- OCEANS,TOPOGRAPHY)
+ OCEANS,TOPOGRAPHY,CEM_REQUEST,CEM_ACCEPT)
! sets radius for each discontinuity and ocean density values
call get_model_parameters_radii(REFERENCE_1D_MODEL,ROCEAN,RMIDDLE_CRUST, &
@@ -83,7 +84,7 @@
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ATTENUATION_3D, &
CASE_3D,CRUSTAL,HETEROGEN_3D_MANTLE,HONOR_1D_SPHERICAL_MOHO, &
ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY, &
- OCEANS,TOPOGRAPHY)
+ OCEANS,TOPOGRAPHY,CEM_REQUEST,CEM_ACCEPT)
use constants
@@ -95,7 +96,8 @@
logical ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ATTENUATION_3D, &
CASE_3D,CRUSTAL,HETEROGEN_3D_MANTLE,HONOR_1D_SPHERICAL_MOHO,&
- ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY
+ ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY,&
+ CEM_REQUEST,CEM_ACCEPT
logical OCEANS,TOPOGRAPHY
! local parameters
@@ -177,6 +179,10 @@
ISOTROPIC_3D_MANTLE = .false.
HONOR_1D_SPHERICAL_MOHO = .false.
+ ! no CEM by default
+ CEM_REQUEST = .false.
+ CEM_ACCEPT = .false.
+
! no 3D model by default
THREE_D_MODEL = 0
TRANSVERSE_ISOTROPY = .false.
@@ -375,6 +381,14 @@
THREE_D_MODEL = THREE_D_MODEL_S362ANI
TRANSVERSE_ISOTROPY = .true.
+ else if (MODEL_ROOT == 'CEM_REQUEST') then
+ CEM_REQUEST = .true.
+ TRANSVERSE_ISOTROPY = .true.
+
+ else if (MODEL_ROOT == 'CEM_ACCEPT') then
+ CEM_ACCEPT = .true.
+ TRANSVERSE_ISOTROPY = .true.
+
else if (MODEL_ROOT == 'PPM') then
! superimposed based on isotropic-prem
CASE_3D = .true.
diff --git a/src/shared/read_compute_parameters.f90 b/src/shared/read_compute_parameters.f90
index 4014ff7..8b6257a 100644
--- a/src/shared/read_compute_parameters.f90
+++ b/src/shared/read_compute_parameters.f90
@@ -100,7 +100,8 @@
OCEANS,TOPOGRAPHY, &
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)
+ R80_FICTITIOUS_IN_MESHER,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS, &
+ CEM_REQUEST,CEM_ACCEPT)
! sets time step size and number of layers
! right distribution is determined based upon maximum value of NEX
diff --git a/src/shared/shared_par.f90 b/src/shared/shared_par.f90
index 3cf22b5..f1c8711 100644
--- a/src/shared/shared_par.f90
+++ b/src/shared/shared_par.f90
@@ -165,7 +165,8 @@
! model flags
integer :: REFERENCE_1D_MODEL,THREE_D_MODEL
logical :: TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- CRUSTAL,ONE_CRUST,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE
+ CRUSTAL,ONE_CRUST,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE, &
+ CEM_REQUEST,CEM_ACCEPT
logical :: ATTENUATION_3D
logical :: INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE
logical :: EMULATE_ONLY = .false.
More information about the CIG-COMMITS
mailing list