[cig-commits] [commit] devel, master: added a comment referring to eq (14.21) in Dahlen and Tromp (1998) regarding ellipticity (49bb85f)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem3d_globe

On branches: devel,master
Link       : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f

>---------------------------------------------------------------

commit 49bb85f785105dcec155da95898bd6c8e089a5ac
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date:   Wed May 7 22:55:56 2014 +0200

    added a comment referring to eq (14.21) in Dahlen and Tromp (1998) regarding ellipticity


>---------------------------------------------------------------

49bb85f785105dcec155da95898bd6c8e089a5ac
 src/auxiliaries/combine_vol_data_shared.f90 | 20 +++++++++-----------
 src/meshfem3D/finalize_mesher.f90           |  2 +-
 src/shared/make_ellipticity.f90             | 11 +++--------
 3 files changed, 13 insertions(+), 20 deletions(-)

diff --git a/src/auxiliaries/combine_vol_data_shared.f90 b/src/auxiliaries/combine_vol_data_shared.f90
index 57dc6b3..3f07dfd 100644
--- a/src/auxiliaries/combine_vol_data_shared.f90
+++ b/src/auxiliaries/combine_vol_data_shared.f90
@@ -25,9 +25,6 @@
 !
 !=====================================================================
 
-!
-! ------------------------------------------------------------------------------------------------
-!
   subroutine reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
 
   use constants
@@ -184,25 +181,26 @@
 
   do i=2,NR
     call intgrl(i_rho,r,1,i,rho,s1,s2,s3)
+! Radau approximation of Clairaut's equation for first-order terms of ellipticity, see e.g. Jeffreys H.,
+! The figures of rotating planets, Mon. Not. R. astr. Soc., vol. 113, p. 97-105 (1953).
+! The Radau approximation is mentioned on page 97.
+! For more details see Section 14.1.2 in Dahlen and Tromp (1998)
+! (see also in file ellipticity_equations_from_Dahlen_Tromp_1998.pdf in the "doc" directory of the code).
     call intgrl(i_radau,r,1,i,radau,s1,s2,s3)
     z=(2.0d0/3.0d0)*i_radau/(i_rho*r(i)*r(i))
+! this comes from equation (14.19) in Dahlen and Tromp (1998)
     eta(i)=(25.0d0/4.0d0)*((1.0d0-(3.0d0/2.0d0)*z)**2.0d0)-1.0d0
     k(i)=eta(i)/(r(i)**3.0d0)
   enddo
 
-  g_a=4.0D0*i_rho
-
   bom=TWO_PI/(HOURS_PER_DAY*SECONDS_PER_HOUR)
 
 ! non-dimensionalized version
   bom=bom/sqrt(PI*GRAV*RHOAV)
 
-!! DK DK I think 24.d0 below stands for HOURS_PER_DAY and thus I replace it here
-!! DK DK in order to be consistent if someone uses the code one day for other planets
-!! DK DK with a different rotation rate, or for the Earth in the past of in the future i.e. with a different rate as well.
-!! DK DK Please do not hesitate to fix it back if my assumption below was wrong.
-!! DK DK  epsilonval(NR)=15.0d0*(bom**2.0d0)/(24.0d0*i_rho*(eta(NR)+2.0d0))
-  epsilonval(NR)=15.0d0*(bom**2.0d0)/(HOURS_PER_DAY*i_rho*(eta(NR)+2.0d0))
+  g_a = 4.0d0*i_rho
+! this is the equation right above (14.21) in Dahlen and Tromp (1998)
+  epsilonval(NR) = (5.0d0/2.d0)*(bom**2.0d0)*R_UNIT_SPHERE / (g_a * (eta(NR)+2.0d0))
 
   do i=1,NR-1
     call intgrl(exponentval,r,i,NR,k,s1,s2,s3)
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index 5870b77..79fde9d 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -112,7 +112,7 @@
     if(NCHUNKS == 6) then
       write(IMAIN,*)
       write(IMAIN,*) 'computed total Earth mass for this density model and mesh: ',Earth_mass_total,' kg'
-      write(IMAIN,*) '   (should be not too far from 5.97E+24 kg)'
+      write(IMAIN,*) '   (should be not too far from 5.974E+24 kg)'
       write(IMAIN,*)
       ! take into account the fact that dimensions have been non-dimensionalized by dividing them by R_EARTH
       write(IMAIN,*) 'average density for this density model and mesh: ',Earth_mass_total / (volume_total * R_EARTH**3),' kg/m3'
diff --git a/src/shared/make_ellipticity.f90 b/src/shared/make_ellipticity.f90
index 6766fdf..4b9aa51 100644
--- a/src/shared/make_ellipticity.f90
+++ b/src/shared/make_ellipticity.f90
@@ -155,19 +155,14 @@
     k(i)=eta(i)/(r(i)**3.0d0)
   enddo
 
-  g_a=4.0D0*i_rho
-
   bom=TWO_PI/(HOURS_PER_DAY*SECONDS_PER_HOUR)
 
 ! non-dimensionalized value
   bom=bom/sqrt(PI*GRAV*RHOAV)
 
-!! DK DK I think 24.d0 below stands for HOURS_PER_DAY and thus I replace it here
-!! DK DK in order to be consistent if someone uses the code one day for other planets
-!! DK DK with a different rotation rate, or for the Earth in the past of in the future i.e. with a different rate as well.
-!! DK DK Please do not hesitate to fix it back if my assumption below was wrong.
-!! DK DK  epsilonval(NR)=15.0d0*(bom**2.0d0)/(24.0d0*i_rho*(eta(NR)+2.0d0))
-  epsilonval(NR)=15.0d0*(bom**2.0d0)/(HOURS_PER_DAY*i_rho*(eta(NR)+2.0d0))
+  g_a = 4.0d0*i_rho
+! this is the equation right above (14.21) in Dahlen and Tromp (1998)
+  epsilonval(NR) = (5.0d0/2.d0)*(bom**2.0d0)*R_UNIT_SPHERE / (g_a * (eta(NR)+2.0d0))
 
   do i=1,NR-1
     call intgrl(exponentval,r,i,NR,k,s1,s2,s3)