Rescaled pressure in find_eta routines

Rescaled the pressure used to calculate density integrals in find_eta_3d and find_eta_2d to [R L2 T-2] and used the new pres_scale and SV_scale arguments to int_specific_vol_dp. All answers are bitwise identical.
NOAA-GFDL · Apr 6, 2020 · 94f94ed · 94f94ed
1 parent 3b3c34a
commit 94f94ed
Showing 1 changed file with 17 additions and 12 deletions.
diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90
@@ -47,12 +47,13 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
              !! the units of eta to m; by default this is US%Z_to_m.
 
   ! Local variables
-  real :: p(SZI_(G),SZJ_(G),SZK_(G)+1)    ! Hydrostatic pressure at each interface [Pa]
+  real :: p(SZI_(G),SZJ_(G),SZK_(G)+1)    ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
   real :: dz_geo(SZI_(G),SZJ_(G),SZK_(G)) ! The change in geopotential height
-                                          ! across a layer [m2 s-2].
+                                          ! across a layer [L2 T-2 ~> m2 s-2].
   real :: dilate(SZI_(G))                 ! non-dimensional dilation factor
   real :: htot(SZI_(G))                   ! total thickness [H ~> m or kg m-2]
-  real :: I_gEarth
+  real :: I_gEarth          ! The inverse of the gravitational acceleration times the
+                            ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
   real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names.
   integer i, j, k, isv, iev, jsv, jev, nz, halo
 
@@ -67,7 +68,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
   Z_to_eta = 1.0 ; if (present(eta_to_m)) Z_to_eta = US%Z_to_m / eta_to_m
   H_to_eta = GV%H_to_Z * Z_to_eta
   H_to_rho_eta =  GV%H_to_RZ * Z_to_eta
-  I_gEarth = Z_to_eta /  (US%Z_to_m * GV%mks_g_Earth)
+  I_gEarth = Z_to_eta / GV%g_Earth
 
 !$OMP parallel default(shared) private(dilate,htot)
 !$OMP do
@@ -99,13 +100,14 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
         ! ### THIS SHOULD BE P_SURF, IF AVAILABLE.
         do i=isv,iev ; p(i,j,1) = 0.0 ; enddo
         do k=1,nz ; do i=isv,iev
-          p(i,j,K+1) = p(i,j,K) + GV%H_to_Pa*h(i,j,k)
+          p(i,j,K+1) = p(i,j,K) + GV%g_Earth*GV%H_to_RZ*h(i,j,k)
         enddo ; enddo
       enddo
 !$OMP do
       do k=1,nz
         call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), &
-                                 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
+                                 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, &
+                                 SV_scale=US%R_to_kg_m3, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
       enddo
 !$OMP do
       do j=jsv,jev
@@ -159,11 +161,12 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
              !! the units of eta to m; by default this is US%Z_to_m.
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
-    p     ! The pressure at interfaces [Pa].
+    p          ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
-    dz_geo     ! The change in geopotential height across a layer [m2 s-2].
+    dz_geo     ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2].
   real :: htot(SZI_(G))  ! The sum of all layers' thicknesses [H ~> m or kg m-2].
-  real :: I_gEarth
+  real :: I_gEarth          ! The inverse of the gravitational acceleration times the
+                            ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
   real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names.
   integer i, j, k, is, ie, js, je, nz, halo
 
@@ -174,7 +177,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
   Z_to_eta = 1.0 ; if (present(eta_to_m)) Z_to_eta = US%Z_to_m / eta_to_m
   H_to_eta = GV%H_to_Z * Z_to_eta
   H_to_rho_eta =  GV%H_to_RZ * Z_to_eta
-  I_gEarth = Z_to_eta / (US%Z_to_m * GV%mks_g_Earth)
+  I_gEarth = Z_to_eta / GV%g_Earth
 
 !$OMP parallel default(shared) private(htot)
 !$OMP do
@@ -196,16 +199,18 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
     if (associated(tv%eqn_of_state)) then
 !$OMP do
       do j=js,je
+        ! ### THIS SHOULD BE P_SURF, IF AVAILABLE.
         do i=is,ie ; p(i,j,1) = 0.0 ; enddo
 
         do k=1,nz ; do i=is,ie
-          p(i,j,k+1) = p(i,j,k) + GV%H_to_Pa*h(i,j,k)
+          p(i,j,k+1) = p(i,j,k) + GV%g_Earth*GV%H_to_RZ*h(i,j,k)
         enddo ; enddo
       enddo
 !$OMP do
       do k = 1, nz
         call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, &
-                                 G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
+                                 G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, &
+                                 SV_scale=US%R_to_kg_m3, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
       enddo
 !$OMP do
       do j=js,je ; do k=1,nz ; do i=is,ie