*Non-Boussinesq revision of diabatic_driver

  This commit completes the non-Boussinesq revision of the diabatic_driver code.
The revised code uses thickness_to_dz and works with internal variables in
vertical distances in the denominator of diffusive flux calculations in
diabatic_ALE_legacy, diabatic_ALE and layered_diabatic.  The code now uses
find_ustar to extract the friction velocities passed to KPP_compute.  The (tiny)
boundary layer tracer diffusivity is also rescaled to [H2 T-1].  With this set
of changes, all implicit references to Boussinesq reference density are
eliminated from the calculations in diabatic_driver when in non-Boussinesq mode.

  A total of 14 thickness rescaling factors were cancelled out, and there are 15
new or renamed variables in the diabatic routines.  4 new checksum calls are
also added to diabatic_ALE_legacy to assist in debugging.

  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 the use of the layer
averaged specific volume to convert thicknesses to vertical distances leads to
changing answers.