*+Non-Boussinesq revision of set_viscosity

  This commit revises set_viscous_BBL and set_viscous_ML to work in fully
non-Boussinesq mode, eliminating all dependencies on the Boussinesq reference
density when in non-Boussinesq mode.  These changes include using the specific
volume interpolated to velocity points for the conversion between thickness and
height changes when tv%SpV_avg is allocated, including in the calculation of the
rescaled ustar in the top and bottom boundary layers and in the Rayleigh drag
calculation.  In set_viscous_ML, the derivatives of specific volume are used to
determine the reduced gravity across interfaces when in non-Boussinesq mode.

  Both set_viscous_BBL and set_viscous_ML now work extensively in height units
rather than thickness units in various places where it is more appropriate to do
so when in non-Boussinesq mode, and both use calls to thickness_to_dz to convert
between thicknesses and vertical distances.  In some places this leads to extra
calculations using separate arrays that are rescaled duplicates of each other in
Boussinesq mode, which will probably slow the model down a little.

  There is one rescaling bug due to a hard-coded unrescaled dimensional constant
that was fixed.  It occurs when setting a tiny floor in a thickness in one
inverse calculation.  However, this appears to be so small that it does not
change any answers in the MOM6-examples test suite, and could only do so if
ANGSTROM is set to be many orders of magnitude less than the typical value of
1e-10 m, so it was decided that no runtime bug-flag is needed in this case.

  A call to find_star is now used to specify the friction velocity used in
set_viscous_ML.  When in non-Boussinesq mode, this has the effect of using
forces%tau_mag and tv%SpV_avg instead of forces%ustar and GV%Rho0 to determine
the friction velocity.  In set_viscous_BBL the layer specific volumes are used
in non-Boussinesq mode to find the average velocities used in the linearization
of the nonlinear bottom drag, with similar changes in set_viscous_ML for the
drag under ice-shelves.

  When in non-Boussinesq mode, the units of Kd_shear, Kv_shear and Kv_shear_Bu
in the restart files are scaled to their natural MKS units of [Pa s] or [kg m-1
s-1] rather than [m2 s-1] to avoid round-off level changes across restarts, with
the units in the files documenting these changes.  These altered units do not
apply to Boussinesq mode restart files.

  In set_visc_init, the rescaling factor used to set CS%Hbbl was changed from
GV%Z_to_H to (US%Z_to_m*GV%m_to_H) so that it does not depend directly on RHO_0.

  These changes are extensive but localized to this file, with one new element
of the opaque set_visc_CS type and 6 elements with altered units, 37 new or
renamed internal variables, and 29 existing internal variables that have revised
units.  In particular one extensively used thickness curvature variable that had
been called just "a" was renamed "crv" to more accurately reflect its purpose
and to make it easier to find when rescaling it, with similar changes to two
other closely related variables.  These changes include the elimination of 29
rescaling factors between thickness and height units, as were two spots that
inappropriately (and unnecessarily, as it turns out) used GV%Rho0 even when in
non-Boussinesq mode.  There are also 4 new debugging checksum output calls.

  In Boussinesq or semiBoussinesq mode, no public interfaces are changed and all
answers are bitwise identical, but in non-Boussinesq mode answers will change
and there are changes to the units of 3 variables in restart files.