*+Add BT_RHO_LINEARIZED to MOM_barotropic.F90

  Added the new runtime parameter BT_RHO_LINEARIZED to specify the density that
is used to convert total water column thicknesses into mass in non-Boussinesq
mode with linearized options in the barotropic solver or when estimating the
stable barotropic timestep without access to the full baroclinic model state.
The default is set to RHO_0 and answers do not change by default.  This new
parameter is used in non-Boussinesq mode with some options in btcalc and
find_face_areas, when LINEARIZED_BT_CORIOLIS = True or BT_NONLIN_STRESS = False,
and in the unit conversion of the ice strength with dynamic pressure.

  Also cancelled out factors of GV%Z_to_H in MOM_barotropic.F90 to simplify the
code and reduce the dependence on the value of GV%Rho_0 in non-Boussinesq mode.
This involved changing the units of 4 variables in the barotropic_CS type, 3
internal variables in btstep and an internal variable in barotropic_init to use
thickness units.  The rescaled internal variable mass_to_Z was also replaced
with the equivalent GV%RZ_to_H.  There are also 4 new debugging messages.  Also
modified the units of the gtot_est argument to match those of pbce.  There is a
new element in barotropic_CS.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode.  Additionally there is a new runtime parameter that will
appear in some MOM_parameter_doc files.