[cig-commits] [commit] devel, master: clarified the definition of the reference ellipsoid and left some notes about this in the users manual and in the conversion routine (719b654)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:12:18 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 719b654b3f2fd8bc0cfe220b255a4528de9c3846
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed Apr 30 18:46:24 2014 +0200
clarified the definition of the reference ellipsoid and left some notes about this in the users manual and in the conversion routine
>---------------------------------------------------------------
719b654b3f2fd8bc0cfe220b255a4528de9c3846
doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf | Bin 24341435 -> 24343297 bytes
doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex | 11 +++++++++++
setup/constants.h.in | 24 ++++++++++++++++++++++++
src/shared/rthetaphi_xyz.f90 | 10 ++--------
4 files changed, 37 insertions(+), 8 deletions(-)
diff --git a/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf b/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf
index f9f0127..ee9e3ce 100644
Binary files a/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf and b/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf differ
diff --git a/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex b/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex
index 63bf3c8..72af062 100644
--- a/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex
+++ b/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex
@@ -780,6 +780,17 @@ is on.
\item [{\texttt{ELLIPTICITY}}] Set to \texttt{.true.} if the mesh should
make the Earth model elliptical in shape according to Clairaut's equation
\citep{DaTr98}. This feature adds no cost to the simulation.
+%
+From \cite{DaTr98}: "Spherically-symmetric Earth models all have the same hydrostatic surface ellipticity 1/299.8. This is 0.5 percent smaller than observed flattening of best-fitting ellipsoid 1/298.3. The discrepancy is referred to as the "excess equatorial bulge of the Earth", an early discovery of artificial satellite geodesy."
+%
+From Paul Melchior, IUGG General Assembly, Vienna, Austria, August 1991 Union lecture, available at http://www.agu.org/books: "It turns out that the spheroidal models constructed on the basis of the spherically-symmetric models (PREM, 1066A) by using the Clairaut differential equation to calculate the flattening
+in function of the radius vector imply hydrostaticity. These have surface ellipticity 1/299.8 and
+a corresponding dynamical flattening of .0033 (PREM). The actual ellipticty of the Earth for a best-fitting ellipsoid is 1/298.3 with a corresponding dynamical flattening of .0034."
+%
+Thus, flattening f = 1/299.8 is what is used in SPECFEM3D\_GLOBE, as it should. And thus eccentricity squared $e^2 = 1 - (1-f)^2 = 1 - (1 - 1/299.8)^2 = 0.00665998813529$, and the correction factor used in the code to convert geographic latitudes to geocentric is $1 - e^2 = (1-f)^2 = (1 - 1/299.8)^2 = 0.9933400118647$.
+%
+As a comparison, the classical World Geodetic System reference ellipsoid WGS 84 (see e.g. http://en.wikipedia.org/wiki/World\_Geodetic\_System) has $f = 1/298.2572236$.
+%
\item [{\texttt{TOPOGRAPHY}}] Set to \texttt{.true.} if topography and
bathymetry should be incorporated based upon model ETOPO5 \citep{Etopo5}.
This feature adds no cost to the simulation.
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 43e100b..9309125 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -131,6 +131,30 @@
!! pathname of the topography file (un-smoothed)
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/ETOPO1.xyz'
+! For the reference ellipsoid to convert geographic latitudes to geocentric:
+!
+! From Dahlen and Tromp (1998): "Spherically-symmetric Earth models all have the same hydrostatic
+! surface ellipticity 1/299.8. This is 0.5 percent smaller than observed flattening of best-fitting ellipsoid 1/298.3.
+! The discrepancy is referred to as the "excess equatorial bulge of the Earth",
+! an early discovery of artificial satellite geodesy."
+!
+! From Paul Melchior, IUGG General Assembly, Vienna, Austria, August 1991 Union lecture,
+! available at http://www.agu.org/books/sp/v035/SP035p0047/SP035p0047.pdf :
+! "It turns out that the spheroidal models constructed on the basis of the spherically-symmetric models (PREM, 1066A)
+! by using the Clairaut differential equation to calculate the flattening in function of the radius vector imply hydrostaticity.
+! These have surface ellipticity 1/299.8 and a corresponding dynamical flattening of .0033 (PREM).
+! The actual ellipticty of the Earth for a best-fitting ellipsoid is 1/298.3 with a corresponding dynamical flattening of .0034."
+!
+! Thus, flattening f = 1/299.8 is what is used in SPECFEM3D_GLOBE, as it should.
+! And thus eccentricity squared e^2 = 1 - (1-f)^2 = 1 - (1 - 1/299.8)^2 = 0.00665998813529,
+! and the correction factor used in the code to convert geographic latitudes to geocentric
+! is 1 - e^2 = (1-f)^2 = (1 - 1/299.8)^2 = 0.9933400118647.
+!
+! As a comparison, the classical World Geodetic System reference ellipsoid WGS 84
+! (see e.g. http://en.wikipedia.org/wiki/World_Geodetic_System) has f = 1/298.2572236.
+ double precision, parameter :: FLATTENING_F = 1.d0 / 299.8d0
+ double precision, parameter :: ONE_MINUS_F_SQUARED = (1.d0 - FLATTENING_F)**2
+
! Use GLL points to capture TOPOGRAPHY and ELLIPTICITY (experimental feature, currently does not work, DO NOT USE)
logical,parameter :: USE_GLL = .false.
diff --git a/src/shared/rthetaphi_xyz.f90 b/src/shared/rthetaphi_xyz.f90
index f56bbbc..821d43c 100644
--- a/src/shared/rthetaphi_xyz.f90
+++ b/src/shared/rthetaphi_xyz.f90
@@ -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_OLD_VERSION_5_1_5_FORMAT
+ use constants, only: PI_OVER_TWO,DEGREES_TO_RADIANS,ASSUME_PERFECT_SPHERE,ONE_MINUS_F_SQUARED
implicit none
@@ -291,13 +291,7 @@
if( .not. ASSUME_PERFECT_SPHERE ) then
! converts geographic (lat_prime) to geocentric latitude and converts to co-latitude (theta)
- 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
- ! note: factor 0.99333999308198295 = 1 - e**2 with eccentricity e**2 = 0.0066600069180170474
- theta = PI_OVER_TWO - atan( 0.99333999308198295d0*dtan(lat_prime * DEGREES_TO_RADIANS) )
- endif
+ theta = PI_OVER_TWO - atan( ONE_MINUS_F_SQUARED*dtan(lat_prime * DEGREES_TO_RADIANS) )
else
! for perfect sphere, geocentric and geographic latitudes are the same
! converts latitude (in degrees to co-latitude (in radians)
More information about the CIG-COMMITS
mailing list