*Non-Boussinesq revision of set_diffusivity

  This commit revises the internal routines called by set_diffusivity to work in
fully non-Boussinesq mode, eliminating all dependencies on the Boussinesq
reference density when in non-Boussinesq mode.  The publicly visible interfaces
to this module and the external routines it calls have already been revised, so
only this file needs to be updated.  The specific changes include:

 - Refactored add_LOTW_BBL_diffusivity, add_MLrad_diffusivity, and set_BBL_TKE
   to work in units of layer vertical extent rather than layer thickness to give
   results in non-Boussinesq mode that avoid dependence on the Boussinesq
   reference density.

 - Work with internal variables in vertical distances in the denominator of
   diffusive flux calculations in find_TKE_to_Kd, while other expressions in
   this routine are recast in terms of thicknesses to avoid conversions.

 - Use tv%SpV_avg instead of 1/RHO_0 in find_TKE_to_Kd when in fully
   non-Boussinesq mode.

 - Use layer target density differences in place of g_prime in
   set_density_ratios in mathematically equivalent expressions when
   non-Boussinesq.

 - Use thickness_to_dz in 3 places to convert layer thicknesses into vertical
   distances.

 - The thickness argument to add_MLrad_diffusivity (in [H ~> m or kg m-2]) has
   been replaced with a vertical extent argument (in [Z ~> m]).

 - Use fluxes%tau_mag in place of fluxes%ustar in add_MLrad_diffusivity when in
   non-Boussinesq or semi-Boussinesq mode.  There is a new thermo_var_ptrs type
   argument to the internal routine add_MLrad_diffusivity to permit this
   changes.

 - Use the in situ near-bottom density when adding certain contributions to
   non-Boussinesq diffusivities.  This change includes the addition of a new
   bottom density (rho_bot) argument to add_int_tide_diffusivity,
   add_LOTW_BBL_diffusivity and add_drag_diffusivity.

 - Use GV%dZ_subroundoff in 4 spots in place of GV%H_to_Z*GV%H_subroundoff.

 - A long-standing comment questioning whether there is double-counting of tidal
   mixing has been addressed (there is not) and the comment has been revised
   accordingly.

  These changes involved changing the units of 21 internal variables and 1
element in the set_diffusivity_CS type.  There are 11 new internal variables,
while 9 internal variables were eliminated.  A total of 44 thickness rescaling
factors were eliminated, and there are 2 places where GV%Rho_0 was being used
explicitly that were changed into equivalent rescaling factors.

  All answers are bitwise identical in Boussinesq mode, but some solutions will
change in non-Boussinesq mode with this commit.