diag_send_complete() is removed from disable_averaging(), due to
potential performance costs over sweeping through files and variables
for synchronization across nodes or other thread-like backround
operations.

It also prevents potential errors in drivers where
diag_manager_set_time_end was unset.

We will come back to this issue and work out a better strategy for
diagnostic synchronization.

Three FMAs were introduced when the BBL convex L(:) computation were
moved to a separate function.  Parentheses have been added to disable
these FMAs.

Calculation of the cell width fractions `L` was consolidated and moved
into a new function, `find_L_open_concave_analytic` (later renamed to
`find_open_L_open_concave_trigonometric`).  When compiled with Intel
and using the following flags,

    -O2 -march=core-avx2 -fp-model source -qno-opt-dynamic-align

additional FMAs are added to the function, specifically the following
expressions:

    crv_3 = (Dp + Dm - 2.0*D_vel)

    L(K) = L0*(1.0 + ax2_3apb*L0)

    L(K) = apb_4a * (1.0 - 2.0 * cos(
        C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)
    )

(Third expression FMA was added inside of `acos()` by the compiler)

The crucial flag seems to be `-march=core-avx2` which controls the
introduction of FMAs and explains why the problem went undetected.

The legacy answers were restored by adding new parentheses to the above
equations.

Since this function is not generally considered reproducible and is
provided for legacy support, the loss of performance is not considered
to be an issue.

Five fields associated with KPP were previously initialized during
allocation (using `source=0.`).  It was found that the allocations were
unnecessary, and the fields were redefined as local variables within
functions, without the initialization to zero.

Although the uninitialized points were masked and unused in the rest of
the model, they were still able to change the checksums and were being
reported as regressions.  This patch restores the original checksums by
re-implementing the zero initialization.

Thanks to Bob Hallberg (@Hallberg-NOAA) for diagnosing this issue.

Due to answer changes in multiple production experiments, this patch
introduces a new parameter, USE_HUYNH_STENCIL_BUG, which sets the tracer
advection stencil width to its previous incorrect value.  This parameter
replicates the previous method, which keeps `stencil` at 2 when PPM:H3
is used.

----

The logical parameters are

* P = CS%usePPM
* H = CS%useHuynh
* E = USE_HUYNH_STENCIL_BUG
* R = CS%useHuynhStencilBug

The three expressions of interest:

1. P & ~H
2. P
3. P & ~R

(1) is the original incorrect expression, (2) is the fixed expression,
and (3) is the proposed update.

What we want:

If E is false, then (3) should reduce to (2).  If E is true, then (3)
should reduce to (1).

R is computed as follows:

* R = False  (from derived type initialization)
* if (H) R = E  (from get_param call)

This is equivalent to R = H & E, and (3) is P & ~(H & E).

* If E is False, then P & ~(H & False) = P.
* If E is True, then P & ~(H & True) = P & ~H

So this flag should replicate both the previous and current behavior.

  This commit restores the effectiveness of the runtime parameter
USE_WRIGHT_2ND_DERIV_BUG in determining whether a bug is corrected in the
calculation of two of the second derivative terms returned by
calculate_density_second_derivs_elem() with the "WRIGHT" equation of state,
recreating the behavior (and answers) that are currently on the main branch of
MOM6.  To do this, it adds and calls the new routine set_params_buggy_Wright()
when appropriate, and adds the new element "three" to the buggy_Wright_EOS type.
When the bug is fixed, buggy_Wright_EOS%three = 3, but ...%three = 2 to
recreate the bug.  This commit does change answers for cases using the "WRIGHT"
equation of state and one of the "USE_STANLEY_..." parameterizations from those
on the dev/gfdl branch of MOM6, but in so doing it restores the answers on the
main branch of MOM6.  There is also a new publicly visible subroutine.

  Moved a check for whether to initialize an uninitialized visc%h_ML field from
visc%MLD outside of the test for the value of CS%mixedlayer_restart.  This
change will prevent answer changes between code versions for certain cases that
initialize from a restart file, use the melt_potential field and have
MLE_USE_PBL_MLD = True but also have MIXEDLAYER_RESTRAT = False.  This commit
corrects a very specific oversight that was introduced when the roles of
visc%h_ML and visc%MLD were separated, and it can change answers (reverting them
to answers from older versions of the main branch of MOM6) in very specific
cases where MOM6 is initialized from a restart file from an older versions of
the code, including the 5x5-degree test case in the dev/emc regression suite.

  Added checksum calls for the melt_potential, ocean_mass, ocean_heat and
ocean_salt elements of the surface state in MOM_surface_chksum if these fields
are allocated for more comprehensive debugging.  Also added the symmetric
optional argument to the call to MOM_surface_chksum form extract_surface_state
so that all of the surface velocity values that could contribute to the ocean
surface velocities that are seen by the coupler are checksummed.  All solutions
are bitwise identical, but there are enhancements to the MOM6 debugging
capabilities.

GFDL to main (2024-05-31)

Resolved conflicts in src/core/MOM_checksum_packages.F90 with
collision between switch to unscale= and the addition of new
checksums with previous scale= argument.