[cig-commits] [commit] devel, master: added constant SECONDS_PER_HOUR = 3600.d0 (f0ab85a)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:12:07 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit f0ab85a89fe59ceed13a4c0726a60a7dc50c86e6
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed Apr 30 17:06:07 2014 +0200
added constant SECONDS_PER_HOUR = 3600.d0
>---------------------------------------------------------------
f0ab85a89fe59ceed13a4c0726a60a7dc50c86e6
setup/constants.h.in | 1 +
src/auxiliaries/combine_vol_data_shared.f90 | 13 +++++++++++--
src/meshfem3D/create_mass_matrices.f90 | 8 ++++----
src/meshfem3D/model_atten3D_QRFSI12.f90 | 8 +++++---
src/shared/make_ellipticity.f90 | 15 ++++++++++++---
src/specfem3D/prepare_timerun.f90 | 12 ++++++------
6 files changed, 39 insertions(+), 18 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 67c7d07..43e100b 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -597,6 +597,7 @@
! number of hours per day for rotation rate of the Earth
double precision, parameter :: HOURS_PER_DAY = 24.d0
+ double precision, parameter :: SECONDS_PER_HOUR = 3600.d0
! for lookup table for gravity every 100 m in radial direction of Earth model
integer, parameter :: NRAD_GRAVITY = 70000
diff --git a/src/auxiliaries/combine_vol_data_shared.f90 b/src/auxiliaries/combine_vol_data_shared.f90
index a0e90ee..90af0b4 100644
--- a/src/auxiliaries/combine_vol_data_shared.f90
+++ b/src/auxiliaries/combine_vol_data_shared.f90
@@ -189,9 +189,18 @@
enddo
g_a=4.0D0*i_rho
- bom=TWO_PI/(24.0d0*3600.0d0)
+
+ bom=TWO_PI/(HOURS_PER_DAY*SECONDS_PER_HOUR)
+
+! non-dimensionalized version
bom=bom/sqrt(PI*GRAV*RHOAV)
- epsilonval(NR)=15.0d0*(bom**2.0d0)/(24.0d0*i_rho*(eta(NR)+2.0d0))
+
+!! 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))
do i=1,NR-1
call intgrl(exponentval,r,i,NR,k,s1,s2,s3)
diff --git a/src/meshfem3D/create_mass_matrices.f90 b/src/meshfem3D/create_mass_matrices.f90
index a539435..88a7928 100644
--- a/src/meshfem3D/create_mass_matrices.f90
+++ b/src/meshfem3D/create_mass_matrices.f90
@@ -269,15 +269,15 @@
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- two_omega_earth_dt = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat)
+ two_omega_earth_dt = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv) * deltat)
else
- two_omega_earth_dt = 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat
+ two_omega_earth_dt = 2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv) * deltat
endif
if(CUSTOM_REAL == SIZE_REAL) then
- b_two_omega_earth_dt = - sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat)
+ b_two_omega_earth_dt = - sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv) * deltat)
else
- b_two_omega_earth_dt = - 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat
+ b_two_omega_earth_dt = - 2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv) * deltat
endif
! definition depends if region is fluid or solid
diff --git a/src/meshfem3D/model_atten3D_QRFSI12.f90 b/src/meshfem3D/model_atten3D_QRFSI12.f90
index f8ac371..8b62bb0 100644
--- a/src/meshfem3D/model_atten3D_QRFSI12.f90
+++ b/src/meshfem3D/model_atten3D_QRFSI12.f90
@@ -205,8 +205,9 @@
real(kind=4) :: shdep(NSQ)
real(kind=4) :: wk1(NSQ),wk2(NSQ),wk3(NSQ)
real(kind=4) :: xlmvec(NSQ**2)
- double precision, parameter :: rmoho_prem = 6371.0-24.4
- double precision, parameter :: rcmb = 3480.0
+
+ double precision, parameter :: rmoho_prem = 6371.0d0 - 24.4d0
+ double precision, parameter :: rcmb = 3480.0d0
! in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
! to agree with the rest of specfem
@@ -265,7 +266,8 @@
smallq_ref=QRFSI12_Q_refqmu(ifnd)
smallq = smallq_ref
- if(depth < 650.d0) then !Colleen's model is only defined between depths of 24.4 and 650km
+! Colleen's model is only defined between depths of 24.4 and 650km
+ if(depth < 650.d0) then
do j=1,NSQ
shdep(j)=0.
enddo
diff --git a/src/shared/make_ellipticity.f90 b/src/shared/make_ellipticity.f90
index df0f60c..bbe960c 100644
--- a/src/shared/make_ellipticity.f90
+++ b/src/shared/make_ellipticity.f90
@@ -30,7 +30,7 @@
! creates a spline for the ellipticity profile in PREM
! radius and density are non-dimensional
- use constants,only: NR,TWO_PI,PI,GRAV,RHOAV
+ use constants,only: NR,TWO_PI,PI,GRAV,RHOAV,HOURS_PER_DAY,SECONDS_PER_HOUR
implicit none
@@ -150,9 +150,18 @@
enddo
g_a=4.0D0*i_rho
- bom=TWO_PI/(24.0d0*3600.0d0)
+
+ bom=TWO_PI/(HOURS_PER_DAY*SECONDS_PER_HOUR)
+
+! non-dimensionalized value
bom=bom/sqrt(PI*GRAV*RHOAV)
- epsilonval(NR)=15.0d0*(bom**2.0d0)/(24.0d0*i_rho*(eta(NR)+2.0d0))
+
+!! 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))
do i=1,NR-1
call intgrl(exponentval,r,i,NR,k,s1,s2,s3)
diff --git a/src/specfem3D/prepare_timerun.f90 b/src/specfem3D/prepare_timerun.f90
index 0230905..0a6030e 100644
--- a/src/specfem3D/prepare_timerun.f90
+++ b/src/specfem3D/prepare_timerun.f90
@@ -700,23 +700,23 @@
! distinguish between single and double precision for reals
if (SIMULATION_TYPE == 1) then
if(CUSTOM_REAL == SIZE_REAL) then
- two_omega_earth = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv))
+ two_omega_earth = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv))
else
- two_omega_earth = 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv)
+ two_omega_earth = 2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv)
endif
else
if(CUSTOM_REAL == SIZE_REAL) then
- two_omega_earth = - sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv))
+ two_omega_earth = - sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv))
else
- two_omega_earth = - 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv)
+ two_omega_earth = - 2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv)
endif
endif
if (SIMULATION_TYPE == 3) then
if(CUSTOM_REAL == SIZE_REAL) then
- b_two_omega_earth = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv))
+ b_two_omega_earth = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv))
else
- b_two_omega_earth = 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv)
+ b_two_omega_earth = 2.d0 * TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR * scale_t_inv)
endif
endif
else
More information about the CIG-COMMITS
mailing list