+Optionally use SSH in calculate density for PGF

  Added the option of including the atmospheric or ice pressure and sea surface
height displacements from the global reference height in the pressures used in
the density calculations for Boussinesq pressure gradient calculations.  Note
that the full pressures were already being used everywhere apart from the
calculation of the equation of state.  This capability is controlled by the new
runtime parameter SSH_IN_EOS_PRESSURE_FOR_PGF.  This commit changes the Z_0p
argument to int_density_dz and 8 other routines (int_density_dz_generic_pcm,
int_density_dz_generic_plm, int_density_dz_generic_ppm, analytic_int_density_dz,
int_density_dz_wright, int_density_dz_wright_full, int_density_dz_wright_red
and int_density_dz_linear) from scalars into 2-d arrays, as were the internal
z0pres arrays in most of these routines.  By default, all answers are bitwise
identical, but there is a new runtime parameter in the MOM_parameter_doc files
for Boussinesq cases.

src/core/MOM_PressureForce_FV.F90

@@ -63,6 +63,9 @@ module MOM_PressureForce_FV
                             !! for the finite volume pressure gradient calculation.
                             !! By the default (1) is for a piecewise linear method
 
+  logical :: use_SSH_in_Z0p !< If true, adjust the height at which the pressure used in the
+                            !! equation of state is 0 to account for the displacement of the sea
+                            !! surface including adjustments for atmospheric or sea-ice pressure.
   logical :: use_stanley_pgf  !< If true, turn on Stanley parameterization in the PGF
   integer :: tides_answer_date !< Recover old answers with tides in Boussinesq mode
   integer :: id_e_tide = -1 !< Diagnostic identifier
@@ -522,6 +525,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
                   ! [Z ~> m].
     e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
                   ! specific to tides [Z ~> m].
+    Z_0p, &       ! The height at which the pressure used in the equation of state is 0 [Z ~> m]
     SSH, &      ! The sea surface height anomaly, in depth units [Z ~> m].
     dM          ! The barotropic adjustment to the Montgomery potential to
                 ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
@@ -577,6 +581,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
                              ! in roundoff and can be neglected [H ~> m].
   real :: I_Rho0             ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
   real :: G_Rho0             ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
+  real :: I_g_rho            ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]
   real :: rho_ref            ! The reference density [R ~> kg m-3].
   real :: dz_neglect         ! A minimal thickness [Z ~> m], like e.
   real :: H_to_RL2_T2        ! A factor to convert from thickness units (H) to pressure
@@ -753,6 +758,21 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
     enddo ; enddo
   endif
 
+  if (CS%use_SSH_in_Z0p .and. use_p_atm) then
+    I_g_rho = 1.0 / (CS%rho0*GV%g_Earth)
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho
+    enddo ; enddo
+  elseif (CS%use_SSH_in_Z0p) then
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      Z_0p(i,j) = e(i,j,1)
+    enddo ; enddo
+  else
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      Z_0p(i,j) = G%Z_ref
+    enddo ; enddo
+  endif
+
   do k=1,nz
     ! Calculate 4 integrals through the layer that are required in the
     ! subsequent calculation.
@@ -769,19 +789,19 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
                     G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
                     intx_dpa(:,:,k), inty_dpa(:,:,k), &
                     MassWghtInterp=CS%MassWghtInterp, &
-                    use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
+                    use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=Z_0p)
         elseif ( CS%Recon_Scheme == 2 ) then
           call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
                     rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
                     G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
                     intx_dpa(:,:,k), inty_dpa(:,:,k), &
-                    MassWghtInterp=CS%MassWghtInterp, Z_0p=G%Z_ref)
+                    MassWghtInterp=CS%MassWghtInterp, Z_0p=Z_0p)
         endif
       else
         call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
                   rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), &
                   intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, dz_neglect, &
-                  CS%MassWghtInterp, Z_0p=G%Z_ref)
+                  CS%MassWghtInterp, Z_0p=Z_0p)
       endif
       if (GV%Z_to_H /= 1.0) then
         !$OMP parallel do default(shared)
@@ -1042,6 +1062,12 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
   endif
   call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
                  "If true, calculate self-attraction and loading.", default=CS%tides)
+  call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", CS%use_SSH_in_Z0p, &
+                 "If true, include contributions from the sea surface height in the height-based "//&
+                 "pressure used in the equation of state calculations for the Boussinesq pressure "//&
+                 "gradient forces, including adjustments for atmospheric or sea-ice pressure.", &
+                 default=.false., do_not_log=.not.GV%Boussinesq)
+
   call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
                  "If True, use the ALE algorithm (regridding/remapping). "//&
                  "If False, use the layered isopycnal algorithm.", default=.false. )

src/core/MOM_density_integrals.F90

@@ -80,7 +80,8 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa,
   real,       optional, intent(in)  :: dz_neglect !< A minuscule thickness change [Z ~> m]
   integer,    optional, intent(in)  :: MassWghtInterp !< A flag indicating whether and how to use
                                            !! mass weighting to interpolate T/S in integrals
-  real,       optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
 
   if (EOS_quadrature(EOS)) then
     call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
@@ -139,7 +140,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
                                            !! mass weighting to interpolate T/S in integrals
   logical,    optional, intent(in)  :: use_inaccurate_form !< If true, uses an inaccurate form of
                                           !! density anomalies, as was used prior to March 2018.
-  real,       optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(in)  :: Z_0p  !< The height at which the pressure is 0 [Z ~> m]
 
   ! Local variables
   real :: T5((5*HI%iscB+1):(5*(HI%iecB+2)))  ! Temperatures along a line of subgrid locations [C ~> degC]
@@ -157,7 +159,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
   real :: dz         ! The layer thickness [Z ~> m]
   real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]
   real :: dz_y(5,HI%isc:HI%iec)   ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]
-  real :: z0pres     ! The height at which the pressure is zero [Z ~> m]
+  real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
   real :: hWght      ! A pressure-thickness below topography [Z ~> m]
   real :: hL, hR     ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
   real :: iDenom     ! The inverse of the denominator in the weights [Z-2 ~> m-2]
@@ -184,7 +186,13 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
 
   GxRho = G_e * rho_0
   I_Rho = 1.0 / rho_0
-  z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
+  if (present(Z_0p)) then
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      z0pres(i,j) = Z_0p(i,j)
+    enddo ; enddo
+  else
+    z0pres(:,:) = 0.0
+  endif
   use_rho_ref = .true.
   if (present(use_inaccurate_form)) then
     if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form
@@ -209,7 +217,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
       dz = z_t(i,j) - z_b(i,j)
       do n=1,5
         T5(i*5+n) = T(i,j) ; S5(i*5+n) = S(i,j)
-        p5(i*5+n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz)
+        p5(i*5+n) = -GxRho*((z_t(i,j) - z0pres(i,j)) - 0.25*real(n-1)*dz)
       enddo
     enddo
 
@@ -260,7 +268,8 @@ 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))) - z0pres)
+        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))) )
         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)
@@ -326,7 +335,8 @@ 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))) - z0pres)
+        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))) )
         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)
@@ -414,7 +424,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
                                            !! mass weighting to interpolate T/S in integrals
   logical,    optional, intent(in)  :: use_inaccurate_form !< If true, uses an inaccurate form of
                                            !! density anomalies, as was used prior to March 2018.
-  real,       optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
 
 ! This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
@@ -464,7 +475,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
   real :: massWeightToggle          ! A non-dimensional toggle factor (0 or 1) [nondim]
   real :: Ttl, Tbl, Ttr, Tbr        ! Temperatures at the velocity cell corners [C ~> degC]
   real :: Stl, Sbl, Str, Sbr        ! Salinities at the velocity cell corners [S ~> ppt]
-  real :: z0pres                    ! The height at which the pressure is zero [Z ~> m]
+  real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
   real :: hWght                     ! A topographically limited thickness weight [Z ~> m]
   real :: hL, hR                    ! Thicknesses to the left and right [Z ~> m]
   real :: iDenom                    ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
@@ -480,7 +491,13 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
 
   GxRho = G_e * rho_0
   I_Rho = 1.0 / rho_0
-  z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
+  if (present(Z_0p)) then
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      z0pres(i,j) = Z_0p(i,j)
+    enddo ; enddo
+  else
+    z0pres(:,:) = 0.0
+  endif
   massWeightToggle = 0.
   if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif
   use_rho_ref = .true.
@@ -517,7 +534,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
     do i = Isq,Ieq+1
       dz(i) = e(i,j,K) - e(i,j,K+1)
       do n=1,5
-        p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz(i))
+        p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres(i,j)) - 0.25*real(n-1)*dz(i))
         ! Salinity and temperature points are linearly interpolated
         S5(i*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * S_b(i,j,k)
         T5(i*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * T_b(i,j,k)
@@ -610,7 +627,8 @@ 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))) - z0pres)
+        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))) )
 
         ! Pressure
         do n=2,5
@@ -706,7 +724,8 @@ 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))) - z0pres)
+        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))) )
 
         ! Pressure
         do n=2,5
@@ -812,7 +831,8 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
                                            !! divided by the y grid spacing [R L2 T-2 ~> Pa]
   integer,    optional, intent(in)  :: MassWghtInterp !< A flag indicating whether and how to use
                                            !! mass weighting to interpolate T/S in integrals
-  real,       optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
 
 ! This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
@@ -864,7 +884,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
   real :: t6 ! PPM curvature coefficient for T [C ~> degC]
   real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T [C ~> degC]
   real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S [S ~> ppt]
-  real :: z0pres ! The height at which the pressure is zero [Z ~> m]
+  real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
   real :: hWght  ! A topographically limited thickness weight [Z ~> m]
   real :: hL, hR ! Thicknesses to the left and right [Z ~> m]
   real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
@@ -879,7 +899,13 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
 
   GxRho = G_e * rho_0
   I_Rho = 1.0 / rho_0
-  z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
+  if (present(Z_0p)) then
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      z0pres(i,j) = Z_0p(i,j)
+    enddo ; enddo
+  else
+    z0pres(:,:) = 0.0
+  endif
   massWeightToggle = 0.
   if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif
 
@@ -924,7 +950,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
       endif
       dz = e(i,j,K) - e(i,j,K+1)
       do n=1,5
-        p5(I*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz)
+        p5(I*5+n) = -GxRho*((e(i,j,K) - z0pres(i,j)) - 0.25*real(n-1)*dz)
         ! Salinity and temperature points are reconstructed with PPM
         S5(I*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) )
         T5(I*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) )
@@ -1011,7 +1037,8 @@ 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))) - z0pres)
+        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))) )
         do n=2,5
           p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i)
         enddo
@@ -1116,7 +1143,8 @@ 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))) - z0pres)
+        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))) )
         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.F90

@@ -1345,7 +1345,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS,
   real,       optional, intent(in)  :: dz_neglect !< A miniscule thickness change [Z ~> m]
   integer,    optional, intent(in)  :: MassWghtInterp !< A flag indicating whether and how to use
                                           !! mass weighting to interpolate T/S in integrals
-  real,       optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(in)  :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
 
   ! Local variables
   real :: rho_scale  ! A multiplicative factor by which to scale density from kg m-3 to the