Add entrainment rule BLD calculation to KPP

This commit introduces an alternative method for computing the boundary layer depth
(BLD) in KPP, based on the entrainment rule. This method is activated when the
non-solar surface buoyancy flux (surfBuoy_NS) is negative and STOKES_MOST = True.
In all other cases, the standard Richardson Number-based method is used.

The non-solar surface buoyancy flux is estimated internally using an exponential
attenuation function, since only the total surface buoyancy flux (solar + non-solar)
is currently available to the module. Ideally, the non-solar component should be
provided explicitly.

Several new diagnostics have been added, and some of the existing ones have been
updated. Specifically, the parameterized shear, buoyancy, and Stokes TKE production
terms are calculated. We conducted multiple tests where these terms are passed to
the MLE module and used in the denominator of the Bodner streamfunction.
However, because of the small Cr value we adopted (0.01), this had minimal effect
on the solutions, so we decided to keep the original simpler implementation
and not pursue this option for now.

This commit should be evaluated together with CVMix/CVMix-src#106.

Author: gustavo-marques