[cig-commits] [commit] : Modified Vp-rho scaling to use that of Brocher et al. (2005). (756e7e1)

[cig_noreply at geodynamics.org]

Repository : ssh://geoshell/specfem3d

On branch  : 
Link       : https://github.com/geodynamics/specfem2d/compare/1e201257d91c794056b990a43329e05d04f77454...0000000000000000000000000000000000000000

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

commit 756e7e19dd1a14c6673129e37eeac085a080827a
Author: Carl Tape <carltape at gi.alaska.edu>
Date:   Thu Jun 3 14:41:26 2010 +0000

    Modified Vp-rho scaling to use that of Brocher et al. (2005).


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

756e7e19dd1a14c6673129e37eeac085a080827a
 compute_rho_estimate.f90 | 25 ++++++++++++++++++++-----
 constants.h.in           |  6 +++---
 2 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/compute_rho_estimate.f90 b/compute_rho_estimate.f90
index 5fb34b7..50ce76a 100644
--- a/compute_rho_estimate.f90
+++ b/compute_rho_estimate.f90
@@ -25,17 +25,32 @@
 
   subroutine compute_rho_estimate(rho,vp)
 
-! compute rho estimate in Gocad block and in Hauksson's model
-! based upon Vp
+! compute rho estimate in Gocad block and in Hauksson model from Vp
 
   implicit none
 
   include "constants.h"
 
   double precision rho,vp
-
-! scale density - use empirical rule from Christiane
-  rho = 0.33d0 * vp + 1280.d0
+  double precision :: P0, P1, P2, P3, P4, P5
+
+!!$! Vp-rho empirical rule from Stidham et al. (2001)
+!!$! -- used in Komatitsch et al. (2004) simulations
+!!$  P1 = 0.33d0 
+!!$  P0 = 1280.d0
+!!$  rho = P1*vp + P0
+!!$!  rho = 0.33d0 * vp + 1280.d0
+
+! Vp-rho empirical rule from Ludwig-Nafe-Drake (1970),
+! which is listed in Brocher (2005a)
+! -- used in Tape et al. (2009) simulations
+   P5 =  1.0600d-16
+   P4 = -4.3000d-12
+   P3 =  6.7100d-08
+   P2 = -4.7210d-04
+   P1 =  1.6612d0
+   P0 =  0.0d0
+   rho = P5*vp**5 + P4*vp**4 + P3*vp**3 + P2*vp**2 + P1*vp + P0
 
 ! make sure density estimate is reasonable
   if(rho > DENSITY_MAX) rho = DENSITY_MAX
diff --git a/constants.h.in b/constants.h.in
index 6ca6ca9..66c075a 100644
--- a/constants.h.in
+++ b/constants.h.in
@@ -74,9 +74,9 @@
 ! to avoid taking into account spurious oscillations in topography model
   double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 10.d0
 
-! min and max density in the model
-  double precision, parameter :: DENSITY_MAX = 3000.d0
-  double precision, parameter :: DENSITY_MIN = 2000.d0
+! min and max density in the model (based on Brocher, 2005, eq. 1)
+  double precision, parameter :: DENSITY_MAX = 3476.d0  ! 3000.d0
+  double precision, parameter :: DENSITY_MIN = 1635.d0  ! 2000.d0
 
 ! density of sea water
   real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0