[cig-commits] r22955 - in seismo/3D/SPECFEM3D/trunk: doc src/specfem3D
ampuero at geodynamics.org
ampuero at geodynamics.org
Wed Oct 9 00:18:53 PDT 2013
Author: ampuero
Date: 2013-10-09 00:18:53 -0700 (Wed, 09 Oct 2013)
New Revision: 22955
Modified:
seismo/3D/SPECFEM3D/trunk/doc/README_gravityPerturbation bis.txt
seismo/3D/SPECFEM3D/trunk/src/specfem3D/gravity_perturbation.f90
Log:
gravity perturbations: cleaned up module and its documentation, added time to output file
Modified: seismo/3D/SPECFEM3D/trunk/doc/README_gravityPerturbation bis.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/doc/README_gravityPerturbation bis.txt 2013-10-09 00:32:51 UTC (rev 22954)
+++ seismo/3D/SPECFEM3D/trunk/doc/README_gravityPerturbation bis.txt 2013-10-09 07:18:53 UTC (rev 22955)
@@ -4,10 +4,19 @@
modifications.
+Theory
+------
+
+The computation of gravity perturbations induced by deformation is based on
+equation 4 of:
+J. Harms, R. DeSalvo, S. Dorsher and V. Mandic (2009), "Simulation of underground
+gravity gradients from stochastic seismic fields", Phys. Rev. D, 80, 122001
+
+
Enabling gravity computations
-----------------------------
-Gravity perturbation can be computed for any SPECFEM3D simulation by
+Gravity perturbation can be computed in any SPECFEM3D simulation by
placing a file called "gravity_stations" in the DATA directory, in
addition to the regular input files.
@@ -15,7 +24,7 @@
Input file format
-----------------
-The format of the "gravity_stations" input file is as follows:
+The format of the "gravity_stations" input file is:
n dt_gap
x1 y1 z1
@@ -29,6 +38,16 @@
SPECFEM3D simulation
+Output file format
+-----------------
+
+Time series of gravity are output in files named "OUTPUT_FILES/stat*.grav",
+where * is the station index (one file per station). Their format is four
+columns:
+ t ax ay az
+(time and acceleration along x, y and z, respectively).
+
+
Code modifications
------------------
@@ -39,7 +58,7 @@
1. during the initialization,
2. during the iterative time stepping scheme and
3. at the output stage.
-In #1 the code checks for the presence of the input file gravity_stations.
+In #1 the code checks for the presence of the input file "gravity_stations".
If the file exists, the flag "GRAVITY_SIMULATION" is turned on and the
subroutines #2 and #3 are invoked.
@@ -48,4 +67,4 @@
------
Surendra Somala (Caltech) surendra at caltech.edu - 2013
-with advice from Jan Harms and Pablo Ampuero (ampuero at gps.caltech.edu)
\ No newline at end of file
+with advice from Jan Harms and Pablo Ampuero (ampuero at gps.caltech.edu)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/gravity_perturbation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/gravity_perturbation.f90 2013-10-09 00:32:51 UTC (rev 22954)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/gravity_perturbation.f90 2013-10-09 07:18:53 UTC (rev 22955)
@@ -42,8 +42,8 @@
double precision :: Jac3D
! coordinates of the control points
double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
- integer ia,iax,iay,iaz
- integer iaddx(NGNOD),iaddy(NGNOD),iaddz(NGNOD)
+ integer ia
+ integer, dimension(NGNOD) :: iaddx,iaddy,iaddz,iax,iay,iaz
integer nstep_grav
open(unit=IIN_G,file='../DATA/gravity_stations',status='old',iostat=ier)
@@ -92,53 +92,50 @@
call usual_hex_nodes(NGNOD,iaddx,iaddy,iaddz)
- do ispec=1,NSPEC_AB
- where( rho_vs(:,:,:,ispec) > TINYVAL )
- vs_elem(:,:,:) = mustore(:,:,:,ispec) / rho_vs(:,:,:,ispec)
- elsewhere
- vs_elem(:,:,:) = 0.0
- endwhere
+ ! define coordinates of the control points of the element
+ do ia=1,NGNOD
- rho_elem = rho_vs(:,:,:,ispec)/vs_elem
+ if(iaddx(ia) == 0) then
+ iax(ia) = 1
+ else if(iaddx(ia) == 1) then
+ iax(ia) = (NGLLX+1)/2
+ else if(iaddx(ia) == 2) then
+ iax(ia) = NGLLX
+ else
+ call exit_MPI(myrank,'incorrect value of iaddx')
+ endif
- ! define coordinates of the control points of the element
- do ia=1,NGNOD
+ if(iaddy(ia) == 0) then
+ iay(ia) = 1
+ else if(iaddy(ia) == 1) then
+ iay(ia) = (NGLLY+1)/2
+ else if(iaddy(ia) == 2) then
+ iay(ia) = NGLLY
+ else
+ call exit_MPI(myrank,'incorrect value of iaddy')
+ endif
- if(iaddx(ia) == 0) then
- iax = 1
- else if(iaddx(ia) == 1) then
- iax = (NGLLX+1)/2
- else if(iaddx(ia) == 2) then
- iax = NGLLX
- else
- call exit_MPI(myrank,'incorrect value of iaddx')
- endif
+ if(iaddz(ia) == 0) then
+ iaz(ia) = 1
+ else if(iaddz(ia) == 1) then
+ iaz(ia) = (NGLLZ+1)/2
+ else if(iaddz(ia) == 2) then
+ iaz(ia) = NGLLZ
+ else
+ call exit_MPI(myrank,'incorrect value of iaddz')
+ endif
- if(iaddy(ia) == 0) then
- iay = 1
- else if(iaddy(ia) == 1) then
- iay = (NGLLY+1)/2
- else if(iaddy(ia) == 2) then
- iay = NGLLY
- else
- call exit_MPI(myrank,'incorrect value of iaddy')
- endif
+ enddo
- if(iaddz(ia) == 0) then
- iaz = 1
- else if(iaddz(ia) == 1) then
- iaz = (NGLLZ+1)/2
- else if(iaddz(ia) == 2) then
- iaz = NGLLZ
- else
- call exit_MPI(myrank,'incorrect value of iaddz')
- endif
+ do ispec=1,NSPEC_AB
- iglob = ibool(iax,iay,iaz,ispec)
+ rho_elem = rho_vs(:,:,:,ispec)*rho_vs(:,:,:,ispec)/mustore(:,:,:,ispec)
+
+ do ia=1,NGNOD
+ iglob = ibool(iax(ia),iay(ia),iaz(ia),ispec)
xelm(ia) = dble(xstore(iglob))
yelm(ia) = dble(ystore(iglob))
zelm(ia) = dble(zstore(iglob))
-
enddo
do k = 1,NGLLZ
@@ -150,6 +147,7 @@
enddo
enddo
enddo
+
enddo
end subroutine gravity_init
@@ -329,7 +327,7 @@
subroutine gravity_output()
- use specfem_par, only : myrank,NPROC,NSTEP
+ use specfem_par, only : myrank,NPROC,NSTEP,DT
implicit none
integer :: isample,istat,nstep_grav
@@ -343,7 +341,7 @@
write(sisname,"('../OUTPUT_FILES/stat',I0,'.grav')") istat
open(unit=IOUT,file=sisname,status='replace')
do isample = 1,nstep_grav
- write(IOUT,*) accE(isample,istat),accN(isample,istat),accZ(isample,istat)
+ write(IOUT,*) isample*DT*ntimgap, accE(isample,istat),accN(isample,istat),accZ(isample,istat)
enddo
close(IOUT)
enddo
@@ -353,7 +351,7 @@
write(sisname,"('../OUTPUT_FILES/stat',I0,'.grav')") istat
open(unit=IOUT,file=sisname,status='replace')
do isample = 1,nstep_grav
- write(IOUT,*) accE(isample,istat),accN(isample,istat),accZ(isample,istat)
+ write(IOUT,*) isample*DT*ntimgap, accE(isample,istat),accN(isample,istat),accZ(isample,istat)
enddo
close(IOUT)
enddo
More information about the CIG-COMMITS
mailing list