[cig-commits] [commit] devel, master: adds scalar moment and moment magnitude information to source output (8c5d25d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:29:46 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 8c5d25d24d616fc8fc61643fc384077ad4c7ba49
Author: daniel peter <peterda at ethz.ch>
Date: Wed Sep 3 15:05:11 2014 +0200
adds scalar moment and moment magnitude information to source output
>---------------------------------------------------------------
8c5d25d24d616fc8fc61643fc384077ad4c7ba49
src/specfem3D/get_cmt.f90 | 107 +++++++++++++++++++++++++++++++++++++--
src/specfem3D/locate_sources.f90 | 16 ++++--
2 files changed, 116 insertions(+), 7 deletions(-)
diff --git a/src/specfem3D/get_cmt.f90 b/src/specfem3D/get_cmt.f90
index f895d87..7e12553 100644
--- a/src/specfem3D/get_cmt.f90
+++ b/src/specfem3D/get_cmt.f90
@@ -95,11 +95,15 @@
! debug
!print*,'line ----',string,'----'
- ! read header line with event information
+ ! reads header line with event information (assumes fixed format)
! old line: read(string,"(a4,i5,i3,i3,i3,i3,f6.2)") datasource,yr,mo,da,ho,mi,sec
- ! reads in chunks of the first line (which contain numbers)
- ! to get rid of the first datasource qualifyer string like PDE,PDEQ,.. which can have variable length)
+ ! reads header line with event information (free format)
+ ! gets rid of the first datasource qualifyer string which can have variable length, like:
+ ! "PDE 2014 9 3 .."
+ ! " PDEQ2014 9 3 .."
+ ! " MLI 1971 1 1 .."
+ ! note: globalcmt.org solutions might have missing spaces after datasource qualifier
!
! reads in year,month,day,hour,minutes,seconds
istart = 1
@@ -167,7 +171,7 @@
! checks time information
- if (yr <= 0 .or. yr > 10000) then
+ if (yr <= 0 .or. yr > 3000) then
write(IMAIN,*) 'Error reading year: ',yr,' in source ',isource,'is invalid'
stop 'Error reading year out of header line in CMTSOLUTION file'
endif
@@ -243,6 +247,8 @@
endif
read(string(7:len_trim(string)),*) depth(isource)
+ ! seismic moment tensor
+ ! CMTSOLUTION: components given in dyne-cm
! read Mrr
read(IIN,"(a)",iostat=ier) string
if (ier /= 0) then
@@ -370,3 +376,96 @@
end function
end subroutine get_cmt
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ double precision function get_cmt_scalar_moment(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! calculates scalar moment (M0)
+
+ use constants,only: RHOAV,R_EARTH,PI,GRAV
+
+ implicit none
+
+ double precision, intent(in) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! local parameters
+ double precision :: scalar_moment,scaleM
+
+ ! scalar moment:
+ ! see equation (1.4) in P.G. Silver and T.H. Jordan, 1982,
+ ! "Optimal estiamtion of scalar seismic moment",
+ ! Geophys. J.R. astr. Soc., 70, 755 - 787
+ !
+ ! or see equation (5.91) in Dahlen & Tromp (1998)
+ !
+ ! moment tensor M is a symmetric 3x3 tensor, and has six independent components
+ !
+ ! the euclidean matrix norm is invariant under rotation.
+ ! thus, input can be:
+ ! Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! or
+ ! Mrr,Mtt,Mpp,Mrt,Mrp,Mtp
+ !
+ ! euclidean (or Frobenius) norm of a matrix: M0**2 = sum( Mij**2 )
+ scalar_moment = Mxx**2 + Myy**2 + Mzz**2 + 2.d0 * Mxy**2 + 2.d0 * Mxz**2 + 2.d0 * Myz**2
+
+ ! adds 1/2 to be coherent with double couple or point sources
+ scalar_moment = dsqrt(scalar_moment/2.0d0)
+
+ ! note: moment tensor is non-dimensionalized
+ !
+ ! re-adds scale factor for the moment tensor
+ ! CMTSOLUTION file values are in dyne.cm
+ ! 1 dyne is 1 gram * 1 cm / (1 second)^2
+ ! 1 Newton is 1 kg * 1 m / (1 second)^2
+ ! thus 1 Newton = 100,000 dynes
+ ! therefore 1 dyne.cm = 1e-7 Newton.m
+ scaleM = 1.d7 * RHOAV * (R_EARTH**5) * PI * GRAV * RHOAV
+
+ ! return value (in dyne-cm)
+ get_cmt_scalar_moment = scalar_moment * scaleM
+
+ end function
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ double precision function get_cmt_moment_magnitude(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! calculates scalar moment (M0)
+
+ implicit none
+
+ double precision, intent(in) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! local parameters
+ double precision :: M0,Mw
+ double precision,external :: get_cmt_scalar_moment
+
+ ! scalar moment
+ M0 = get_cmt_scalar_moment(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! moment magnitude by Hanks & Kanamori, 1979
+ ! Mw = 2/3 log( M0 ) - 10.7 (dyne-cm)
+ !
+ ! alternative forms:
+ ! Mw = 2/3 ( log( M0 ) - 16.1 ) (N-m) "moment magnitude" by Hanks & Kanamori(1979) or "energy magnitude" by Kanamori (1977)
+ !
+ ! Aki & Richards ("Quantitative Seismology",2002):
+ ! Mw = 2/3 ( log( M0 ) - 9.1 ) (N-m)
+ !
+ ! conversion: dyne-cm = 10**-7 N-m
+ !
+ ! we follow here the USGS magnitude policy:
+ ! "All USGS statements of moment magnitude should use M = (log M0)/1.5-10.7
+ ! for converting from scalar moment M0 to moment magnitude. (..)"
+ ! see: http://earthquake.usgs.gov/aboutus/docs/020204mag_policy.php
+
+ Mw = 2.d0/3.d0 * log10( M0 ) - 10.7
+
+ ! return value
+ get_cmt_moment_magnitude = Mw
+
+ end function
diff --git a/src/specfem3D/locate_sources.f90 b/src/specfem3D/locate_sources.f90
index e2df265..9dc1e0c 100644
--- a/src/specfem3D/locate_sources.f90
+++ b/src/specfem3D/locate_sources.f90
@@ -124,6 +124,9 @@
integer :: imin,imax,jmin,jmax,kmin,kmax
double precision :: f0,t0_ricker
+ double precision, external :: get_cmt_scalar_moment
+ double precision, external :: get_cmt_moment_magnitude
+
! mask source region (mask values are between 0 and 1, with 0 around sources)
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: mask_source
@@ -665,6 +668,13 @@
write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
endif
write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
+ write(IMAIN,*)
+ write(IMAIN,*) 'magnitude of the source:'
+ write(IMAIN,*) ' scalar moment M0 = ', &
+ get_cmt_scalar_moment(Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource)),' dyne-cm'
+ write(IMAIN,*) ' moment magnitude Mw = ', &
+ get_cmt_moment_magnitude(Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource))
+ write(IMAIN,*)
! writes out actual source position to VTK file
write(IOUT_VTK,'(3e18.6)') sngl(x_found_source(isource_in_this_subset)), &
@@ -896,6 +906,7 @@
double precision, external :: comp_source_time_function,comp_source_spectrum
double precision, external :: comp_source_time_function_rickr
+ double precision, external :: get_cmt_scalar_moment
character(len=150) :: plot_file
@@ -925,9 +936,8 @@
status='unknown',iostat=ier)
if (ier /= 0 ) call exit_mpi(0,'Error opening plot_source_time_function file')
- scalar_moment = Mxx(isource)**2 + Myy(isource)**2 + Mzz(isource)**2 &
- + Mxy(isource)**2 + Mxz(isource)**2 + Myz(isource)**2
- scalar_moment = dsqrt(scalar_moment/2.0d0)
+ ! calculates scalar moment M0
+ scalar_moment = get_cmt_scalar_moment(Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource))
! define t0 as the earliest start time
! note: this calculation here is only used for outputting the plot_source_time_function file
More information about the CIG-COMMITS
mailing list