+(*)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.

  Within ePBL_column, the hp_a variable was converted from a scalar into a 1-d
array to emphasize where these variables are located and to bring ePBL_column
closer to what would be required if the new diffusivities were being applied on
top of diffusivities from other processes, as is the case in ePBL_BBL_column.
Several comments in find_PE_change and ePBL_column were corrected.  Two
allocatable three-dimensional diagnostic arrays that had been in the
energetic_PBL control structure were converted into local arrays to reduce the
dedicated memory requirements when these diagnostics are enabled.

  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)