[cig-commits] [commit] devel: - added a routine to compute the total mass of the Earth - fixed double precision compilation support, which was broken - changed USE_VERSION_5_1_5 to USE_OLD_VERSION_5_1_5_FORMAT (in the process it seems my script also removed all white spaces at the end of lines, if any) - added support for ROLAND_SYLVAIN simulations (86b5dab)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Apr 29 16:45:30 PDT 2014
Repository : ssh://geoshell/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/73e65c6efc9d20cf5318d73e3a5b9d5c16b817df...bd213f2df89745397090a39fef05f3632538204d
>---------------------------------------------------------------
commit 86b5dab8841a73c51583185724ce337ff37ffd26
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed Apr 30 01:38:50 2014 +0200
- added a routine to compute the total mass of the Earth
- fixed double precision compilation support, which was broken
- changed USE_VERSION_5_1_5 to USE_OLD_VERSION_5_1_5_FORMAT (in the process it seems my script also removed all white spaces at the end of lines, if any)
- added support for ROLAND_SYLVAIN simulations
>---------------------------------------------------------------
86b5dab8841a73c51583185724ce337ff37ffd26
setup/constants.h.in | 13 +++---
src/auxiliaries/combine_vol_data.F90 | 3 +-
src/cuda/assemble_MPI_scalar_cuda.cu | 0
src/cuda/assemble_MPI_vector_cuda.cu | 0
src/cuda/check_fields_cuda.cu | 0
src/cuda/compute_add_sources_elastic_cuda.cu | 0
src/cuda/compute_coupling_cuda.cu | 0
src/cuda/compute_forces_crust_mantle_cuda.cu | 0
src/cuda/compute_forces_inner_core_cuda.cu | 0
src/cuda/compute_forces_outer_core_cuda.cu | 0
src/cuda/compute_kernels_cuda.cu | 0
src/cuda/compute_stacey_acoustic_cuda.cu | 0
src/cuda/compute_stacey_elastic_cuda.cu | 0
src/cuda/initialize_cuda.cu | 0
src/cuda/noise_tomography_cuda.cu | 0
src/cuda/prepare_mesh_constants_cuda.cu | 0
src/cuda/save_and_compare_cpu_vs_gpu.c | 0
src/cuda/specfem3D_gpu_cuda_method_stubs.c | 0
src/cuda/transfer_fields_cuda.cu | 0
src/cuda/update_displacement_cuda.cu | 0
src/cuda/write_seismograms_cuda.cu | 0
src/meshfem3D/add_topography_410_650.f90 | 14 +-----
src/meshfem3D/compute_area.f90 | 29 +++++++++++-
src/meshfem3D/compute_element_properties.f90 | 2 +-
src/meshfem3D/compute_volumes.f90 | 69 +++++++++++++++++++++++++++-
src/meshfem3D/create_mass_matrices.f90 | 8 ++--
src/meshfem3D/create_regions_mesh.F90 | 20 ++++++--
src/meshfem3D/finalize_mesher.f90 | 12 ++++-
src/meshfem3D/initialize_mesher.f90 | 5 ++
src/meshfem3D/meshfem3D_par.f90 | 3 ++
src/meshfem3D/moho_stretching.f90 | 10 ++--
src/shared/adios_method_stubs.c | 0
src/shared/force_ftz.c | 0
src/shared/parallel.f90 | 24 ++++++++++
src/shared/param_reader.c | 0
src/shared/read_parameter_file.f90 | 2 +-
src/shared/rthetaphi_xyz.f90 | 9 ++--
src/shared/write_c_binary.c | 0
src/specfem3D/asdf_data.f90 | 0
src/specfem3D/file_io_threads.c | 0
src/specfem3D/iterate_time.F90 | 6 +--
src/specfem3D/prepare_timerun.f90 | 31 ++++++-------
src/specfem3D/specfem3D_par.F90 | 6 +--
src/specfem3D/visual_vtk_stubs.c | 0
44 files changed, 201 insertions(+), 65 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 210d8c4..67c7d07 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -48,6 +48,9 @@
!*********************************************************************************************************
+!! DK DK for Roland_Sylvain
+ logical, parameter :: ROLAND_SYLVAIN = .false.
+
! if files on a local path on each node are also seen as global with same path
! set to .true. typically on a machine with a common (shared) file system, e.g. LUSTRE, GPFS or NFS-mounted /home
! set to .false. typically on a cluster of nodes with local disks used to store the mesh (not very common these days)
@@ -128,7 +131,7 @@
!! pathname of the topography file (un-smoothed)
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/ETOPO1.xyz'
-! Use GLL points to capture TOPOGRAPHY and ELLIPTICITY (experimental feature)
+! Use GLL points to capture TOPOGRAPHY and ELLIPTICITY (experimental feature, currently does not work, DO NOT USE)
logical,parameter :: USE_GLL = .false.
! maximum depth of the oceans in trenches and height of topo in mountains
@@ -306,13 +309,9 @@
!! new version compatibility
!!
!!-----------------------------------------------------------
-!! (uncomment desired resolution)
-! old version
! old version 5.1.5 uses full 3d attenuation arrays (set to .true.), custom_real for attenuation arrays (Qmu_store, tau_e_store)
-! logical, parameter :: USE_VERSION_5_1_5 = .true.
-! new version
! new version uses optional full 3d attenuation
- logical, parameter :: USE_VERSION_5_1_5 = .false.
+ logical, parameter :: USE_OLD_VERSION_5_1_5_FORMAT = .false.
!!-----------------------------------------------------------
@@ -406,7 +405,7 @@
integer, parameter :: MIDY = (NGLLY+1)/2
integer, parameter :: MIDZ = (NGLLZ+1)/2
-! gravitational constant
+! gravitational constant in S.I. units i.e. in m3 kg-1 s-2
double precision, parameter :: GRAV = 6.6723d-11
! a few useful constants
diff --git a/src/auxiliaries/combine_vol_data.F90 b/src/auxiliaries/combine_vol_data.F90
index be9847d..8118c11 100644
--- a/src/auxiliaries/combine_vol_data.F90
+++ b/src/auxiliaries/combine_vol_data.F90
@@ -61,7 +61,8 @@ program combine_vol_data_vtk
integer num_ibool(NGLOB_CRUST_MANTLE)
logical mask_ibool(NGLOB_CRUST_MANTLE), HIGH_RESOLUTION_MESH
- real x, y, z, dat
+ real(kind=CUSTOM_REAL) x, y, z
+ real dat
integer numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
integer iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
! instead of taking the first value which appears for a global point, average the values
diff --git a/src/meshfem3D/add_topography_410_650.f90 b/src/meshfem3D/add_topography_410_650.f90
index 0bc6e84..efb0430 100644
--- a/src/meshfem3D/add_topography_410_650.f90
+++ b/src/meshfem3D/add_topography_410_650.f90
@@ -51,11 +51,6 @@
! note: adding topography to 410 and 660 strongly affects PcP, PKiKP, etc. phases,
! we leave it in and check whether the stretching makes simulation unstable
-! if( .not. USE_VERSION_5_1_5 ) then
-! !! DK DK added this safety test for now
-! !stop 'there is a currently a bug in this routine, it makes the mesher crash'
-! endif
-
! we loop on all the points of the element
do ia = 1,NGNOD
@@ -63,7 +58,7 @@
y = yelm(ia)
z = zelm(ia)
- if( USE_VERSION_5_1_5) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT) then
! convert to r theta phi
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
call reduce(theta,phi)
@@ -153,11 +148,6 @@
! note: adding topography to 410 and 660 strongly affects PcP, PKiKP, etc. phases,
! we leave it in and check whether the stretching makes simulation unstable
-! if( .not. USE_VERSION_5_1_5 ) then
-! !! DK DK added this safety test for now
-! !stop 'there is a currently a bug in this routine, it makes the mesher crash'
-! endif
-
! we loop on all GLL points of the element
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -167,7 +157,7 @@
y = ystore(i,j,k,ispec)
z = zstore(i,j,k,ispec)
- if( USE_VERSION_5_1_5) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT) then
! convert to r theta phi
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
call reduce(theta,phi)
diff --git a/src/meshfem3D/compute_area.f90 b/src/meshfem3D/compute_area.f90
index a1eb286..4605b4b 100644
--- a/src/meshfem3D/compute_area.f90
+++ b/src/meshfem3D/compute_area.f90
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine compute_area(myrank,NCHUNKS,iregion_code, &
area_local_bottom,area_local_top,&
volume_local,volume_total, &
@@ -95,3 +94,31 @@
endif
end subroutine compute_area
+
+!=====================================================================
+
+ ! computes total Earth mass
+
+ subroutine compute_total_Earth_mass(myrank,Earth_mass_local,Earth_mass_total)
+
+ use meshfem3D_models_par
+
+ implicit none
+
+ integer :: myrank
+
+ double precision :: Earth_mass_local
+ double precision :: Earth_mass_total
+
+ ! local parameters
+ double precision :: Earth_mass_total_region
+
+ ! use an MPI reduction to compute the total Earth mass
+ Earth_mass_total_region = ZERO
+ call sum_all_dp(Earth_mass_local,Earth_mass_total_region)
+
+ ! sum volume over all the regions
+ if(myrank == 0) Earth_mass_total = Earth_mass_total + Earth_mass_total_region
+
+ end subroutine compute_total_Earth_mass
+
diff --git a/src/meshfem3D/compute_element_properties.f90 b/src/meshfem3D/compute_element_properties.f90
index 43d86ff..cad46e3 100644
--- a/src/meshfem3D/compute_element_properties.f90
+++ b/src/meshfem3D/compute_element_properties.f90
@@ -158,7 +158,7 @@
! THREE_D_MODEL_S29EA
! THREE_D_MODEL_GLL
! which show significant transverse isotropy also below 220km depth
- if( USE_VERSION_5_1_5) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT) then
! assignes TI only to elements below (2-layer) ficticious moho down to 670
if( idoubling(ispec)==IFLAG_220_80 &
.or. idoubling(ispec)==IFLAG_80_MOHO &
diff --git a/src/meshfem3D/compute_volumes.f90 b/src/meshfem3D/compute_volumes.f90
index 0ae22c0..7610df7 100644
--- a/src/meshfem3D/compute_volumes.f90
+++ b/src/meshfem3D/compute_volumes.f90
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine compute_volumes(volume_local,area_local_bottom,area_local_top, &
nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
@@ -109,3 +108,71 @@
end subroutine compute_volumes
+!=====================================================================
+
+ ! compute Earth mass of that part of the slice
+
+ subroutine compute_Earth_mass(Earth_mass_local, &
+ nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore)
+
+ use constants
+
+ implicit none
+
+ double precision :: Earth_mass_local
+
+ integer :: nspec
+ double precision :: wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore
+
+ ! local parameters
+ double precision :: weight
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ integer :: i,j,k,ispec
+
+ ! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
+ double precision, parameter :: non_dimensionalizing_factor = RHOAV*R_EARTH**3
+
+ ! initializes
+ Earth_mass_local = ZERO
+
+ ! calculates volume of all elements in mesh
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+
+ ! compute the jacobian
+ xixl = xixstore(i,j,k,ispec)
+ xiyl = xiystore(i,j,k,ispec)
+ xizl = xizstore(i,j,k,ispec)
+ etaxl = etaxstore(i,j,k,ispec)
+ etayl = etaystore(i,j,k,ispec)
+ etazl = etazstore(i,j,k,ispec)
+ gammaxl = gammaxstore(i,j,k,ispec)
+ gammayl = gammaystore(i,j,k,ispec)
+ gammazl = gammazstore(i,j,k,ispec)
+
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ Earth_mass_local = Earth_mass_local + dble(jacobianl)*rhostore(i,j,k,ispec)*weight
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
+ Earth_mass_local = Earth_mass_local * non_dimensionalizing_factor
+
+ end subroutine compute_Earth_mass
+
diff --git a/src/meshfem3D/create_mass_matrices.f90 b/src/meshfem3D/create_mass_matrices.f90
index b90e13f..a539435 100644
--- a/src/meshfem3D/create_mass_matrices.f90
+++ b/src/meshfem3D/create_mass_matrices.f90
@@ -832,7 +832,7 @@
! note: old version (5.1.5)
! only for models where 3D crustal stretching was used (even without topography?)
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
if( CASE_3D ) then
do_ocean_load = .true.
endif
@@ -873,9 +873,9 @@
! slightly move points to avoid roundoff problem when exactly on the polar axis
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
- if( .not. USE_VERSION_5_1_5 ) then
+ if( .not. USE_OLD_VERSION_5_1_5_FORMAT ) then
! adds small margins
- !! DK DK Jul 2013: added a test to only do this if we are on the axis
+ !! DK DK: added a test to only do this if we are on the axis
if(abs(theta) > 89.99d0) then
theta = theta + 0.0000001d0
phi = phi + 0.0000001d0
@@ -888,7 +888,7 @@
! note: bathymetry is given in geographic lat/lon
! (i.e., latitutde with respect to reference ellipsoid)
! we will need convert the geocentric positions here to geographic ones
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
! always converts
theta = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
else
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index 57d9d22..bc805c2 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -47,7 +47,7 @@
use meshfem3D_par,only: &
ibool,idoubling,xstore,ystore,zstore, &
- IMAIN,volume_total,myrank,LOCAL_PATH, &
+ IMAIN,volume_total,Earth_mass_total,myrank,LOCAL_PATH, &
IREGION_CRUST_MANTLE,IREGION_OUTER_CORE,IREGION_INNER_CORE, &
IFLAG_IN_FICTITIOUS_CUBE, &
NCHUNKS,SAVE_MESH_FILES,ABSORBING_CONDITIONS, &
@@ -56,7 +56,7 @@
NGLOB1D_RADIAL_CORNER, &
NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
ADIOS_ENABLED,ADIOS_FOR_ARRAYS_SOLVER, &
- ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION
+ ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,ROLAND_SYLVAIN
use meshfem3D_models_par,only: &
SAVE_BOUNDARY_MESH,SUPPRESS_CRUSTAL_MESH,REGIONAL_MOHO_MESH, &
@@ -104,6 +104,7 @@
! check area and volume of the final mesh
double precision :: area_local_bottom,area_local_top
double precision :: volume_local
+ double precision :: Earth_mass_local
! user output
if(myrank == 0 ) then
@@ -358,6 +359,8 @@
! save the binary files
call synchronize_all()
+!! DK DK for Roland_Sylvain
+ if(.not. ROLAND_SYLVAIN) then
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...saving binary files'
@@ -375,6 +378,7 @@
iregion_code,xstore,ystore,zstore, &
NSPEC2D_TOP,NSPEC2D_BOTTOM)
endif
+ endif
! frees memory
deallocate(rmassx,rmassy,rmassz)
@@ -413,7 +417,15 @@
! computes total area and volume
call compute_area(myrank,NCHUNKS,iregion_code, area_local_bottom, &
- area_local_top, volume_local,volume_total, RCMB,RICB,R_CENTRAL_CUBE)
+ area_local_top,volume_local,volume_total,RCMB,RICB,R_CENTRAL_CUBE)
+
+ ! compute Earth mass of that part of the slice
+ call compute_Earth_mass(Earth_mass_local, &
+ nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore)
+
+ ! computes total Earth mass
+ call compute_total_Earth_mass(myrank,Earth_mass_local,Earth_mass_total)
! create AVS or DX mesh data for the slices
if(SAVE_MESH_FILES) then
@@ -1023,7 +1035,7 @@
! user output
write(IMAIN,'(a,f5.1,a,a,i2.2,a,i2.2,a,i2.2,a)') &
" ",(ilayer_loop-ifirst_region+1.0)/(ilast_region-ifirst_region+1.0) * 100.0,"%", &
- " current clock time is: ",tval(5),"h ",tval(6),"min ",tval(7),"sec"
+ " current clock (NOT elapsed) time is: ",tval(5),"h ",tval(6),"min ",tval(7),"sec"
! flushes I/O buffer
call flush_IMAIN()
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index 745f750..d717443 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine finalize_mesher()
use meshfem3D_par
@@ -56,6 +55,17 @@
endif
endif
+ ! check total Earth mass
+ if(NCHUNKS == 6) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed total Earth mass for this density model and mesh: ',Earth_mass_total
+ write(IMAIN,*) ' (should be not too far from 5.97E+24 kg)'
+ write(IMAIN,*)
+ ! take into account the fact that dimensions have been non-dimensionalized by dividing them by R_EARTH
+ write(IMAIN,*) 'average density for this density model and mesh: ',Earth_mass_total / (volume_total * R_EARTH**3)
+ write(IMAIN,*) ' (should be not too far from 5514 kg/m3)'
+ endif
+
! infos output
numelem_crust_mantle = NSPEC(IREGION_CRUST_MANTLE)
numelem_outer_core = NSPEC(IREGION_OUTER_CORE)
diff --git a/src/meshfem3D/initialize_mesher.f90 b/src/meshfem3D/initialize_mesher.f90
index 04fb743..249fd73 100644
--- a/src/meshfem3D/initialize_mesher.f90
+++ b/src/meshfem3D/initialize_mesher.f90
@@ -75,6 +75,11 @@
! check that the code is running with the requested number of processes
if(sizeprocs /= NPROCTOT) call exit_MPI(myrank,'wrong number of MPI processes')
+!! DK DK for Roland_Sylvain
+ ! in the case of ROLAND_SYLVAIN we should always use double precision
+ if(ROLAND_SYLVAIN .and. CUSTOM_REAL /= SIZE_DOUBLE) &
+ call exit_MPI(myrank,'for ROLAND_SYLVAIN use double precision i.e. configure the code with --enable-double-precision')
+
! synchronizes processes
call synchronize_all()
diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90
index d9afcca..a078adc 100644
--- a/src/meshfem3D/meshfem3D_par.f90
+++ b/src/meshfem3D/meshfem3D_par.f90
@@ -156,6 +156,9 @@
! check area and volume of the final mesh
double precision :: volume_total
+ ! check Earth mass computed in the final mesh
+ double precision :: Earth_mass_total
+
! for loop on all the slices
integer :: iregion_code
integer :: iproc_xi,iproc_eta,ichunk
diff --git a/src/meshfem3D/moho_stretching.f90 b/src/meshfem3D/moho_stretching.f90
index 9fdbca9..e1e6e34 100644
--- a/src/meshfem3D/moho_stretching.f90
+++ b/src/meshfem3D/moho_stretching.f90
@@ -34,7 +34,7 @@
use constants,only: &
NGNOD,R_EARTH_KM,R_EARTH,R_UNIT_SPHERE, &
- PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,USE_VERSION_5_1_5
+ PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,USE_OLD_VERSION_5_1_5_FORMAT
use meshfem3D_par,only: &
RMOHO_FICTITIOUS_IN_MESHER,R220,RMIDDLE_CRUST
@@ -74,7 +74,7 @@
z = zelm(ia)
! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
! note: at this point, the mesh is still perfectly spherical, thus no need to
! convert the geocentric colatitude to a geographic colatitude
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
@@ -105,7 +105,7 @@
! nevertheless its moho depth should be set and will be used in linear stretching
if( moho < TINYVAL ) call exit_mpi(myrank,'error moho depth to honor')
- if( .not. USE_VERSION_5_1_5 ) then
+ if( .not. USE_OLD_VERSION_5_1_5_FORMAT ) then
! limits moho depth to a threshold value to avoid stretching problems
if( moho < MOHO_MINIMUM ) then
print*,'moho value exceeds minimum (in km): ',moho*R_EARTH_KM,MOHO_MINIMUM*R_EARTH_KM,'lat/lon:',lat,lon
@@ -236,7 +236,7 @@
use constants,only: &
NGNOD,R_EARTH_KM,R_EARTH,R_UNIT_SPHERE, &
- PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,HONOR_DEEP_MOHO,USE_VERSION_5_1_5
+ PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,HONOR_DEEP_MOHO,USE_OLD_VERSION_5_1_5_FORMAT
use meshfem3D_par,only: &
R220
@@ -267,7 +267,7 @@
z = zelm(ia)
! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
! note: at this point, the mesh is still perfectly spherical, thus no need to
! convert the geocentric colatitude to a geographic colatitude
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
diff --git a/src/shared/parallel.f90 b/src/shared/parallel.f90
index 179246f..f35d5e9 100644
--- a/src/shared/parallel.f90
+++ b/src/shared/parallel.f90
@@ -1066,6 +1066,30 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine gatherv_all_cr_single_prec(sendbuf, sendcnt, recvbuf, recvcount, recvoffset,recvcounttot, NPROC)
+
+ use constants
+ use mpi
+
+ implicit none
+
+ integer :: sendcnt,recvcounttot,NPROC
+ integer, dimension(NPROC) :: recvcount,recvoffset
+ real, dimension(sendcnt) :: sendbuf
+ real, dimension(recvcounttot) :: recvbuf
+
+ integer :: ier
+
+ call MPI_GATHERV(sendbuf,sendcnt,MPI_REAL, &
+ recvbuf,recvcount,recvoffset,MPI_REAL, &
+ 0,MPI_COMM_WORLD,ier)
+
+ end subroutine gatherv_all_cr_single_prec
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine world_size(sizeval)
diff --git a/src/shared/read_parameter_file.f90 b/src/shared/read_parameter_file.f90
index 684bc93..cc80cdb 100644
--- a/src/shared/read_parameter_file.f90
+++ b/src/shared/read_parameter_file.f90
@@ -259,7 +259,7 @@
if(.not. ROTATION) EXACT_MASS_MATRIX_FOR_ROTATION = .false.
! produces simulations compatible with old globe version 5.1.5
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
print*
print*,'**************'
print*,'using globe version 5.1.5 compatible simulation parameters'
diff --git a/src/shared/rthetaphi_xyz.f90 b/src/shared/rthetaphi_xyz.f90
index 7425064..f56bbbc 100644
--- a/src/shared/rthetaphi_xyz.f90
+++ b/src/shared/rthetaphi_xyz.f90
@@ -213,7 +213,7 @@
! = PI/2 - atan( 1/(1 - e**2) * cos(theta)/sin(theta) )
! = PI/2 - atan( 1.00670466 * cos(theta)/sin(theta) )
- use constants, only: PI_OVER_TWO,TINYVAL,ASSUME_PERFECT_SPHERE,USE_VERSION_5_1_5
+ use constants, only: PI_OVER_TWO,TINYVAL,ASSUME_PERFECT_SPHERE,USE_OLD_VERSION_5_1_5_FORMAT
implicit none
@@ -224,7 +224,7 @@
if(.not. ASSUME_PERFECT_SPHERE) then
! mesh is elliptical
- if( USE_VERSION_5_1_5 ) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT ) then
theta_prime = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
else
! converts geocentric colatitude theta to geographic colatitude theta_prime
@@ -280,7 +280,7 @@
! converts geographic latitude (lat_prime) (in degrees) to geocentric colatitude (theta) (in radians)
- use constants, only: PI_OVER_TWO,DEGREES_TO_RADIANS,ASSUME_PERFECT_SPHERE,USE_VERSION_5_1_5
+ use constants, only: PI_OVER_TWO,DEGREES_TO_RADIANS,ASSUME_PERFECT_SPHERE,USE_OLD_VERSION_5_1_5_FORMAT
implicit none
@@ -289,10 +289,9 @@
! co-latitude (in radians)
double precision,intent(out) :: theta
-
if( .not. ASSUME_PERFECT_SPHERE ) then
! converts geographic (lat_prime) to geocentric latitude and converts to co-latitude (theta)
- if( USE_VERSION_5_1_5) then
+ if( USE_OLD_VERSION_5_1_5_FORMAT) then
! note: factor 0.99329534 = 1 - e**2 with eccentricity e**2 = 0.00670466
theta = PI_OVER_TWO - atan( 0.99329534d0*dtan(lat_prime * DEGREES_TO_RADIANS) )
else
diff --git a/src/specfem3D/iterate_time.F90 b/src/specfem3D/iterate_time.F90
index 907e54c..c908241 100644
--- a/src/specfem3D/iterate_time.F90
+++ b/src/specfem3D/iterate_time.F90
@@ -267,7 +267,7 @@
real :: currenttime
integer :: iglob,inum,data_size
- real(kind=CUSTOM_REAL),dimension(1):: dummy
+ real, dimension(1) :: dummy
! vtk rendering at frame interval
if( mod(it,NTSTEP_BETWEEN_FRAMES) == 0 ) then
@@ -303,14 +303,14 @@
data_size = size(vtkdata)
if( myrank == 0 ) then
! gather data
- call gatherv_all_cr(vtkdata,data_size,&
+ call gatherv_all_cr_single_prec(vtkdata,data_size,&
vtkdata_all,vtkdata_points_all,vtkdata_offset_all, &
vtkdata_numpoints_all,NPROCTOT_VAL)
! updates vtk window
call visualize_vtkdata(it,currenttime,vtkdata_all)
else
! all other process just send data
- call gatherv_all_cr(vtkdata,data_size,&
+ call gatherv_all_cr_single_prec(vtkdata,data_size,&
dummy,vtkdata_points_all,vtkdata_offset_all, &
1,NPROCTOT_VAL)
endif
diff --git a/src/specfem3D/prepare_timerun.f90 b/src/specfem3D/prepare_timerun.f90
index fc4fa4b..8b24e42 100644
--- a/src/specfem3D/prepare_timerun.f90
+++ b/src/specfem3D/prepare_timerun.f90
@@ -2257,8 +2257,8 @@
integer, dimension(:,:),allocatable :: vol_conn_all
integer, dimension(:),allocatable :: vol_conn_offset_all,vol_conn_nspec_all
integer :: vol_nspec_all,ispec_start,ispec_end
- real(kind=CUSTOM_REAL),dimension(1):: dummy
- integer,dimension(1):: dummy_i
+ real,dimension(1) :: dummy
+ integer,dimension(1) :: dummy_i
real(kind=CUSTOM_REAL) :: x,y,z
@@ -2480,13 +2480,13 @@
if( myrank == 0 ) then
! locations
!if( myrank == 0 ) print*," locations..."
- call gatherv_all_cr(free_x,free_np, &
+ call gatherv_all_cr_single_prec(free_x,free_np, &
free_x_all,free_points_all,free_offset_all, &
free_np_all,NPROC)
- call gatherv_all_cr(free_y,free_np, &
+ call gatherv_all_cr_single_prec(free_y,free_np, &
free_y_all,free_points_all,free_offset_all, &
free_np_all,NPROC)
- call gatherv_all_cr(free_z,free_np, &
+ call gatherv_all_cr_single_prec(free_z,free_np, &
free_z_all,free_points_all,free_offset_all, &
free_np_all,NPROC)
@@ -2513,15 +2513,14 @@
free_nspec_all,free_conn_all)
else
- ! all other process just send data
- ! locations
- call gatherv_all_cr(free_x,free_np, &
+ ! all other process just send data locations
+ call gatherv_all_cr_single_prec(free_x,free_np, &
dummy,free_points_all,free_offset_all, &
1,NPROC)
- call gatherv_all_cr(free_y,free_np, &
+ call gatherv_all_cr_single_prec(free_y,free_np, &
dummy,free_points_all,free_offset_all, &
1,NPROC)
- call gatherv_all_cr(free_z,free_np, &
+ call gatherv_all_cr_single_prec(free_z,free_np, &
dummy,free_points_all,free_offset_all, &
1,NPROC)
! connectivity
@@ -2756,13 +2755,13 @@
if( myrank == 0 ) then
! locations
!if( myrank == 0 ) print*," locations..."
- call gatherv_all_cr(vol_x,vol_np, &
+ call gatherv_all_cr_single_prec(vol_x,vol_np, &
vol_x_all,vtkdata_points_all,vtkdata_offset_all, &
vtkdata_numpoints_all,NPROC)
- call gatherv_all_cr(vol_y,vol_np, &
+ call gatherv_all_cr_single_prec(vol_y,vol_np, &
vol_y_all,vtkdata_points_all,vtkdata_offset_all, &
vtkdata_numpoints_all,NPROC)
- call gatherv_all_cr(vol_z,vol_np, &
+ call gatherv_all_cr_single_prec(vol_z,vol_np, &
vol_z_all,vtkdata_points_all,vtkdata_offset_all, &
vtkdata_numpoints_all,NPROC)
@@ -2791,13 +2790,13 @@
else
! all other process just send data
! locations
- call gatherv_all_cr(vol_x,vol_np, &
+ call gatherv_all_cr_single_prec(vol_x,vol_np, &
dummy,vtkdata_points_all,vtkdata_offset_all, &
1,NPROC)
- call gatherv_all_cr(vol_y,vol_np, &
+ call gatherv_all_cr_single_prec(vol_y,vol_np, &
dummy,vtkdata_points_all,vtkdata_offset_all, &
1,NPROC)
- call gatherv_all_cr(vol_z,vol_np, &
+ call gatherv_all_cr_single_prec(vol_z,vol_np, &
dummy,vtkdata_points_all,vtkdata_offset_all, &
1,NPROC)
! connectivity
diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90
index 99be077..6aa9b6c 100644
--- a/src/specfem3D/specfem3D_par.F90
+++ b/src/specfem3D/specfem3D_par.F90
@@ -460,9 +460,9 @@ module specfem_par_crustmantle
integer :: npoints_slice
integer, dimension(NM_KL_REG_PTS_VAL) :: points_slice
integer, dimension(NM_KL_REG_PTS_VAL) :: ispec_reg
- real, dimension(NGLLX, NM_KL_REG_PTS_VAL) :: hxir_reg
- real, dimension(NGLLY, NM_KL_REG_PTS_VAL) :: hetar_reg
- real, dimension(NGLLZ, NM_KL_REG_PTS_VAL) :: hgammar_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLX, NM_KL_REG_PTS_VAL) :: hxir_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLY, NM_KL_REG_PTS_VAL) :: hetar_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLZ, NM_KL_REG_PTS_VAL) :: hgammar_reg
! NOISE_TOMOGRAPHY
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: noise_sourcearray
More information about the CIG-COMMITS
mailing list