Add MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP option for the top boundary, using
top interface height, analogously to what is done near the bathymetry when
MASS_WEIGHT_IN_PRESSURE_GRADIENT.  The information from the is parameter is
encoded in the MassWghtInterp value passed to the various int_density_... and
int_spec_vol_... routines, so the routine interfaces do not change.  For now
this new bit of information about whether to do mass weighting near the surface
under ice shelves is only used in int_density_dz_generic_plm.  By default this
new option is set to false and no answers are changed.

  Implemented the MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP related capabilities to
the 13 density or specific volume integral routines that did not have them yet,
mirroring what was already done in int_density_dz_generic_plm.  This includes
both the density integral routines that are used primarily for Boussinesq cases
and the specific volume integral routines that are primarily used with
non-Boussinesq simulations.

  This change includes the addition of an optional SSH argument to
int_density_dz and 6 other related routines (int_density_dz_generic_pcm,
analytic_int_density_dz, int_density_dz_wright, int_density_dz_wright_full,
int_density_dz_wright_red and int_density_dz_linear).  It also includes the
addition of an optional P_surf argument to int_specific_vol_dp and 7 other
related routines (int_spec_vol_dp_generic_pcm, int_spec_vol_dp_generic_plm,
analytic_int_specific_vol_dp int_spec_vol_dp_wright,
int_spec_vol_dp_wright_full, int_spec_vol_dp_wright_red and
int_spec_vol_dp_linear).  Both of these new optional arguments are required
when the new near-surface mass weighting is activated, and there are test that
would issue a fatal error if they are not provided in such cases.  (Note that
these routines can be called from other places than the pressure gradient force
calculation, in which case the calculations that need these fields are not
done.)

  The diagnose_mass_weight_Z and diagnose_mass_weight_p diagnostic routines were
similarly revised to mirror the expanded range of capabilities in the integral
routines.

  This commit can change answers when MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP is
true, but all answers are bitwise identical otherwise.  There are new arguments
to 15 publicly visible routines.

  Add CORRECTION_INTXPA which makes a quadratic correction to surface pressure
integrals (intx_pa and inty_pa) under ice based on the horizontal gradients of
the in-situ density  anomaly along the surface.  The non-Boussinesq version
corrects the topmost value of intx_za and inty_za and uses in-situ specific
volume gradients along the ocean surface.  By default new the runtime parameter
CORRECTION_INTXPA is false, and answers are unchanged.

  Add CORRECTION_INTXPA_5PT which uses 5 point quadrature to calculate surface
pressure integral and therefore could work with a nonlinear EOS, or (if added)
different subgrid distributions.  This option requires CORRECTION_INTXPA = True.
By default CORRECTION_INTXPA_5PT is false and no answers are changed.

  Add RESET_INTXPA_INTEGRAL which resets intxpa and intypa at a trusted cell in
the interior (non-vanished and non-MWIPG-affected), and then integrates both up
and down to update intxpa and intypa for the interfaces above and below.  This
also adds the new runtime parameter RESET_INTXPA_H_NONVANISHED and to determine
when a cell is trusted.  This option is recommended with
MASS_WEIGHT_IN_PRESSURE_GRADIENT_IS for a quiet zstar ice shelf.  By default,
this option is not on and no answers change, but there are new parameters in
some MOM_parameter_doc files.

  Refactored and revised the recently added code in MOM_PressureForce_FV.F90 to
reduce the number of calls to the equation of state routines, and corrected a
number of minor bugs in the original implementation.  Answers are bitwise
identical unless the new options to reset the pressure gradient calculations are
actively selected.

  Refactored the CORRECTION_INTX_PA calculations to avoid multiple intermediate
steps, while also adding comments documenting the derivation of the final
expression.  Also calculate the pressures used in the equation of state calls
with the Boussinesq CORRECTION_INTX_PA and RESET_INTXPA_INTEGRAL options
consistently with the other terms.  These changes will change the answers when
either of those options are in use, but are bitwise identical when they are not.

  Added non-Boussinesq versions of the code that is exercised when
CORRECTION_INTXPA_5PT and RESET_INTXPA_INTEGRAL are set to true.  These options
use the same names as in their Boussinesq forms, even though the arrays that are
being adjusted are actually intx_za and inty_za.  As a part of the testing of
this commit, several checksums were added to the PressureForce_FV routines;
these are enabled when DEBUG = True.  The changes in this commit will change the
answers in non-Boussinesq cases when either of those options are in use, but are
bitwise identical when they are not, and all Boussinesq answers are bitwise
identical.

  Added code to do 5-point quadrature integrals when RESET_INTXPA_INTEGRAL and
CORRECTION_INTXPA_5PT are true.  Also extended the ranged of interfaces that
can be used for the corrections to include the bottom interface and use the
surface interface for the nonlinear pressure gradient force corrections either
when the ocean surface interface is within the pressure or height range of the
top cell, or when no appropriate interior interface has been found (for lack of
a better idea).  This latter case should probably be reconsidered later.  This
commit will change answers if RESET_INTXPA_INTEGRAL is true and one description
of a runtime parameter in the MOM_parameter_doc files has been revised.

  Eliminated the runtime parameter CORRECTION_INTXPA_5PT by always selecting the
code that would have been used if it were true.  This decision was based on the
relative performance of the options with this set to true or false, and on the
desire to maintain consistency with the 5-point Boole's rule quadrature used
elsewhere in the pressure gradient force calculations.  CORRECTION_INTXPA_5PT
was only added in the same PR as this commit, so it can been simply removed and
there is no need to obsolete it.  This will change answers if CORRECTION_INTXPA
or RESET_INTXPA_INTEGRAL are true and CORRECTION_INTXPA_5PT was set to false,
but otherwise answers are bitwise identical.  There are fewer entries in the
MOM_parameter_doc files as a result of this commit.

  Set MASS_WEIGHT_IN_PRESSURE_GRADIENT = True in the MOM_input files for the
tc1, tc2 and tc4 test cases in .testing, to permit the code to pass this testing
while dealing with inconsistent settings for two diagnostics related to the
mass-weighting in the pressure gradient forces.  This commit redefines these 3
test cases.