[cig-commits] [commit] devel, master: added flag REUSE_EXISTING_OBSERVATION_SURF to be able to reuse an existing observation surface from a previous run. Also added directory utils/Roland_Sylvain_gravity. (a378188)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:17:44 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit a378188b4f7269fe7a05cec7457d0a91a7c67db2
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu May 22 02:18:07 2014 +0200
added flag REUSE_EXISTING_OBSERVATION_SURF to be able to reuse an existing observation surface from a previous run.
Also added directory utils/Roland_Sylvain_gravity.
>---------------------------------------------------------------
a378188b4f7269fe7a05cec7457d0a91a7c67db2
setup/constants.h.in | 5 +
src/meshfem3D/compute_coordinates_grid.f90 | 53 +++++-
src/specfem3D/write_output_ASDF.F90 | 22 ++-
utils/Roland_Sylvain_gravity/GMTglobe.cpt | 53 ++++++
utils/Roland_Sylvain_gravity/Par_file | 203 +++++++++++++++++++++
.../Roland_Sylvain_gravity}/constants.h.in | 39 ++--
.../convert_all_result_files_to_GRD_format.sh | 16 ++
...G_tensor_components_SPECFEM3D_GLOBE_with_GMT.sh | 138 ++++++++++++++
...g_vector_components_SPECFEM3D_GLOBE_with_GMT.sh | 80 ++++++++
...splay_g_vector_norm_SPECFEM3D_GLOBE_with_GMT.sh | 38 ++++
...splay_topography_of_SPECFEM3D_GLOBE_with_GMT.sh | 39 ++++
utils/Roland_Sylvain_gravity/script_MPI_128.sh | 34 ++++
utils/Roland_Sylvain_gravity/submit_all_MPI_128.sh | 19 ++
13 files changed, 712 insertions(+), 27 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 9f06255..1d3e3dc 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -411,6 +411,11 @@
!! DK DK for Roland_Sylvain
logical, parameter :: ROLAND_SYLVAIN = .false.
+! reuse an existing observation surface created in another run and stored to disk,
+! so that we are sure that they are exactly the same (for instance when comparing results for a reference ellipsoidal Earth
+! and results for a 3D Earth with topography)
+ logical, parameter :: REUSE_EXISTING_OBSERVATION_SURF = .false.
+
! altitude of the observation points in km
double precision, parameter :: altitude_of_observation_points = 255.d0
diff --git a/src/meshfem3D/compute_coordinates_grid.f90 b/src/meshfem3D/compute_coordinates_grid.f90
index 2afe74e..badbe03 100644
--- a/src/meshfem3D/compute_coordinates_grid.f90
+++ b/src/meshfem3D/compute_coordinates_grid.f90
@@ -331,15 +331,15 @@
use constants
- use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,lon_observation,lat_observation, &
- ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES
+ use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,x_observation1D,y_observation1D,z_observation1D, &
+ lon_observation,lat_observation,ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES
use meshfem3D_models_par, only: nspl,rspl,espl,espl2,ibathy_topo
implicit none
! local variables
- integer :: ix,iy,ichunk
+ integer :: ix,iy,ichunk,ier
double precision :: gamma,x,y,rgt
double precision :: x_top,y_top,z_top
@@ -350,9 +350,42 @@
! the observation surface is always above the whole Earth, thus each mesh chunk has a size of PI / 2 in each direction
double precision, parameter :: ANGULAR_WIDTH_XI_RAD = PI_OVER_TWO, ANGULAR_WIDTH_ETA_RAD = PI_OVER_TWO
+! reuse an existing observation surface created in another run and stored to disk,
+! so that we are sure that they are exactly the same (for instance when comparing results for a reference ellipsoidal Earth
+! and results for a 3D Earth with topography)
+ if(REUSE_EXISTING_OBSERVATION_SURF) then
+
+ if(myrank == 0) then
+ open(unit=9965,file=trim(OUTPUT_FILES)//'/saved_observation_grid_real_x_y_z_used_by_the_code.txt',status='old', &
+ action='read',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening file for REUSE_EXISTING_OBSERVATION_SURF')
+
+! loop on all the chunks and then on all the observation nodes in each chunk
+ do ichunk = 1,NCHUNKS_MAX
+ do iy = 1,NY_OBSERVATION
+ do ix = 1,NX_OBSERVATION
+ ! read the saved values
+ read(9965,*) x_observation(ix,iy,ichunk),y_observation(ix,iy,ichunk),z_observation(ix,iy,ichunk)
+ enddo
+ enddo
+ enddo
+ close(unit=9965)
+
+ endif
+
+! the 3D and 1D versions of these arrays are equivalenced in the module in which they are declared
+ call bcast_all_dp(x_observation1D, NTOTAL_OBSERVATION)
+ call bcast_all_dp(y_observation1D, NTOTAL_OBSERVATION)
+ call bcast_all_dp(z_observation1D, NTOTAL_OBSERVATION)
+
+ else
+
! for future GMT display
- if(myrank == 0) open(unit=IOUT,file=trim(OUTPUT_FILES)//'/observation_grid_long_lat_topo_for_GMT.txt', &
- status='unknown',action='write')
+ if(myrank == 0) then
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/observation_grid_long_lat_topo_for_GMT.txt',status='unknown',action='write')
+ open(unit=9965,file=trim(OUTPUT_FILES)//'/saved_observation_grid_real_x_y_z_used_by_the_code.txt',status='unknown', &
+ action='write')
+ endif
! loop on all the chunks and then on all the observation nodes in each chunk
do ichunk = 1,NCHUNKS_MAX
@@ -454,12 +487,20 @@
y_observation(ix,iy,ichunk) = y_top * observation_elevation_ratio
z_observation(ix,iy,ichunk) = z_top * observation_elevation_ratio
+ ! store the values obtained so that they can be reused from other runs
+ if(myrank == 0) write(9965,*) x_observation(ix,iy,ichunk),y_observation(ix,iy,ichunk),z_observation(ix,iy,ichunk)
+
enddo
enddo
enddo
! for future GMT display
- if(myrank == 0) close(unit=IOUT)
+ if(myrank == 0) then
+ close(unit=IOUT)
+ close(unit=9965)
+ endif
+
+ endif ! of if(REUSE_EXISTING_OBSERVATION_SURF)
end subroutine compute_observation_surface
diff --git a/src/specfem3D/write_output_ASDF.F90 b/src/specfem3D/write_output_ASDF.F90
index 523ae4c..71c0f28 100644
--- a/src/specfem3D/write_output_ASDF.F90
+++ b/src/specfem3D/write_output_ASDF.F90
@@ -849,8 +849,13 @@ subroutine gather_offset_info(local_dim, global_dim, offset,&
if (ierr /= 0) call exit_MPI (rank, 'Allocate failed.')
call synchronize_all()
- call MPI_Gather(local_dim, 1, MPI_INTEGER, dim_all_proc, 1, &
- MPI_INTEGER, 0, comm, ierr)
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+ call MPI_Gather(local_dim, 1, MPI_INTEGER, dim_all_proc, 1, MPI_INTEGER, 0, comm, ierr)
if(rank==0)then
offset_proc(1)=0
@@ -860,8 +865,17 @@ subroutine gather_offset_info(local_dim, global_dim, offset,&
global_dim=sum(dim_all_proc(1:nproc))
endif
- call MPI_Scatter(offset_proc, 1, MPI_INTEGER, offset, &
- 1, MPI_INTEGER, 0, comm, ierr)
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+ call MPI_Scatter(offset_proc, 1, MPI_INTEGER, offset, 1, MPI_INTEGER, 0, comm, ierr)
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!! DK DK you should use the MPI subroutines in src/shared/parallel.f90 here instead !!!!!!!!!!!!!!!!!
call MPI_Bcast(global_dim, 1, MPI_INTEGER, 0, comm, ierr)
deallocate(dim_all_proc, STAT=ierr)
diff --git a/utils/Roland_Sylvain_gravity/GMTglobe.cpt b/utils/Roland_Sylvain_gravity/GMTglobe.cpt
new file mode 100644
index 0000000..615ac94
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/GMTglobe.cpt
@@ -0,0 +1,53 @@
+# $Id: GMT_globe.cpt,v 1.1 2001/09/23 23:11:20 pwessel Exp $
+#
+# Colormap using in global relief maps
+# Bathymetry colours manually redefined for blue-shade effect and
+# new topography colour scheme for use with Generic Mapping Tools.
+# Designed by Lester M. Anderson (CASP, UK) lester.anderson at casp.cam.ac.uk
+# COLOR_MODEL = RGB
+-10000 153 0 255 -9500 153 0 255
+-9500 153 0 255 -9000 153 0 255
+-9000 153 0 255 -8500 153 0 255
+-8500 136 17 255 -8000 136 17 255
+-8000 119 34 255 -7500 119 34 255
+-7500 102 51 255 -7000 102 51 255
+-7000 85 68 255 -6500 85 68 255
+-6500 68 85 255 -6000 68 85 255
+-6000 51 102 255 -5500 51 102 255
+-5500 34 119 255 -5000 34 119 255
+-5000 17 136 255 -4500 17 136 255
+-4500 0 153 255 -4000 0 153 255
+-4000 27 164 255 -3500 27 164 255
+-3500 54 175 255 -3000 54 175 255
+-3000 81 186 255 -2500 81 186 255
+-2500 108 197 255 -2000 108 197 255
+-2000 134 208 255 -1500 134 208 255
+-1500 161 219 255 -1000 161 219 255
+-1000 188 230 255 -500 188 230 255
+-500 215 241 255 -200 215 241 255
+-200 241 252 255 0 241 252 255
+0 51 102 0 100 51 204 102
+100 51 204 102 200 187 228 146
+200 187 228 146 500 255 220 185
+500 255 220 185 1000 243 202 137
+1000 243 202 137 1500 230 184 88
+1500 230 184 88 2000 217 166 39
+2000 217 166 39 2500 168 154 31
+2500 168 154 31 3000 164 144 25
+3000 164 144 25 3500 162 134 19
+3500 162 134 19 4000 159 123 13
+4000 159 123 13 4500 156 113 7
+4500 156 113 7 5000 153 102 0
+5000 153 102 0 5500 162 89 89
+5500 162 89 89 6000 178 118 118
+6000 178 118 118 6500 183 147 147
+6500 183 147 147 7000 194 176 176
+7000 194 176 176 7500 204 204 204
+7500 204 204 204 8000 229 229 229
+8000 229 229 229 8500 242 242 242
+8500 242 242 242 9000 255 255 255
+9000 255 255 255 9500 255 255 255
+9500 255 255 255 10000 255 255 255
+F 255 255 255
+B 0 0 0
+N 128 128 128
diff --git a/utils/Roland_Sylvain_gravity/Par_file b/utils/Roland_Sylvain_gravity/Par_file
new file mode 100644
index 0000000..2d9b885
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/Par_file
@@ -0,0 +1,203 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # set to 1 for forward simulations, 2 for adjoint simulations for sources, and 3 for kernel simulations
+NOISE_TOMOGRAPHY = 0 # flag of noise tomography, three steps (1,2,3). If earthquake simulation, set it to 0.
+SAVE_FORWARD = .false. # save last frame of forward simulation or not
+
+# number of chunks (1,2,3 or 6)
+NCHUNKS = 6
+
+# angular width of the first chunk (not used if full sphere with six chunks)
+ANGULAR_WIDTH_XI_IN_DEGREES = 20.d0 # angular size of a chunk
+ANGULAR_WIDTH_ETA_IN_DEGREES = 20.d0
+CENTER_LATITUDE_IN_DEGREES = 40.d0
+CENTER_LONGITUDE_IN_DEGREES = 25.d0
+GAMMA_ROTATION_AZIMUTH = 0.d0
+
+# number of elements at the surface along the two sides of the first chunk
+# (must be multiple of 16 and 8 * multiple of NPROC below)
+NEX_XI = 160
+NEX_ETA = 160
+
+# number of MPI processors along the two sides of the first chunk
+NPROC_XI = 10
+NPROC_ETA = 10
+
+# 1D models with real structure:
+# 1D_isotropic_prem, 1D_transversely_isotropic_prem, 1D_iasp91, 1D_1066a, 1D_ak135f_no_mud, 1D_ref, 1D_ref_iso, 1D_jp3d,1D_sea99
+#
+# 1D models with only one fictitious averaged crustal layer:
+# 1D_isotropic_prem_onecrust, 1D_transversely_isotropic_prem_onecrust, 1D_iasp91_onecrust, 1D_1066a_onecrust, 1D_ak135f_no_mud_onecrust
+#
+# fully 3D models:
+# transversely_isotropic_prem_plus_3D_crust_2.0, 3D_anisotropic, 3D_attenuation,
+# s20rts, s40rts, s362ani, s362iso, s362wmani, s362ani_prem, s362ani_3DQ, s362iso_3DQ,
+# s29ea, s29ea,sea99_jp3d1994,sea99,jp3d1994,heterogen
+#
+# 3D models with 1D crust: append "_1Dcrust" the the 3D model name
+# 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 = s40rts # 1D_isotropic_prem
+
+# parameters describing the Earth model
+OCEANS = .false.
+ELLIPTICITY = .true. # .false.
+TOPOGRAPHY = .true. # .false.
+GRAVITY = .false.
+ROTATION = .false.
+ATTENUATION = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .false.
+
+# record length in minutes
+RECORD_LENGTH_IN_MINUTES = 2.5d0
+
+# to undo attenuation for sensitivity kernel calculations or forward runs with SAVE_FORWARD
+# use one (and only one) of the two flags below. UNDO_ATTENUATION is much better (it is exact)
+# but requires a significant amount of disk space for temporary storage.
+ATTENUATION_1D_WITH_3D_STORAGE = .false.
+PARTIAL_PHYS_DISPERSION_ONLY = .false.
+UNDO_ATTENUATION = .false.
+NT_DUMP_ATTENUATION = 100 # how often we dump restart files to undo attenuation, only needed when using UNDO_ATTENUATION
+
+# three mass matrices instead of one are needed to handle rotation very accurately;
+# otherwise rotation is handled slightly less accurately (but still reasonably well);
+# set to .true. if you are interested in precise effects related to rotation;
+# set to .false. if you are solving very large inverse problems at high frequency and also undoing attenuation exactly
+# using the UNDO_ATTENUATION flag above, in which case saving as much memory as possible can be a good idea.
+# You can also safely set it to .false. if you are not in a period range in which rotation matters, e.g. if you are targetting very short-period body waves.
+# if in doubt, set to .true.
+# Set it to .true. if you have ABSORBING_CONDITIONS above, because in that case the code will use the three mass matrices anyway
+# and thus there is no additional cost.
+# this flag is of course unused if ROTATION above is set to .false.
+EXACT_MASS_MATRIX_FOR_ROTATION = .false.
+
+# this for LDDRK high-order time scheme instead of Newmark
+USE_LDDRK = .false.
+
+# the maximum CFL of LDDRK is significantly higher than that of the Newmark scheme,
+# in a ratio that is theoretically 1.327 / 0.697 = 1.15 / 0.604 = 1.903 for a solid with Poisson's ratio = 0.25
+# and for a fluid (see the manual of the 2D code, SPECFEM2D, Tables 4.1 and 4.2, and that ratio does not
+# depend on whether we are in 2D or in 3D). However in practice a ratio of about 1.5 to 1.7 is often safer
+# (for instance for models with a large range of Poisson's ratio values).
+# Since the code computes the time step using the Newmark scheme, for LDDRK we will simply
+# multiply that time step by this ratio when LDDRK is on and when flag INCREASE_CFL_FOR_LDDRK is true.
+INCREASE_CFL_FOR_LDDRK = .true.
+RATIO_BY_WHICH_TO_INCREASE_IT = 1.5d0
+
+# save AVS or OpenDX movies
+#MOVIE_COARSE saves movie only at corners of elements (SURFACE OR VOLUME)
+#MOVIE_COARSE does not work with create_movie_AVS_DX
+MOVIE_SURFACE = .false.
+MOVIE_VOLUME = .false.
+MOVIE_COARSE = .true.
+NTSTEP_BETWEEN_FRAMES = 50
+HDUR_MOVIE = 0.d0
+
+# save movie in volume. Will save element if center of element is in prescribed volume
+# top/bottom: depth in KM, use MOVIE_TOP = -100 to make sure the surface is stored.
+# west/east: longitude, degrees East [-180/180] top/bottom: latitute, degrees North [-90/90]
+# start/stop: frames will be stored at MOVIE_START + i*NSTEP_BETWEEN_FRAMES, where i=(0,1,2..) and iNSTEP_BETWEEN_FRAMES <= MOVIE_STOP
+# movie_volume_type: 1=strain, 2=time integral of strain, 3=\mu*time integral of strain
+# type 4 saves the trace and deviatoric stress in the whole volume, 5=displacement, 6=velocity
+MOVIE_VOLUME_TYPE = 2
+MOVIE_TOP_KM = -100.0
+MOVIE_BOTTOM_KM = 1000.0
+MOVIE_WEST_DEG = -90.0
+MOVIE_EAST_DEG = 90.0
+MOVIE_NORTH_DEG = 90.0
+MOVIE_SOUTH_DEG = -90.0
+MOVIE_START = 0
+MOVIE_STOP = 40000
+
+# save mesh files to check the mesh
+SAVE_MESH_FILES = .false.
+
+# restart files (number of runs can be 1 or higher, choose 1 for no restart files)
+NUMBER_OF_RUNS = 1
+NUMBER_OF_THIS_RUN = 1
+
+# path to store the local database files on each node
+LOCAL_PATH = ./DATABASES_MPI
+# temporary wavefield/kernel/movie files
+LOCAL_TMP_PATH = ./DATABASES_MPI
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 500
+
+# interval in time steps for temporary writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000
+NTSTEP_BETWEEN_READ_ADJSRC = 1000
+
+# output format for the seismograms (one can use either or all of the three formats)
+OUTPUT_SEISMOS_ASCII_TEXT = .true.
+OUTPUT_SEISMOS_SAC_ALPHANUM = .false.
+OUTPUT_SEISMOS_SAC_BINARY = .false.
+OUTPUT_SEISMOS_ASDF = .false.
+
+# rotate seismograms to Radial-Transverse-Z or use default North-East-Z reference frame
+ROTATE_SEISMOGRAMS_RT = .false.
+
+# decide if master process writes all the seismograms or if all processes do it in parallel
+WRITE_SEISMOGRAMS_BY_MASTER = .false.
+
+# save all seismograms in one large combined file instead of one file per seismogram
+# to avoid overloading shared non-local file systems such as GPFS for instance
+SAVE_ALL_SEISMOS_IN_ONE_FILE = .false.
+USE_BINARY_FOR_LARGE_FILE = .false.
+
+# flag to impose receivers at the surface or allow them to be buried
+RECEIVERS_CAN_BE_BURIED = .true.
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
+#-----------------------------------------------------------
+#
+# adjoint kernel outputs
+#
+#-----------------------------------------------------------
+
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL = .false.
+
+# output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels when ANISOTROPIC_KL above is set to .true.
+# means to save radial anisotropic kernels, i.e., sensitivity kernels for beta_v, beta_h, etc.
+SAVE_TRANSVERSE_KL_ONLY = .false.
+
+# output approximate Hessian in crust mantle region.
+# means to save the preconditioning for gradients, they are cross correlations between forward and adjoint accelerations.
+APPROXIMATE_HESS_KL = .false.
+
+# forces transverse isotropy for all mantle elements
+# (default is to use transverse isotropy only between MOHO and 220)
+# means we allow radial anisotropy between the bottom of the crust to the bottom of the transition zone, i.e., 660~km depth.
+USE_FULL_TISO_MANTLE = .false.
+
+# output kernel mask to zero out source region
+# to remove large values near the sources in the sensitivity kernels
+SAVE_SOURCE_MASK = .false.
+
+# output kernels on a regular grid instead of on the GLL mesh points (a bit expensive)
+SAVE_REGULAR_KL = .false.
+
+#-----------------------------------------------------------
+
+# set to true to use GPUs
+GPU_MODE = .false.
+
+# set to true to use the ADIOS library for I/Os
+ADIOS_ENABLED = .false.
+ADIOS_FOR_FORWARD_ARRAYS = .true.
+ADIOS_FOR_MPI_ARRAYS = .true.
+ADIOS_FOR_ARRAYS_SOLVER = .true.
+ADIOS_FOR_SOLVER_MESHFILES = .true.
+ADIOS_FOR_AVS_DX = .true.
+ADIOS_FOR_KERNELS = .true.
+ADIOS_FOR_MODELS = .true.
+
diff --git a/setup/constants.h.in b/utils/Roland_Sylvain_gravity/constants.h.in
similarity index 96%
copy from setup/constants.h.in
copy to utils/Roland_Sylvain_gravity/constants.h.in
index 9f06255..edfd73e 100644
--- a/setup/constants.h.in
+++ b/utils/Roland_Sylvain_gravity/constants.h.in
@@ -63,7 +63,7 @@
integer, parameter :: ISTANDARD_OUTPUT = 6
integer, parameter :: IIN = 40,IOUT = 41
! uncomment this to write messages to a text file
- integer, parameter :: IMAIN = 42
+ integer, parameter :: IMAIN = 6
! uncomment this to write messages to the screen (slows down the code)
! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
! I/O unit for noise data files
@@ -103,12 +103,12 @@
!--- ETOPO4 4-minute model
! size of topography and bathymetry file
- integer, parameter :: NX_BATHY = 5400,NY_BATHY = 2700
+! integer, parameter :: NX_BATHY = 5400,NY_BATHY = 2700
! resolution of topography file in minutes
- integer, parameter :: RESOLUTION_TOPO_FILE = 4
+! integer, parameter :: RESOLUTION_TOPO_FILE = 4
! pathname of the topography file
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_from_etopo2_subsampled.dat'
- character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_smoothed_window_7.dat'
+! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_smoothed_window_7.dat'
!--- ETOPO2 2-minute model
! size of topography and bathymetry file
@@ -125,11 +125,11 @@
!--- ETOPO1 1-minute model
! size of topography and bathymetry file
-! integer, parameter :: NX_BATHY = 21600,NY_BATHY = 10800
+ integer, parameter :: NX_BATHY = 21600,NY_BATHY = 10800
! resolution of topography file in minutes
-! integer, parameter :: RESOLUTION_TOPO_FILE = 1
+ integer, parameter :: RESOLUTION_TOPO_FILE = 1
! pathname of the topography file
-! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo1_ice_c_original_unmodified_unsmoothed.dat'
+ character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo1_ice_c_original_unmodified_unsmoothed.dat'
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo1_ice_c_smoothed_window_3.dat'
@@ -167,7 +167,7 @@
integer, parameter :: TOPO_MAXIMUM = + 9000 ! (height in m, Mount Everest)
! minimum thickness in meters to include the effect of the oceans and topo
- double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 50.d0
+ double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 2.d0 !!!! 100.d0
!!-----------------------------------------------------------
@@ -196,15 +196,15 @@
! use sedimentary layers in crustal model
logical, parameter :: INCLUDE_SEDIMENTS_IN_CRUST = .true.
- logical, parameter :: INCLUDE_ICE_IN_CRUST = .false. ! always set this to false except for Roland_Sylvain gravity calculations
- double precision, parameter :: MINIMUM_SEDIMENT_THICKNESS = 2.d0 ! minimim thickness in km
+ logical, parameter :: INCLUDE_ICE_IN_CRUST = .true.
+ double precision, parameter :: MINIMUM_SEDIMENT_THICKNESS = -100000.d0 ! minimim thickness in km
! default crustal model
! (used as default when CRUSTAL flag is set for simulation)
!-- uncomment for using Crust1.0
-! integer, parameter :: ITYPE_CRUSTAL_MODEL = ICRUST_CRUST1
+ integer, parameter :: ITYPE_CRUSTAL_MODEL = ICRUST_CRUST1
!-- uncomment for using Crust2.0
- integer, parameter :: ITYPE_CRUSTAL_MODEL = ICRUST_CRUST2
+! integer, parameter :: ITYPE_CRUSTAL_MODEL = ICRUST_CRUST2
!-- uncomment for using General Crustmaps instead
! integer, parameter :: ITYPE_CRUSTAL_MODEL = ICRUST_CRUSTMAPS
!-- uncomment for using EPcrust instead (European regional model)
@@ -409,7 +409,12 @@
!!-----------------------------------------------------------
!! DK DK for Roland_Sylvain
- logical, parameter :: ROLAND_SYLVAIN = .false.
+ logical, parameter :: ROLAND_SYLVAIN = .true.
+
+! reuse an existing observation surface created in another run and stored to disk,
+! so that we are sure that they are exactly the same (for instance when comparing results for a reference ellipsoidal Earth
+! and results for a 3D Earth with topography)
+ logical, parameter :: REUSE_EXISTING_OBSERVATION_SURF = .true.
! altitude of the observation points in km
double precision, parameter :: altitude_of_observation_points = 255.d0
@@ -427,13 +432,13 @@
! number of points in each horizontal direction of the observation grid of each cubed-sphere chunk
! at the altitude of the observation point
-!! DK DK 4 is a fictitious value used to save memory when the ROLAND_SYLVAIN option is off
- integer, parameter :: NX_OBSERVATION = 4 ! 500
+!! DK DK 2 is a fictitious value used to save memory when the ROLAND_SYLVAIN option is off
+ integer, parameter :: NX_OBSERVATION = 500 ! 10 !!!!!!!!!!!!! 100
integer, parameter :: NY_OBSERVATION = NX_OBSERVATION
! the code will display sample output values at this particular point as a check
- integer, parameter :: ixr = max(NX_OBSERVATION / 3, 1)
- integer, parameter :: iyr = max(NY_OBSERVATION / 4, 1)
+ integer, parameter :: ixr = NX_OBSERVATION / 3
+ integer, parameter :: iyr = NY_OBSERVATION / 4
integer, parameter :: ichunkr = 3
! how often (every how many spectral elements computed) we print a timestamp to monitor the behavior of the code
diff --git a/utils/Roland_Sylvain_gravity/convert_all_result_files_to_GRD_format.sh b/utils/Roland_Sylvain_gravity/convert_all_result_files_to_GRD_format.sh
new file mode 100644
index 0000000..95a3630
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/convert_all_result_files_to_GRD_format.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+
+ surface results_g_x_for_GMT.txt -Gresults_g_x.grd -Rd -I4m -f0x,1y -V
+ surface results_g_y_for_GMT.txt -Gresults_g_y.grd -Rd -I4m -f0x,1y -V
+ surface results_g_z_for_GMT.txt -Gresults_g_z.grd -Rd -I4m -f0x,1y -V
+ surface results_norm_of_g_for_GMT.txt -Gresults_norm_of_g.grd -Rd -I4m -f0x,1y -V
+ surface results_G_xx_for_GMT.txt -Gresults_G_xx.grd -Rd -I4m -f0x,1y -V
+ surface results_G_yy_for_GMT.txt -Gresults_G_yy.grd -Rd -I4m -f0x,1y -V
+ surface results_G_zz_for_GMT.txt -Gresults_G_zz.grd -Rd -I4m -f0x,1y -V
+ surface results_G_xy_for_GMT.txt -Gresults_G_xy.grd -Rd -I4m -f0x,1y -V
+ surface results_G_xz_for_GMT.txt -Gresults_G_xz.grd -Rd -I4m -f0x,1y -V
+ surface results_G_yz_for_GMT.txt -Gresults_G_yz.grd -Rd -I4m -f0x,1y -V
+
+# Clean up
+rm -f .gmt*
+
diff --git a/utils/Roland_Sylvain_gravity/display_G_tensor_components_SPECFEM3D_GLOBE_with_GMT.sh b/utils/Roland_Sylvain_gravity/display_G_tensor_components_SPECFEM3D_GLOBE_with_GMT.sh
new file mode 100644
index 0000000..78bbaab
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/display_G_tensor_components_SPECFEM3D_GLOBE_with_GMT.sh
@@ -0,0 +1,138 @@
+#!/bin/bash
+
+# The first plotting command needs - and _only_ needs "-K"
+# All following but the last need both, "-K" and "-O"
+# The last one needs _only_ "-O"
+
+central_meridian=0
+
+gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+# DK DK for the difference between the 3D calculation and the 1D reference calculation
+makecpt -T-10/10/0.01 -Z > color2.cpt
+
+################################
+
+file=results_G_xx
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_xx_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_G_yy
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_yy_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_G_zz
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_zz_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_G_xy
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_xy_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_G_xz
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_xz_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_G_yz
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."G_yz_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B2.0:"Difference with reference ellipsoid (in Eotvos)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+# Clean up
+rm -f .gmt* gmt.conf gmt.history color2.cpt
+
diff --git a/utils/Roland_Sylvain_gravity/display_g_vector_components_SPECFEM3D_GLOBE_with_GMT.sh b/utils/Roland_Sylvain_gravity/display_g_vector_components_SPECFEM3D_GLOBE_with_GMT.sh
new file mode 100644
index 0000000..bdf4e3c
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/display_g_vector_components_SPECFEM3D_GLOBE_with_GMT.sh
@@ -0,0 +1,80 @@
+#!/bin/bash
+
+# The first plotting command needs - and _only_ needs "-K"
+# All following but the last need both, "-K" and "-O"
+# The last one needs _only_ "-O"
+
+central_meridian=0
+
+gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+# DK DK for absolute norm of the g vector
+#makecpt -Chaxby -T-10/10/0.001 > color2.cpt
+# DK DK for the difference between the 3D calculation and the 1D reference calculation
+makecpt -T-0.03/0.03/0.0001 -Z > color2.cpt
+
+################################
+
+file=results_g_x
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."g_x_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B0.006:"Difference with reference ellipsoid (in m/s2)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_g_y
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."g_y_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B0.006:"Difference with reference ellipsoid (in m/s2)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+file=results_g_z
+
+## DK DK subtract the 1D result for the reference ellipsoid calculation from the 3D result with topography and s20rts
+grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."g_z_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B0.006:"Difference with reference ellipsoid (in m/s2)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+################################
+
+# Clean up
+rm -f .gmt* gmt.conf gmt.history color2.cpt
+
diff --git a/utils/Roland_Sylvain_gravity/display_g_vector_norm_SPECFEM3D_GLOBE_with_GMT.sh b/utils/Roland_Sylvain_gravity/display_g_vector_norm_SPECFEM3D_GLOBE_with_GMT.sh
new file mode 100644
index 0000000..7592438
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/display_g_vector_norm_SPECFEM3D_GLOBE_with_GMT.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+# The first plotting command needs - and _only_ needs "-K"
+# All following but the last need both, "-K" and "-O"
+# The last one needs _only_ "-O"
+
+central_meridian=0
+
+gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+# DK DK for absolute norm of the g vector
+#makecpt -T9.0/9.2/0.001 -Z > color2.cpt
+# DK DK for the difference between the 3D calculation and the 1D reference calculation
+makecpt -T-0.03/0.03/0.0001 -Z > color2.cpt
+
+################################
+
+file=results_norm_of_g
+
+## DK DK subtract the 3D result with topography and s20rts from the 1D result for the reference ellipsoid calculation
+#grdmath new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_topo_and_s40rts/${file}.grd new_SPECFEM3D_GLOBE_NEX_160_on_600_cores_whole_Earth_with_ellipticity_only_and_PREM1D/${file}.grd NEG ADD = ${file}_difference.grd
+
+## DK DK offset de 4.5cm
+## DK DK -Rd signifie -R-180/180/-90/90
+grdimage ${file}_difference.grd -Ccolor2.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > ${file}_difference.ps
+
+## DK DK offset de -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."norm_of_g_difference": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE gravity calculations by Dimitri Komatitsch" -V -O -K >> ${file}_difference.ps
+
+psscale -Ccolor2.cpt -D12.5/-1.5/16/0.25h -B0.006:"Difference with reference ellipsoid (in m/s2)": -V -O >> ${file}_difference.ps
+
+## DK DK convert the final file to PDF
+ps2pdf ${file}_difference.ps
+rm -f ${file}_difference.ps
+
+# Clean up
+rm -f .gmt* gmt.conf gmt.history color2.cpt
+
diff --git a/utils/Roland_Sylvain_gravity/display_topography_of_SPECFEM3D_GLOBE_with_GMT.sh b/utils/Roland_Sylvain_gravity/display_topography_of_SPECFEM3D_GLOBE_with_GMT.sh
new file mode 100644
index 0000000..4eb3154
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/display_topography_of_SPECFEM3D_GLOBE_with_GMT.sh
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+### Dimitri Komatitsch, CNRS Marseille, France, May 2014
+
+# The first plotting command needs and _only_ needs "-K"
+# All following but the last need both "-K" and "-O"
+# The last one needs _only_ "-O"
+
+if [[ ! -f GMTglobe.cpt ]]; then
+ echo "Color palette file GMTglobe.cpt not found!"
+ exit 1
+fi
+
+ps=topography_of_SPECFEM_mesh.ps
+
+central_meridian=0
+
+rm -f $ps
+
+gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+# convert the non-regular cubed-sphere grid output of SPECFEM to a regular grid for GMT
+ surface observation_grid_long_lat_topo_for_GMT.txt -Gtopography_of_SPECFEM_mesh.grd -Rd -I4m -f0x,1y -V
+
+## DK DK offset of 4.5cm
+## DK DK -Rd is an alias for -R-180/180/-90/90
+grdimage topography_of_SPECFEM_mesh.grd -CGMTglobe.cpt -Rd -JK$central_meridian/9i -Y4.5c -K -V > $ps
+
+## DK DK offset of -1.5 inch
+pscoast -Rd -JK$central_meridian/9i -B45g30:."Topography of SPECFEM3D_GLOBE mesh": -W -Dc -A1000 -U/-0.75i/-1.5i/"SPECFEM3D_GLOBE calculations by Dimitri Komatitsch" -V -K -O >> $ps
+
+psscale -CGMTglobe.cpt -D12.5/-1.5/16/0.25h -B2000.:"Elevation (m)": -V -O >> $ps
+
+## DK DK convert the final file to PDF
+ps2pdf $ps
+
+# Clean up
+rm -f .gmt* gmt.conf gmt.history topography_of_SPECFEM_mesh.grd $ps
+
diff --git a/utils/Roland_Sylvain_gravity/script_MPI_128.sh b/utils/Roland_Sylvain_gravity/script_MPI_128.sh
new file mode 100755
index 0000000..8be03dc
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/script_MPI_128.sh
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+########### -n : number of tasks that will be used in parallel mode (default=1)
+########### -N : number of nodes to allocate for parallel usage (default is chosen by the underlying system)
+########### -c : number of cores per parallel task to allocate (default=1)
+#MSUB -n 600 # Reservation des coeurs de processeurs
+#MSUB -N 38 # nombre total de noeuds a utiliser (mettre n/16 car il y a 16 coeurs par noeud)
+#MSUB -c 1 # 1 coeur par tache MPI
+
+#################### #MSUB -E '-R select[rank!=1004&rank!=1005]'
+
+##################MSUB -A gen7165
+#MSUB -A gen6351
+#MSUB -r SPECFEM128 # Nom du job
+#MSUB -T 86400 # Limite de temps elapsed du job en secondes
+#BSUB -x # mode exclusif
+#MSUB -o output_night_%I.o # Sortie standard et %I est le job_ID
+#MSUB -e output_night_%I.e # Sortie d'erreur et %I est le job_ID
+##############MSUB -j oe # join output and error
+#MSUB -k n # do not keep older output and error from previous jobs with the same name
+
+set -x
+
+# ${BRIDGE_MSUB_PWD} est une variable d'environnement representant le repertoire de soumission
+cd ${BRIDGE_MSUB_PWD}
+
+#export CUDA_PROFILE=1
+#export CUDA_PROFILE_LOG=dimitri_CUDA_profile.log
+#export CUDA_PROFILE_CONFIG=dimitri_CUDA_profile_config.txt
+
+##########################
+
+ccc_mprun ./bin/xmeshfem3D
+
diff --git a/utils/Roland_Sylvain_gravity/submit_all_MPI_128.sh b/utils/Roland_Sylvain_gravity/submit_all_MPI_128.sh
new file mode 100755
index 0000000..e5d1676
--- /dev/null
+++ b/utils/Roland_Sylvain_gravity/submit_all_MPI_128.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+/bin/rm -f output_night*.o output_night*.e core.* seismo*txt
+
+/bin/rm -f -r DATABASES_MPI
+/bin/rm -f -r OUTPUT_FILES
+
+mkdir -p DATABASES_MPI OUTPUT_FILES
+mkdir -p OUTPUT_FILES/DATABASES_MPI
+
+if [[ -f saved_observation_grid_real_x_y_z_used_by_the_code.txt ]]; then
+ cp saved_observation_grid_real_x_y_z_used_by_the_code.txt OUTPUT_FILES
+fi
+
+ulimit -S -s unlimited
+
+#ccc_msub -p genb002 -q gpu script_MPI_030.sh
+ccc_msub -q standard script_MPI_128.sh
+
More information about the CIG-COMMITS
mailing list