[cig-commits] [commit] devel, master: updates geocentric_2_geographic_dble() routine to account for new flattening constant (957c705)

[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 957c7059aa194f3c8db9d94c0fdd7a653b749bfe
Author: daniel peter <peterda at ethz.ch>
Date:   Thu Sep 18 13:03:02 2014 +0200

    updates geocentric_2_geographic_dble() routine to account for new flattening constant


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

957c7059aa194f3c8db9d94c0fdd7a653b749bfe
 src/shared/rthetaphi_xyz.f90 | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/shared/rthetaphi_xyz.f90 b/src/shared/rthetaphi_xyz.f90
index f635afe..a82246a 100644
--- a/src/shared/rthetaphi_xyz.f90
+++ b/src/shared/rthetaphi_xyz.f90
@@ -308,13 +308,20 @@
 ! Anyway, I'm sorry that I'm not giving a clear answer, but hopefully this
 ! gives some thoughts.
 
-  use constants, only: PI_OVER_TWO,TINYVAL,ASSUME_PERFECT_SPHERE,USE_OLD_VERSION_5_1_5_FORMAT
+  use constants, only: PI_OVER_TWO,TINYVAL,ASSUME_PERFECT_SPHERE,USE_OLD_VERSION_5_1_5_FORMAT,ONE_MINUS_F_SQUARED
 
   implicit none
 
   double precision,intent(in) :: theta
   double precision,intent(out) :: theta_prime
 
+  ! note: september, 2014
+  ! factor: 1/(1 - e^2) = 1/(1 - (1 - (1-f)^2)) = 1/( (1-f)^2 )
+  !         with eccentricity e^2 = 1 - (1-f)^2
+  ! see about Earth flattening in constants.h: flattening factor changed to 1/299.8
+  !                                            f = 1/299.8 -> 1/( (1-f)^2 ) = 1.0067046409645724
+  double precision, parameter :: FACTOR_TAN = 1.d0 / ONE_MINUS_F_SQUARED
+
   ! note: instead of 1/tan(theta) we take cos(theta)/sin(theta) and avoid division by zero
 
   if (.not. ASSUME_PERFECT_SPHERE) then
@@ -323,7 +330,7 @@
       theta_prime = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
     else
       ! converts geocentric colatitude theta to geographic colatitude theta_prime
-      theta_prime = PI_OVER_TWO - datan(1.00670466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
+      theta_prime = PI_OVER_TWO - datan(FACTOR_TAN*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
     endif
   else
     ! mesh is spherical, thus geocentric and geographic colatitudes are identical