This commit refactors entrainment_diffusive to avoid any dependencies on the
Boussinesq reference density when in non-Boussinesq mode, including using
calculate_specific_vol_derivs for one diagnostic when non-Boussinesq.

  This commit includes making the formerly optional arguments kb_out, Kd_Lay and
Kd_int to entrainment_diffusive non-optional, as they have been used in all
calls to this routine for many years.

  Layer target density differences are now used in place of g_prime in
entrainment_diffusive in mathematically equivalent expressions for the density
difference ratios when non-Boussinesq.

  The non-Boussinesq upper layer buoyancy flux with entrainment_diffusive was
revised to avoid using the Boussinesq reference density when the model is in
layered mode but there is no equation of state or bulk mixed layer in use.

  The units of 3 internal variables were changed and there are 3 new internal
variables as a part of these changes, and 4 thickness rescaling factors were
eliminated.  A default private setting is used to simplify a block of OMP
directives.

  All answers are bitwise identical in Boussinesq mode, but they can change for
some non-Boussinesq configurations, and three previously optional arguments have
been made mandatory.