updates model_prem_iso() routine (for backward version compatibility to match gravity)

Author: danielpeter

src/shared/model_prem.f90

@@ -114,7 +114,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 
   ! local parameters
   double precision :: r,scaleval
-  double precision :: RICB
+  double precision :: ROCEAN,RMIDDLE_CRUST,RMOHO,R80,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB
 
 ! compute real physical radius in meters
   r = x * R_PLANET
@@ -124,11 +124,42 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 
 ! note: using stop statements, not exit_mpi() calls to avoid the need for MPI libraries when linking xcreate_header_file
 
-  ! version compatibility
+  ! default setting
+  ROCEAN = PREM_ROCEAN
+  RMIDDLE_CRUST = PREM_RMIDDLE_CRUST
+  RMOHO = PREM_RMOHO
+  R80 = PREM_R80
+  R220 = PREM_R220
+  R400 = PREM_R400
+  R600 = PREM_R600
+  R670 = PREM_R670
+  R771 = PREM_R771
+  RTOPDDOUBLEPRIME = PREM_RTOPDDOUBLEPRIME
+  RCMB = PREM_RCMB
+  RICB = PREM_RICB
+
+  ! note: different versions use different radii when calling this routine
   if (USE_OLD_VERSION_FORMAT) then
+    ! version compatibility with older versions 7.0 - 8.0
     RICB = PREM_RICB_OLD
-  else
-    RICB = PREM_RICB
+
+    ! version compatibility with special case Berkeley model
+    ! note: old Berkeley implementation used radii different to PREM when calling this routine for gravity evaluation.
+    !       this was inconsistent with the the spline evaluation in the old make_gravity.f90.
+    !       the change here below is to re-introduce this inconsistency to match the old Berkeley model output.
+    !
+    !       new(er) versions use the "original" PREM model for calculating gravity (and ellipticity).
+    !       that is, gravity (acceleration g) and ellipticity (based on density model) is based on PREM, since we know that
+    !       PREM was constructed with corresponding constraints.
+    !       other reference models might better adapt to velocities, but might violate for example total mass constraints.
+    !       the different radii of the reference models are used for meshing purposes only.
+    !
+    if (REFERENCE_1D_MODEL == REFERENCE_MODEL_SEMUCB) then
+      RICB = PREM_RICB
+      RMOHO = 6341000.d0
+      R400 = 5961000.d0
+      R670 = 5721000.d0
+    endif
   endif
 
   if (check_doubling_flag) then
@@ -145,31 +176,31 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
     !
     !--- outer core
     !
-    else if (r > RICB .and. r < PREM_RCMB) then
+    else if (r > RICB .and. r < RCMB) then
       if (idoubling /= IFLAG_OUTER_CORE_NORMAL) &
         stop 'wrong doubling flag for outer core point in model_prem_iso()'
     !
     !--- D" at the base of the mantle
     !
-    else if (r > PREM_RCMB .and. r < PREM_RTOPDDOUBLEPRIME) then
+    else if (r > RCMB .and. r < RTOPDDOUBLEPRIME) then
       if (idoubling /= IFLAG_MANTLE_NORMAL) &
         stop 'wrong doubling flag for D" point in model_prem_iso()'
     !
     !--- mantle: from top of D" to d670
     !
-    else if (r > PREM_RTOPDDOUBLEPRIME .and. r < PREM_R670) then
+    else if (r > RTOPDDOUBLEPRIME .and. r < R670) then
       if (idoubling /= IFLAG_MANTLE_NORMAL) &
         stop 'wrong doubling flag for top D" to d670 point in model_prem_iso()'
     !
     !--- mantle: from d670 to d220
     !
-    else if (r > PREM_R670 .and. r < PREM_R220) then
+    else if (r > R670 .and. r < R220) then
       if (idoubling /= IFLAG_670_220) &
         stop 'wrong doubling flag for d670 to d220 point in model_prem_iso()'
     !
     !--- mantle and crust: from d220 to MOHO and then to surface
     !
-    else if (r > PREM_R220) then
+    else if (r > R220) then
       if (idoubling /= IFLAG_220_80 .and. idoubling /= IFLAG_80_MOHO .and. idoubling /= IFLAG_CRUST) &
         stop 'wrong doubling flag for d220 to Moho to surface point in model_prem_iso()'
     endif
@@ -201,7 +232,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 !
 !--- outer core
 !
-  else if (r > RICB .and. r <= PREM_RCMB) then
+  else if (r > RICB .and. r <= RCMB) then
     drhodr = -1.2638d0 - 2.0d0*3.6426d0*x - 3.0d0*5.5281d0*x*x
     rho = 12.5815d0 - 1.2638d0*x - 3.6426d0*x*x - 5.5281d0*x*x*x
     if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then
@@ -223,7 +254,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 !
 !--- D" at the base of the mantle
 !
-  else if (r > PREM_RCMB .and. r <= PREM_RTOPDDOUBLEPRIME) then
+  else if (r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
     drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x
     rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x
     if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then
@@ -239,7 +270,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 !
 !--- mantle: from top of D" to d670
 !
-  else if (r > PREM_RTOPDDOUBLEPRIME .and. r <= PREM_R771) then
+  else if (r > RTOPDDOUBLEPRIME .and. r <= R771) then
     drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x
     rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x
     if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then
@@ -258,7 +289,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
     vs = 11.1671d0 - 13.7818d0*x + 17.4575d0*x*x - 9.2777d0*x*x*x
     Qmu = 312.0d0
     Qkappa = 57827.0d0
-  else if (r > PREM_R771 .and. r <= PREM_R670) then
+  else if (r > R771 .and. r <= R670) then
     drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x
     rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x
     vp = 29.2766d0 - 23.6027d0*x + 5.5242d0*x*x - 2.5514d0*x*x*x
@@ -268,28 +299,28 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
 !
 !--- mantle: above d670
 !
-  else if (r > PREM_R670 .and. r <= PREM_R600) then
+  else if (r > R670 .and. r <= R600) then
     drhodr = -1.4836d0
     rho = 5.3197d0 - 1.4836d0*x
     vp = 19.0957d0 - 9.8672d0*x
     vs = 9.9839d0 - 4.9324d0*x
     Qmu = 143.0d0
     Qkappa = 57827.0d0
-  else if (r > PREM_R600 .and. r <= PREM_R400) then
+  else if (r > R600 .and. r <= R400) then
     drhodr = -8.0298d0
     rho = 11.2494d0 - 8.0298d0*x
     vp = 39.7027d0 - 32.6166d0*x
     vs = 22.3512d0 - 18.5856d0*x
     Qmu = 143.0d0
     Qkappa = 57827.0d0
-  else if (r > PREM_R400 .and. r <= PREM_R220) then
+  else if (r > R400 .and. r <= R220) then
     drhodr = -3.8045d0
     rho = 7.1089d0 - 3.8045d0*x
     vp = 20.3926d0 - 12.2569d0*x
     vs = 8.9496d0 - 4.4597d0*x
     Qmu = 143.0d0
     Qkappa = 57827.0d0
-  else if (r > PREM_R220 .and. r <= PREM_R80) then
+  else if (r > R220 .and. r <= R80) then
     drhodr = 0.6924d0
     rho = 2.6910d0 + 0.6924d0*x
     vp = 4.1875d0 + 3.9382d0*x
@@ -299,7 +330,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
   else
     if (CRUSTAL .and. .not. SUPPRESS_CRUSTAL_MESH) then
 ! fill with PREM mantle and later add CRUST2.0
-      if (r > PREM_R80) then
+      if (r > R80) then
         ! density/velocity from mantle just below moho
         drhodr = 0.6924d0
         rho = 2.6910d0 + 0.6924d0*x
@@ -311,7 +342,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
       endif
     else
 ! use PREM crust
-      if (r > PREM_R80 .and. r <= PREM_RMOHO) then
+      if (r > R80 .and. r <= RMOHO) then
         drhodr = 0.6924d0
         rho = 2.6910d0 + 0.6924d0*x
         vp = 4.1875d0 + 3.9382d0*x
@@ -322,13 +353,13 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
       else if (SUPPRESS_CRUSTAL_MESH) then
         ! extend the Moho up to the surface instead of the crust
         drhodr = 0.6924d0
-        rho = 2.6910d0 + 0.6924d0*(PREM_RMOHO / R_PLANET)
-        vp = 4.1875d0 + 3.9382d0*(PREM_RMOHO / R_PLANET)
-        vs = 2.1519d0 + 2.3481d0*(PREM_RMOHO / R_PLANET)
+        rho = 2.6910d0 + 0.6924d0*(RMOHO / R_PLANET)
+        vp = 4.1875d0 + 3.9382d0*(RMOHO / R_PLANET)
+        vs = 2.1519d0 + 2.3481d0*(RMOHO / R_PLANET)
         Qmu = 600.0d0
         Qkappa = 57827.0d0
 
-      else if (r > PREM_RMOHO .and. r <= PREM_RMIDDLE_CRUST) then
+      else if (r > RMOHO .and. r <= RMIDDLE_CRUST) then
         drhodr = 0.0d0
         rho = 2.9d0
         vp = 6.8d0
@@ -346,15 +377,15 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_
           Qkappa = 57827.0d0
         endif
 
-      else if (r > PREM_RMIDDLE_CRUST .and. r <= PREM_ROCEAN) then
+      else if (r > RMIDDLE_CRUST .and. r <= ROCEAN) then
         drhodr = 0.0d0
         rho = 2.6d0
         vp = 5.8d0
         vs = 3.2d0
         Qmu = 600.0d0
         Qkappa = 57827.0d0
 ! for density profile for gravity, we do not check that r <= R_PLANET
-      else if (r > PREM_ROCEAN) then
+      else if (r > ROCEAN) then
         drhodr = 0.0d0
         rho = 2.6d0
         vp = 5.8d0

src/specfem3D/prepare_gravity.f90

@@ -230,6 +230,20 @@ subroutine prepare_gravity()
       endif
     enddo
 
+    ! debug - file output
+    if (.false.) then
+      open(21,file='tmp_grav.dat',status='unknown')
+      write(21,*) "# int_radius  # radius  # rho  # minus_g  # minus_dg"
+      do int_radius = 1,NRAD_GRAVITY
+        radius = r(int_radius)
+        rho = density_table(int_radius)
+        minus_g = minus_gravity_table(int_radius)
+        minus_dg = minus_deriv_gravity_table(int_radius)
+        write(21,*) int_radius,radius,rho,minus_g,minus_dg
+      enddo
+      close(21)
+    endif
+
     ! make sure fluid array is only assigned in outer core between 1222 and 3478 km
     ! lookup table is defined every 100 m
     do int_radius = 1,NRAD_GRAVITY
@@ -320,9 +334,9 @@ subroutine prepare_gravity()
       fac = minus_rho_g_over_kappa_fluid(int_radius)
 
       ! gravitational acceleration (integrated and multiply by rho / Kappa) in Cartesian coordinates
-      gxl = dsin(theta) * dcos(phi) * fac
-      gyl = dsin(theta) * dsin(phi) * fac
-      gzl = dcos(theta) * fac
+      gxl = (dsin(theta) * dcos(phi)) * fac
+      gyl = (dsin(theta) * dsin(phi)) * fac
+      gzl = (dcos(theta)) * fac
       gravity_pre_store_outer_core(1,iglob) = real(gxl,kind=CUSTOM_REAL)
       gravity_pre_store_outer_core(2,iglob) = real(gyl,kind=CUSTOM_REAL)
       gravity_pre_store_outer_core(3,iglob) = real(gzl,kind=CUSTOM_REAL)
@@ -437,9 +451,9 @@ subroutine prepare_gravity()
       fac = d_ln_density_dr_table(int_radius)
 
       ! gradient of d ln(rho)/dr in Cartesian coordinates
-      gxl = dsin(theta) * dcos(phi) * fac
-      gyl = dsin(theta) * dsin(phi) * fac
-      gzl = dcos(theta) * fac
+      gxl = (dsin(theta) * dcos(phi)) * fac
+      gyl = (dsin(theta) * dsin(phi)) * fac
+      gzl = (dcos(theta)) * fac
       gravity_pre_store_outer_core(1,iglob) = real(gxl,kind=CUSTOM_REAL)
       gravity_pre_store_outer_core(2,iglob) = real(gyl,kind=CUSTOM_REAL)
       gravity_pre_store_outer_core(3,iglob) = real(gzl,kind=CUSTOM_REAL)
@@ -546,9 +560,9 @@ subroutine prepare_gravity()
 
       ! Cartesian components of the gravitational acceleration
       ! multiplied by common factor rho
-      gxl = minus_g * sin_theta * cos_phi * rho
-      gyl = minus_g * sin_theta * sin_phi * rho
-      gzl = minus_g * cos_theta * rho
+      gxl = (minus_g * sin_theta * cos_phi) * rho
+      gyl = (minus_g * sin_theta * sin_phi) * rho
+      gzl = (minus_g * cos_theta) * rho
       gravity_pre_store_crust_mantle(1,iglob) = real(gxl,kind=CUSTOM_REAL)
       gravity_pre_store_crust_mantle(2,iglob) = real(gyl,kind=CUSTOM_REAL)
       gravity_pre_store_crust_mantle(3,iglob) = real(gzl,kind=CUSTOM_REAL)
@@ -693,9 +707,9 @@ subroutine prepare_gravity()
 
       ! Cartesian components of the gravitational acceleration
       ! multiplied by common factor rho
-      gxl = minus_g * sin_theta * cos_phi * rho
-      gyl = minus_g * sin_theta * sin_phi * rho
-      gzl = minus_g * cos_theta * rho
+      gxl = (minus_g * sin_theta * cos_phi) * rho
+      gyl = (minus_g * sin_theta * sin_phi) * rho
+      gzl = (minus_g * cos_theta) * rho
 
       gravity_pre_store_inner_core(1,iglob) = real(gxl,kind=CUSTOM_REAL)
       gravity_pre_store_inner_core(2,iglob) = real(gyl,kind=CUSTOM_REAL)