*Non-Boussinesq revision of MOM_thickness_diffuse

  Refactored thickness_diffuse when in non-Boussinesq mode to avoid any
dependencies on the Boussinesq reference density, and to translate the volume
streamfunction into the mass streamfunction using an appropriately defined
in-situ density averaged to the interfaces at velocity points.  This form
follows the suggestions of Appendix A.3.2 of Griffies and Greatbatch (Ocean
Modelling, 2012) when in non-Boussinesq mode.  Thickness_diffuse_full was also
revised to work properly in non-Boussinesq mode (and not depend on the
Boussinesq reference density) when no equation of state is used.

  As a part of these changes, the code now uses thickness-based streamfunctions
and other thickness-based internal calculations in MOM_thickness_diffuse.  For
example, the overturning streamfunctions with this change are now in m3/s in
Boussinesq mode, but kg/s in non-Boussinesq mode.  These changes use a call to
thickness_to_dz to set up a separate variable with the vertical distance across
layers, and in non-Boussinesq mode they use tv%SpV_avg to estimate in situ
densities.  Additional debugging checksums were added to thickness_diffuse.

  The code changes are extensive with 15 new or renamed internal variables, and
changes to the units of 9 other internal variables and 3 arguments to the
private routine streamfn_solver.  After this change, GV%Rho, GV%Z_to_H and
GV%H_to_Z are no longer used in any non-Boussinesq calculations (12 such
instances having been elimated).  Because some calculations have to be redone
with the separate thickness and dz variables, this will be more expensive than
the original version.

  No public interfaces are changed, and all answers are bitwise identical in
Boussinesq or semiBoussinesq mode, but they will change in non-Boussinesq mode
when the isopycnal height diffusion parameterization is used.