(*)Refactor p_ave calculation

  Rearranged the calculation of the layer average pressure used in the density
integrals to have one less multiply and be clearer to read.  This could change
answers at roundoff if REFERENCE_HEIGHT is nonzero (it is 0 by default and in
all known cases) or if SSH_IN_EOS_PRESSURE_FOR_PGF is set to true.  Because the
former is mostly used for approximate self-consistency testing and the latter
option has just been added, it is unlikely that answers will change for any
production runs.

Author: Hallberg-NOAA

src/core/MOM_density_integrals.F90

@@ -268,8 +268,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
         pos = i*15+(m-2)*5
         T15(pos+1) = (wtT_L*T(i,j)) + (wtT_R*T(i+1,j))
         S15(pos+1) = (wtT_L*S(i,j)) + (wtT_R*S(i+1,j))
-        p15(pos+1) = -GxRho * ( ((wt_L*z_t(i,j)) + (wt_R*z_t(i+1,j))) - &
-                                ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j))) )
+        p15(pos+1) = -GxRho * ((wt_L*(z_t(i,j)-z0pres(i,j))) + (wt_R*(z_t(i+1,j)-z0pres(i+1,j))))
         do n=2,5
           T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
           p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i)
@@ -335,8 +334,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
         pos = i*15+(m-2)*5
         T15(pos+1) = (wtT_L*T(i,j)) + (wtT_R*T(i,j+1))
         S15(pos+1) = (wtT_L*S(i,j)) + (wtT_R*S(i,j+1))
-        p15(pos+1) = -GxRho * ( ((wt_L*z_t(i,j)) + (wt_R*z_t(i,j+1))) - &
-                                ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1))) )
+        p15(pos+1) = -GxRho * ((wt_L*(z_t(i,j)-z0pres(i,j))) + (wt_R*(z_t(i,j+1)-z0pres(i,j+1))))
         do n=2,5
           T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
           p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i)
@@ -627,8 +625,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
         S15(pos+1) = (w_left*Stl) + (w_right*Str)
         S15(pos+5) = (w_left*Sbl) + (w_right*Sbr)
 
-        p15(pos+1) = -GxRho * ( ((w_left*e(i,j,K)) + (w_right*e(i+1,j,K))) - &
-                                ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) )
+        p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i+1,j,K)-z0pres(i+1,j))))
 
         ! Pressure
         do n=2,5
@@ -724,8 +721,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
         S15(pos+1) = (w_left*Stl) + (w_right*Str)
         S15(pos+5) = (w_left*Sbl) + (w_right*Sbr)
 
-        p15(pos+1) = -GxRho * ( ((w_left*e(i,j,K)) + (w_right*e(i,j+1,K))) - &
-                                ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) )
+        p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i,j+1,K)-z0pres(i,j+1))))
 
         ! Pressure
         do n=2,5
@@ -1037,8 +1033,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
         dz_x(m,i) = (w_left*(e(i,j,K) - e(i,j,K+1))) + (w_right*(e(i+1,j,K) - e(i+1,j,K+1)))
 
         pos = i*15+(m-2)*5
-        p15(pos+1) = -GxRho * ( ((w_left*e(i,j,K)) + (w_right*e(i+1,j,K))) - &
-                                ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) )
+        p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i+1,j,K)-z0pres(i+1,j))))
         do n=2,5
           p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i)
         enddo
@@ -1143,8 +1138,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
         dz_y(m,i) = (w_left*(e(i,j,K) - e(i,j,K+1))) + (w_right*(e(i,j+1,K) - e(i,j+1,K+1)))
 
         pos = i*15+(m-2)*5
-        p15(pos+1) = -GxRho * ( ((w_left*e(i,j,K)) + (w_right*e(i,j+1,K))) - &
-                                ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) )
+        p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i,j+1,K)-z0pres(i,j+1))))
         do n=2,5
           p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i)
         enddo

src/equation_of_state/MOM_EOS_Wright.F90

@@ -592,8 +592,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j))))
+      p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + &
+                      (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
@@ -634,8 +634,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1))))
+      p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + &
+                      (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)

src/equation_of_state/MOM_EOS_Wright_full.F90

@@ -597,8 +597,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j))))
+      p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + &
+                      (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -639,8 +639,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1))))
+      p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + &
+                      (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)

src/equation_of_state/MOM_EOS_Wright_red.F90

@@ -599,8 +599,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j))))
+      p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + &
+                      (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -641,9 +641,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))
 
       dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
-      p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - &
-                      ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1))))
-
+      p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + &
+                      (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1))))
 
       I_al0 = 1.0 / al0
       I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)