+(*)Add ePBL bottom boundary mixing option

  Add the option to do energetically consistent bottom boundary layer mixing
with the new routine ePBL_BBL_column.  ePBL_BBL_column is closely based on the
surface-focused ePBL mixing in ePBL_column, but without adding convective
instability driven mixing or mean-TKE driven mixing to avoid possible
double-counting.  This new option is enabled by setting the new runtime
parameter EPBL_BBL_EFFIC to be positive.

  If both EPBL_BBL_EFFIC and BBL_EFFIC are set to positive values, there is a
risk of double-counting, but this case is not being trapped for now.

  The changes include the addition of a new mandatory vertvisc_type argument to
the publicly visible routine energetic_PBL.

  When this new ePBL bottom boundary layer mixing option is enabled, there are
several new diagnostics available that are related to bottom boundary layer
mixing.  Several new checksum calls were also added with this new option when
DEBUG = True. The MOM_parameter_doc files are altered by the addition of two new
runtime parameters, and by the correction of several spelling errors in the
descriptions of other ePBL parameters.  By default, all answers are bitwise
identical.

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -858,7 +858,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
-    call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
+    call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, &
                        CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
 
     call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US)
@@ -1410,7 +1410,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
-    call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
+    call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, &
                        CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
 
     call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US)