+Optionally use SSH in calculate density for PGF

[Hallberg-NOAA]

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.

[Hallberg-NOAA]

I have re-examined the question of whether the reference height for calculating the pressure used for the equation of state calls is correctly placed relative to the addition of self-attraction and loading adjustments to the geoid, and I am convinced that the order of these two parts of the code is correct.

[MJHarrison-GFDL]

I have reviewed these changes and agree. The new parameter flag SSH_IN_EOS_PRESSURE_FOR_PGF will correctly account for sea surface displacement in the calculation of in situ pressure gradient forces.

[herrwang0]

MOM6/src/core/MOM_density_integrals.F90

Line 250 in 99b95df

hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))

I think z_t also needs to be referenced to z0pres as well, and similarly for thickness weighting in other schemes.

[herrwang0]

For Boussinesq, if we incorporate tidal and SAL potential into interface height, there will be a baroclinic component of these forcing, which means for each layer tidal and SAL (both body forces ) term is different. Maybe we can argue that this is consistent with Boussinesq approximation, but I am not sure.

[Hallberg-NOAA]

MOM6/src/core/MOM_density_integrals.F90

Line 250 in 99b95df

hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))

I think z_t also needs to be referenced to z0pres as well, and similarly for thickness weighting in other schemes.

After re-inspecting the code, I am pretty sure that z_t does not need to be referenced z0pres because it is being taken relative to bathyT. However, I do think that there is a problem here in that z_t is adjusted for tidal and SAL contributions, but bathyT is not. This can be addressed by replacing the (optional) G%bathyT arguments to the various int_density_dz routines with e(:,:,nz+1), but this could change answers in some cases with tides or SAL enabled, so we should handle this with a bug-fix flag to recover the previous answers, and it should be dealt with separately from this PR.

[Hallberg-NOAA]

For Boussinesq, if we incorporate tidal and SAL potential into interface height, there will be a baroclinic component of these forcing, which means for each layer tidal and SAL (both body forces ) term is different. Maybe we can argue that this is consistent with Boussinesq approximation, but I am not sure.

This is a good question, and I think that the answer is ambiguous to when we are making the Boussinesq approximation. However, I also think that this question goes beyond the scope of what is being addressed with this PR.

[marshallward]

OBO @MJHarrison-GFDL

[marshallward]

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/pipelines/24788 ✔️ 🟡