Bugfix for mixing up Rho0 and rho_ref in Boussinesq PGF

[herrwang0]

This PR fixes a bug in Boussinesq finite volume pressure gradient force that GV%Rho0, mean seawater density and CS%Rho0, reference density are incorrectly interchanged in a number of instances. The bug probably has zero effect on any existing runs, as I am unaware of cases where RHO_PGF_REF is set different from RHO_0, otherwise, this bug would have been detected. So the fix here is just for future-proof purpose.

• A new runtime parameter RHO_PGF_REF_BUG is added.
• CS%Rho0 is renamed to CS%rho_ref for clarification.
• A new local variable GxRho_ref is added to distinguish from GxRho0. GxRho_ref is used to offset reference density in interface heights in pa. And GxRho0 is used to calculate surface pressure. Previously GxRho0 was incorrectly using reference density rho_ref (i.e. CS%Rho0).
• Subroutine calls for density integral and baroclinic pressure force pbce were incorrectly using rho_ref instead of GV%rho0. Two local variables rho0_int_density and rho0_set_pbce are added to maintain and fix this issue.

This PR does not change answers.

054354f: Remove unused variable I_Rho in three int_density_dz subroutines
93a86dc: Add runtime flag RHO_PGF_REF_BUG

[Hallberg-NOAA]

Having carefully examined these proposed changes, I am convinced that they are correct. In particular, any choosing any value of RHO_PGF_REF should give mathematically equivalent answers, and these changes should recover this property.

However, in examining these changes, I think that I have spotted a separate bug related to a missing reference geopotential height (G%Z_ref) correction. Please take a look at the relevant line, @herrwang0, and let me know what you think. In particular, answers should be mathematically equivalent when different values of REFERENCE_HEIGHT are used. I leave it up to you, though, @herrwang0, to decide whether this additional bug (if it does exist) should be dealt with in this same PR or whether it should be dealt with separately. Note that in every case I could find, REFERENCE_HEIGHT = 0.0, so such a bug seems very unlikely to change answers.

[Hallberg-NOAA]

With the recent minor update, I think that everything in this PR looks correct. Thank you for this valuable contribution.

[Hallberg-NOAA]

This PR has passed pipeline testing at https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/pipelines/26508 with the expected warnings about a new runtime parameter.