[cig-commits] [commit] devel: added Roland_Sylvain gravity integrals; also simplified run_this_example.sh (d2e6f20)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 11 19:09:37 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/d29a156f26a222b10b5a7129c39318d885ddf3fb...2b6a4170c869377737a04ca3c5561ac3f345d7f9
>---------------------------------------------------------------
commit d2e6f20591bd0144125edfd94a6e6b04185af32b
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu Jun 12 04:07:47 2014 +0200
added Roland_Sylvain gravity integrals; also simplified run_this_example.sh
>---------------------------------------------------------------
d2e6f20591bd0144125edfd94a6e6b04185af32b
examples | 2 +-
run_this_example.sh | 60 ++++++-
setup/constants.h.in | 32 ++++
src/shared/parallel.f90 | 18 +++
src/specfem3D/prepare_timerun.F90 | 331 ++++++++++++++++++++++++++++++++++++++
src/specfem3D/specfem3D_par.f90 | 4 +
6 files changed, 445 insertions(+), 2 deletions(-)
diff --git a/examples b/examples
index e1bcf45..7fff45c 160000
--- a/examples
+++ b/examples
@@ -1 +1 @@
-Subproject commit e1bcf45fba976398509239f8ca07e29f9afd9719
+Subproject commit 7fff45cfbe947337e9e3f5cac572e4f5cea78a38
diff --git a/run_this_example.sh b/run_this_example.sh
deleted file mode 120000
index e0583f3..0000000
--- a/run_this_example.sh
+++ /dev/null
@@ -1 +0,0 @@
-examples/homogeneous_halfspace/run_this_example.sh
\ No newline at end of file
diff --git a/run_this_example.sh b/run_this_example.sh
new file mode 100755
index 0000000..f3d805b
--- /dev/null
+++ b/run_this_example.sh
@@ -0,0 +1,59 @@
+#!/bin/bash
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 15 minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+rm -f -r OUTPUT_FILES
+
+mkdir -p bin
+mkdir -p OUTPUT_FILES
+mkdir -p OUTPUT_FILES/DATABASES_MPI
+
+# note: ./configure will default to gfortran configuration, unless you uncomment the line below
+# ./configure FC=ifort CC=icc MPIFC=mpif90
+./configure
+make clean
+make all
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/CMTSOLUTION OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+
+# get the number of processors
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+# decomposes mesh using the pre-saved mesh files in MESH-default
+echo
+echo " decomposing mesh..."
+echo
+./bin/xdecompose_mesh $NPROC ./MESH-default ./OUTPUT_FILES/DATABASES_MPI/
+
+# runs database generation
+echo
+echo " running database generation on $NPROC processors..."
+echo
+mpirun -np $NPROC ./bin/xgenerate_databases
+
+# runs simulation
+echo
+echo " running solver on $NPROC processors..."
+echo
+mpirun -np $NPROC ./bin/xspecfem3D
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
+
+
diff --git a/setup/constants.h.in b/setup/constants.h.in
index a0bf03e..a98f322 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -125,6 +125,38 @@
!!-----------------------------------------------------------
!!
+!! Roland_Sylvain gravity calculations
+!!
+!!-----------------------------------------------------------
+
+!! DK DK for Roland_Sylvain
+ logical, parameter :: ROLAND_SYLVAIN = .false. !!! .true.
+
+! check for negative Jacobians in the calculation of integrals or not
+! (can safely be done once to check that the mesh is OK at a given resolution, and then permanently
+! turned off in future runs because the mesh does not change)
+ logical, parameter :: CHECK_FOR_NEGATIVE_JACOBIANS = .true.
+
+! file that contains the observation grid (which is a simple ASCII file, with x y z of each point on separate lines)
+ character(len=*), parameter :: OBSERVATION_GRID_FILE = './DATA/observation_grid_to_use_for_gravity.txt'
+
+! number of points in the observation grid
+ integer, parameter :: NTOTAL_OBSERVATION = 4
+
+! a receiver at which we check and output the gravity field (for display purposes only; you can thus safely leave it unchanged)
+ integer, parameter :: iobs_receiver = 1
+
+! how often (every how many spectral elements computed) we print a timestamp to monitor the behavior of the code
+ integer, parameter :: NSPEC_DISPLAY_INTERVAL = 1000
+
+! definition of an Eotvos compared to S.I. units.
+! The unit of gravity gradient is the Eotvos, which is equivalent to 1e-9 s-2 (or 1e-4 mGal/m).
+! A person walking at a distance of 2 meters provides a gravity gradient signal of approximately one Eotvos.
+! Mountains can create signals of several hundred Eotvos.
+ double precision, parameter :: SI_UNITS_TO_EOTVOS = 1.d+9
+
+!!-----------------------------------------------------------
+!!
!! source/receiver setup
!!
!!-----------------------------------------------------------
diff --git a/src/shared/parallel.f90 b/src/shared/parallel.f90
index e4d1191..577f795 100644
--- a/src/shared/parallel.f90
+++ b/src/shared/parallel.f90
@@ -668,6 +668,24 @@
!----
!
+ subroutine sum_all_1Darray_dp(sendbuf, recvbuf, nx)
+
+ use mpi
+
+ implicit none
+
+ integer :: nx
+ double precision, dimension(nx) :: sendbuf, recvbuf
+ integer :: ier
+
+ call MPI_REDUCE(sendbuf,recvbuf,nx,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ end subroutine sum_all_1Darray_dp
+
+!
+!----
+!
+
subroutine sum_all_all_cr(sendbuf, recvbuf)
use mpi
diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90
index 1b6ac40..eb8fa03 100644
--- a/src/specfem3D/prepare_timerun.F90
+++ b/src/specfem3D/prepare_timerun.F90
@@ -100,6 +100,16 @@
call prepare_timerun_OpenMP()
#endif
+ ! compute the Roland_Sylvain gravity integrals if needed
+ if(ROLAND_SYLVAIN) then
+ if( myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) ' ...computing Roland_Sylvain gravity integrals'
+ call flush_IMAIN()
+ endif
+ call compute_Roland_Sylvain_integr()
+ endif
+
! elapsed time since beginning of preparation
if (myrank == 0) then
tCPU = wtime() - time_start
@@ -1442,3 +1452,324 @@
end subroutine prepare_timerun_OpenMP
#endif
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ ! compute Roland_Sylvain integrals of that slice, and then total integrals for the whole mesh
+
+ subroutine compute_Roland_Sylvain_integr()
+
+ use constants
+
+ use specfem_par
+
+ implicit none
+
+ ! local parameters
+ double precision :: weight
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl
+ double precision :: jacobianl
+ integer :: i,j,k,ispec,iglob,ier
+ double precision :: xval,yval,zval
+ double precision :: xval_squared,yval_squared,zval_squared
+ double precision :: x_meshpoint,y_meshpoint,z_meshpoint
+ double precision :: distance_squared,distance_cubed, &
+ three_over_distance_squared,one_over_distance_cubed,three_over_distance_fifth_power
+ double precision :: common_multiplying_factor,common_mult_times_one_over,common_mult_times_three_over
+
+ integer :: iobservation
+
+! read the observation surface
+ x_observation(:) = 0.d0
+ y_observation(:) = 0.d0
+ z_observation(:) = 0.d0
+ if(myrank == 0) then
+ open(unit=IIN,file=OBSERVATION_GRID_FILE(1:len_trim(OBSERVATION_GRID_FILE)),status='old',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening file observation_grid_to_use_for_gravity.txt')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ read(IIN,*) x_observation(iobservation),y_observation(iobservation),z_observation(iobservation)
+ enddo
+ close(unit=IIN)
+ endif
+
+! broadcast the observation surface read
+ call bcast_all_dp(x_observation, NTOTAL_OBSERVATION)
+ call bcast_all_dp(y_observation, NTOTAL_OBSERVATION)
+ call bcast_all_dp(z_observation, NTOTAL_OBSERVATION)
+
+! initialize the gravity arrays
+ g_x(:) = 0.d0
+ g_y(:) = 0.d0
+ g_z(:) = 0.d0
+
+ G_xx(:) = 0.d0
+ G_yy(:) = 0.d0
+ G_zz(:) = 0.d0
+ G_xy(:) = 0.d0
+ G_xz(:) = 0.d0
+ G_yz(:) = 0.d0
+
+ ! calculates volume of all elements in mesh
+ do ispec = 1,NSPEC_AB
+
+ ! print information about number of elements done so far
+ if(myrank == 0 .and. (mod(ispec,NSPEC_DISPLAY_INTERVAL) == 0 .or. ispec == 1 .or. ispec == NSPEC_AB)) then
+ write(IMAIN,*) 'for Roland_Sylvain integrals, ',ispec,' elements computed out of ',NSPEC_AB
+ ! write time stamp file to give information about progression of the calculation of gravity integrals
+ write(outputname,"('/timestamp_gravity_calculations_ispec',i7.7,'_out_of_',i7.7)") ispec,NSPEC_AB
+ ! timestamp file output
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//outputname,status='unknown',action='write')
+ write(IOUT,*) ispec,' elements done for gravity calculations out of ',NSPEC_AB
+ close(unit=IOUT)
+ endif
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+
+ ! compute the Jacobian
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! do this in double precision for accuracy
+ jacobianl = 1.d0 / dble(xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ if(CHECK_FOR_NEGATIVE_JACOBIANS .and. jacobianl <= ZERO) stop 'error: negative Jacobian found in integral calculation'
+
+ iglob = ibool(i,j,k,ispec)
+ x_meshpoint = xstore(iglob)
+ y_meshpoint = ystore(iglob)
+ z_meshpoint = zstore(iglob)
+
+ common_multiplying_factor = jacobianl * weight * rhostore(i,j,k,ispec)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!! beginning of loop on all the data to create
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! loop on all the points in the observation surface
+ do iobservation = 1,NTOTAL_OBSERVATION
+
+ xval = x_meshpoint - x_observation(iobservation)
+ yval = y_meshpoint - y_observation(iobservation)
+ zval = z_meshpoint - z_observation(iobservation)
+
+ xval_squared = xval**2
+ yval_squared = yval**2
+ zval_squared = zval**2
+
+ distance_squared = xval_squared + yval_squared + zval_squared
+ distance_cubed = distance_squared * sqrt(distance_squared)
+
+ three_over_distance_squared = 3.d0 / distance_squared
+ one_over_distance_cubed = 1.d0 / distance_cubed
+ three_over_distance_fifth_power = three_over_distance_squared * one_over_distance_cubed
+
+ common_mult_times_one_over = common_multiplying_factor * one_over_distance_cubed
+ common_mult_times_three_over = common_multiplying_factor * three_over_distance_fifth_power
+
+ g_x(iobservation) = g_x(iobservation) + common_mult_times_one_over * xval
+ g_y(iobservation) = g_y(iobservation) + common_mult_times_one_over * yval
+ g_z(iobservation) = g_z(iobservation) + common_mult_times_one_over * zval
+
+ G_xx(iobservation) = G_xx(iobservation) + common_mult_times_one_over * (xval_squared * three_over_distance_squared - 1.d0)
+ G_yy(iobservation) = G_yy(iobservation) + common_mult_times_one_over * (yval_squared * three_over_distance_squared - 1.d0)
+ G_zz(iobservation) = G_zz(iobservation) + common_mult_times_one_over * (zval_squared * three_over_distance_squared - 1.d0)
+
+ G_xy(iobservation) = G_xy(iobservation) + common_mult_times_three_over * xval*yval
+ G_xz(iobservation) = G_xz(iobservation) + common_mult_times_three_over * xval*zval
+ G_yz(iobservation) = G_yz(iobservation) + common_mult_times_three_over * yval*zval
+
+ enddo
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!! end of loop on all the data to create
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! the result is displayed in Eotvos = 1.e+9 s-2
+ G_xx(:) = G_xx(:) * SI_UNITS_TO_EOTVOS
+ G_yy(:) = G_yy(:) * SI_UNITS_TO_EOTVOS
+ G_zz(:) = G_zz(:) * SI_UNITS_TO_EOTVOS
+ G_xy(:) = G_xy(:) * SI_UNITS_TO_EOTVOS
+ G_xz(:) = G_xz(:) * SI_UNITS_TO_EOTVOS
+ G_yz(:) = G_yz(:) * SI_UNITS_TO_EOTVOS
+
+ ! use an MPI reduction to compute the total value of the integral into a temporary array
+ ! and then copy it back into the original array
+ call sum_all_1Darray_dp(g_x,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) g_x(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(g_y,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) g_y(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(g_z,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) g_z(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_xx,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_xx(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_yy,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_yy(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_zz,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_zz(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_xy,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_xy(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_xz,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_xz(:) = temporary_array_for_sum(:)
+
+ call sum_all_1Darray_dp(G_yz,temporary_array_for_sum,NTOTAL_OBSERVATION)
+ if(myrank == 0) G_yz(:) = temporary_array_for_sum(:)
+
+ !--- print number of points and elements in the mesh for each region
+ if(myrank == 0) then
+
+ temporary_array_for_sum(:) = sqrt(g_x(:)**2 + g_y(:)**2 + g_z(:)**2)
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of norm of g vector on whole observation surface = ',minval(temporary_array_for_sum),' m.s-2'
+ write(IMAIN,*) 'maxval of norm of g vector on whole observation surface = ',maxval(temporary_array_for_sum),' m.s-2'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_xx on whole observation surface = ',minval(G_xx),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xx on whole observation surface = ',maxval(G_xx),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_yy on whole observation surface = ',minval(G_yy),' Eotvos'
+ write(IMAIN,*) 'maxval of G_yy on whole observation surface = ',maxval(G_yy),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_zz on whole observation surface = ',minval(G_zz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_zz on whole observation surface = ',maxval(G_zz),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_xy on whole observation surface = ',minval(G_xy),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xy on whole observation surface = ',maxval(G_xy),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_xz on whole observation surface = ',minval(G_xz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xz on whole observation surface = ',maxval(G_xz),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'minval of G_yz on whole observation surface = ',minval(G_yz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_yz on whole observation surface = ',maxval(G_yz),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'Minval and maxval of trace of G, which in principle should be zero:'
+ write(IMAIN,*)
+ temporary_array_for_sum(:) = abs(G_xx(:) + G_yy(:) + G_zz(:))
+ write(IMAIN,*) 'minval of abs(G_xx + G_yy + G_zz) on whole observation surface = ',minval(temporary_array_for_sum),' Eotvos'
+ write(IMAIN,*) 'maxval of abs(G_xx + G_yy + G_zz) on whole observation surface = ',maxval(temporary_array_for_sum),' Eotvos'
+
+ write(IMAIN,*)
+ write(IMAIN,*) '-----------------------------'
+ write(IMAIN,*)
+ write(IMAIN,*) 'displaying the fields computed at observation point = ',iobs_receiver,' out of ',NTOTAL_OBSERVATION
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed g_x = ',g_x(iobs_receiver),' m.s-2'
+ write(IMAIN,*) 'computed g_y = ',g_y(iobs_receiver),' m.s-2'
+ write(IMAIN,*) 'computed g_z = ',g_z(iobs_receiver),' m.s-2'
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed norm of g vector = ',sqrt(g_x(iobs_receiver)**2 + g_y(iobs_receiver)**2 + &
+ g_z(iobs_receiver)**2),' m.s-2'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed G_xx = ',G_xx(iobs_receiver),' Eotvos'
+ write(IMAIN,*) 'computed G_yy = ',G_yy(iobs_receiver),' Eotvos'
+ write(IMAIN,*) 'computed G_zz = ',G_zz(iobs_receiver),' Eotvos'
+ write(IMAIN,*)
+ write(IMAIN,*) 'G tensor should be traceless, G_xx + G_yy + G_zz = 0.'
+ write(IMAIN,*) 'Actual sum obtained = ',G_xx(iobs_receiver) + G_yy(iobs_receiver) + G_zz(iobs_receiver)
+ if(max(abs(G_xx(iobs_receiver)),abs(G_yy(iobs_receiver)),abs(G_zz(iobs_receiver))) > TINYVAL) &
+ write(IMAIN,*) ' i.e., ',sngl(100.d0*abs(G_xx(iobs_receiver) + G_yy(iobs_receiver) + G_zz(iobs_receiver)) / &
+ max(abs(G_xx(iobs_receiver)),abs(G_yy(iobs_receiver)),abs(G_zz(iobs_receiver)))), &
+ '% of max(abs(G_xx),abs(G_yy),abs(G_zz))'
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed G_xy = ',G_xy(iobs_receiver),' Eotvos'
+ write(IMAIN,*) 'computed G_xz = ',G_xz(iobs_receiver),' Eotvos'
+ write(IMAIN,*) 'computed G_yz = ',G_yz(iobs_receiver),' Eotvos'
+
+ ! save the results
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_g_x_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) g_x(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_g_y_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) g_y(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_g_z_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) g_z(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_norm_of_g_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) sqrt(g_x(iobservation)**2 + g_y(iobservation)**2 + g_z(iobservation)**2)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_xx_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_xx(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_yy_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_yy(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_zz_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_zz(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_xy_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_xy(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_xz_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_xz(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES_PATH)//'/results_G_yz_for_GMT.txt',status='unknown',action='write')
+ do iobservation = 1,NTOTAL_OBSERVATION
+ write(IOUT,*) G_yz(iobservation)
+ enddo
+ close(unit=IOUT)
+
+ endif
+
+ end subroutine compute_Roland_Sylvain_integr
+
diff --git a/src/specfem3D/specfem3D_par.f90 b/src/specfem3D/specfem3D_par.f90
index d4eb3ab..daf0b83 100644
--- a/src/specfem3D/specfem3D_par.f90
+++ b/src/specfem3D/specfem3D_par.f90
@@ -280,6 +280,10 @@ module specfem_par
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
normal_x_noise,normal_y_noise,normal_z_noise, mask_noise
+ ! for Roland_Sylvain integrals
+ double precision, dimension(NTOTAL_OBSERVATION) :: x_observation,y_observation,z_observation, &
+ g_x,g_y,g_z,G_xx,G_yy,G_zz,G_xy,G_xz,G_yz,temporary_array_for_sum
+
end module specfem_par
!=====================================================================
More information about the CIG-COMMITS
mailing list