The autoconf Python interpreter search was slightly modified to search
for Python even if $PYTHON is set to an empty string.  This is done by
unsetting PYTHON if it is set but empty, then following the usual macro.

This was required since `export PYTHON` in a Makefile will create the
`PYTHON` variable but will assign it no value (i.e. empty string).
This causes issues in some build environments.

The backup `configure~` script was also added to the developer
`ac-clean` cleanup rule.

  Refactored 6 files in the ALE directory to calculate the nominal depth in
thickness units in a single place (it is done in regridding_main now) and pass
it to the various places where it is used, in a preparatory step to modify how
this calculation is done in non-Boussinesq mode.  There are new arguments to
several publicly visible routines, including:

- Add non_depth_H arguments to hybgen_regrid, build_zstar_grid,
  build_sigma_grid, build_rho_grid, build_grid_HyCOM1, build_grid_adaptive,
  build_adapt_column and build_grid_arbitrary

- Add optional zScale arguments to build_zstar_grid and build_grid_HyCOM1

- Add unit_scale_type arguments to regridding_main, ALE_regrid_accelerated and
  ALE_offline_inputs

  Also eliminated an incorrect rescaling GV%Z_to_H facto when calculating the
total column thickness from the layer thicknesses when an ice shelf is used
with a Hycom grid.  This would have caused dimensional consistency testing to
fail.

  Added the new runtime parameters HYBGEN_H_THIN, HYBGEN_FAR_FROM_SURFACE
HYBGEN_FAR_FROM_BOTTOM, and HYBGEN_DENSITY_EPSILON to set previously hard-coded
dimensional parameters used in the Hybgen regridding code and store these
values in new variables in hybgen_regrid_CS.  Two of these are no longer passed
to hybgen_column_regrid as separate parameters.  By default these new runtime
parameters recover the previous hard-coded values.

  Also eliminated an unused block of code in build_rho_column.  Several
comments documenting variables or their units were also added.

  All answers are bitwise identical, but there are 4 new runtime parameters
that would appear in some MOM_parameter_doc files and there are changes to the
arguments to 11 routines.

As described in issue mom-ocean#372, I would like to be able to create restart files that
contain information about the particle location. These files will be written at
the same time as other restart files. I cannot add these calls directly to the
driver, because the driver does not have information about the particle location.
We have added save_MOM6_internal_state as a subroutine in MOM.F90, and
we added calls to this subroutine from each of the drivers. We hope this will
allow for more new packages to write restart files in the future.

Co-authored by Spencer Jones <[email protected]>

  Added the integer valid_SpV_halo to the thermo_var_ptrs type to indicate
whether the SpV_array has been updated and its valid halo size, to facilitate
error detection and debugging in non-Boussinesq mode.  Tv%valid_SpV_halo is set
to the halo size in calc_derived_thermo or after a halo update is done to
tv%SpV_avg, and it is set to a negative value right after calls that change
temperatures and salinities (such as by ALE remapping) unless there is a call to
calc_derived_thermo.  Tests for the validity of tv%SpV_avg are added to the
routines behind thickness_to_dz, with fatal errors issued if invalid arrays
would be used, but more tests could perhaps be used in any parameterization
routines where tv%SpV_avg is used directly.  Handling the updates to tv%SpV_avg
this way helps to avoid unnecessary calls to calc_derived_thermo, which in turn
has equation of state calls that can be expensive, while also providing
essential verification of new code related to the non-Boussinesq code.  These
tests can probably be commented out or removed for efficiency once there is a
full suite of regression tests for the fully non-Boussinesq mode of MOM6.  In
addition, a new optional debug argument was added to calc_derived_thermo which
can be used to triggers checksums for the variables used to calculate
tv%SpV_avg.  One call to calc_derived_thermo was also added just before the
initialization call to ALE_regrid that will be needed with the next commit, but
does not change answers yet.  All answers are bitwise identical, but there is a
new element in a transparent and widely used type and a new optional argument to
a public interface.

  Use RHO_KV_CONVERT instead of RHO_0 to set the non-Boussinesq version of
GV%m_to_H, so that there is a mechanism for testing the independence of the
fully non-Boussinesq mode from the Boussinesq reference density.  With this
change, GV%Z_to_H is not guaranteed to be equal to (GV%Z_to_m*GV%m_to_H), with
the latter expression preferred when setting parameters. By default the two
parameters are the same, and they will probably only ever differ in testing the
code.  All Boussinesq solutions are bitwise identical, but there are differences
in the description of RHO_KV_CONVERT that will appear in MOM_parameter_doc
files.

  Add new arguments to 7 routines that will be needed for the non-Boussinesq
capability, but do not use them yet, so that there will be fewer cross file
dependencies as the various changes are being reviewed simultaneously.  The
impacted interfaces are MEKE_int, vertvisc_coef, sumSWoverBands, KPP_calculate,
differential_diffuse_T_S, set_BBL_TKE, and apply_sponge In the three
step_MOM_dyn_... routines and in calculateBuoyancyFlux1d, this change includes
calls to thickness_to_dz to calculate the new vertical distance arrays that will
be passed into vertvisc_coef or sumSWoverBands.  The only place where the new
arguments are actually used is in sumSWoverBands and set_opacity where the
changes are particularly simple.  All answers are bitwise identical, but there
are new non-optional arguments to seven publicly visible routines.

* Restore functionality for reading slices from 3d volumes in MOM_io

 - The recent MOM_io modifications in support of FMS2_io accidentally
   removed support for reading on-grid data (same horizontal grid as
   model) k-slices. This is needed in some configurations in the model
   state initialization.

* Add FMS1 interfaces

* Additional patches to enable reading ongrid state initialization data

 - read local 3d volume rather than attempting to slice ongrid
   data vertically.

 - Related bugfixes in MOM_io

- We were reading KV_ML_INVZ2 without logging, then checking for KVML
  and finally logging based on a combination of the two. This had the side
  affect that we get warnings about not using KVML even if KVML was not
  present.
- The fix checks for KVML first, and then changes the default so that
  when KVML=1e-4 is replaced by KV_ML_INVZ2=1e-4 we end up with no
  warnings and KVML can be obsoleted safely.
  Note: this commit alone does not remove all warnings from the MOM6-examples
  suite because we still need to fix the MOM_input that still use KVML
- KVML needs to be unscaled since it is the default for KV_ML_INVZ2
- tc3 used KVML and has been corrected.

  Use the new runtime parameter RHO_PGF_REF instead of RHO_0 to set the
reference density that is subtracted off from the other densities when
calculating the finite volume pressure gradient forces.  Although the answers
are mathematically equivalent for any value of this parameter, a judicious
choice can reduce the impacts of roundoff errors by about 2 orders of
magnitude.  By default, RHO_PGF_REF is set to RHO_0, and all answers are bitwise
identical.  However, there is a new runtime parameter that appears in many of
the MOM_parameter_doc.all files.

The message that a file is being created was issued as a WARNING
when we all agree it should really be a NOTE. Depth_list.nc is
read if it is present to avoid recomputing a sorted list, but the
absence of the file is not an error and does not warrant a warning.

Changes:
- Changed WARNING to NOTE.
- Removed MOM_mesg from imports since it wasn't being used.

The interpolation scheme for state-dependent diagnostic coordinates
was incorrectly registering as the same parameter as the main model.
This meant it was never possible to change the interpolation scheme
from the default (which was not the same as the main model).

Fix registers the generated parameter name which was always computed
but not used. A typical example of the generated parameter is
"DIAG_COORD_INTERP_SCHEME_RHO2".

  Fixed a bug in which wave_speed_init was effectively discarding any values of
mono_N2_depth passed to it via the optional argument mono_N2_depth, but also
changed the default value of RESOLN_N2_FILTER_DEPTH, which was previously being
discarded, to disable the monotonization and replicate the previous results.
There were also clarifying additions made to the description how to disable
RESOLN_N2_FILTER_DEPTH.  This will change some entries in MOM_parameter_doc
files, and it will change solutions in cases that set RESOLN_N2_FILTER_DEPTH to
a non-default value and have parameter settings that use the resolution function
to scale their horizontal mixing.  There are, however, no known active
simulations where the answers are expected to change.

  Revised the calculation of gprime and the coordinate densities (GV%Rlay) in
fully non-Boussinesq mode to use the arithmetic mean of adjacent coordinate
densities in the denominator of the expression for g_prime in place of RHO_0.
Also use LIGHTEST_DENSITY in place of RHO_0 to specify the top-level coordinate
density in certain coordinate modes.  Also made corresponding changes to the
fully non-Boussinesq APE calculation when CALCULATE_APE is true, and eliminated
an incorrect calculation of the layer volumes in non-Boussinesq mode using the
Boussinesq reference density that was never actually being used when
CALCULATE_APE is false.

  This commit will change answers in some fully non-Boussinesq calculations, and
an existing runtime parameter is used and logged in some new cases, changing
the MOM_parameter_doc file in those cases.

  Refactored thickness_diffuse when in non-Boussinesq mode to avoid any
dependencies on the Boussinesq reference density, and to translate the volume
streamfunction into the mass streamfunction using an appropriately defined
in-situ density averaged to the interfaces at velocity points.  This form
follows the suggestions of Appendix A.3.2 of Griffies and Greatbatch (Ocean
Modelling, 2012) when in non-Boussinesq mode.  Thickness_diffuse_full was also
revised to work properly in non-Boussinesq mode (and not depend on the
Boussinesq reference density) when no equation of state is used.

  As a part of these changes, the code now uses thickness-based streamfunctions
and other thickness-based internal calculations in MOM_thickness_diffuse.  For
example, the overturning streamfunctions with this change are now in m3/s in
Boussinesq mode, but kg/s in non-Boussinesq mode.  These changes use a call to
thickness_to_dz to set up a separate variable with the vertical distance across
layers, and in non-Boussinesq mode they use tv%SpV_avg to estimate in situ
densities.  Additional debugging checksums were added to thickness_diffuse.

  The code changes are extensive with 15 new or renamed internal variables, and
changes to the units of 9 other internal variables and 3 arguments to the
private routine streamfn_solver.  After this change, GV%Rho, GV%Z_to_H and
GV%H_to_Z are no longer used in any non-Boussinesq calculations (12 such
instances having been elimated).  Because some calculations have to be redone
with the separate thickness and dz variables, this will be more expensive than
the original version.

  No public interfaces are changed, and all answers are bitwise identical in
Boussinesq or semiBoussinesq mode, but they will change in non-Boussinesq mode
when the isopycnal height diffusion parameterization is used.

* Salt data structures

* First steps at brine plume: pass info from SIS2

* The brine plume parameterization,

- including now passing the dimensional scaling tests.

* Fix problem when running Tidal_bay case with gnu.

* Avoiding visc_rem issues inside land mask.

Tweaking the brine plume code.

* Using the proper MLD in the brine plumes

- it now works better on restart

* Always including MLD in call to applyBoundary...

- I could move it up and make it not optional.

* Adding some OpenMP directives to brine plumes

This commit brings the drifters interface up-to-date with the current version of the drifters package, which requires h (layer thickness) to calculate the vertical movement of particles. The interfaces in the code and in config_src/external are updated to pass this information to the drifters package.

  Pass dt_kappa_smooth to calc_isoneutral_slopes and vert_fill_TS in units of
[H Z ~> m2 or kg m-1] instead of [Z2 ~> m2] for consistency with the units of
other diffusivities in the code and to reduce the depenency on the  Boussinesq
reference density in non-Boussinesq configurations.  In addition to the changes
to the units of these two arguments, there is a new unit_scale_type argument to
vert_fill_TS and MOM_calc_varT and a new verticalGrid_type argument to
MOM_stoch_eos_init.  The units of 4 vertical diffusivities in the control
structures in 4 different modules are also changed accordingly.

  All answers are bitwise identical in Boussinesq mode, but they can change for
some non-Boussinesq configurations.  There are new mandatory arguments to three
publicly visible routines.

  Added a comment justifying the use of a fixed rescaling factor for the
diffusivity used in vert_fill_TS.  All answers and output are identical.

  Added the new public interface find_ustar to extract the friction velocity
from either a forcing type argument, or a mech_forcing_type argument, either
directly or from tau_mag, and in non-Boussinesq mode by using the time-evolving
surface specific volume.  Find_ustar is an overloaded interface to
find_ustar_fluxes or find_ustar_mech_forcing, which are the same but for the
type of one of their arguments.  For now, the subroutines bulkmixedlayer,
mixedlayer_restrajt_OM4, mixedlayer_restrat_Bodner and mixedlayer_restrat_BML
are calling find_ustar to avoid code duplication during the transition to work
in fully non-Boussinesq mode, but it will eventually be used in about another
half dozen other places.

  All Boussinesq answers are bitwise identical, but non-Boussinesq answers will
change and become less dependent on the Boussinesq reference density, and there
is a new publicly visible interface wrapping two subroutines.

  Changed the units of the optional mono_N2_depth argument to wave_speed,
wave_speed_init and wave_speed_set_param in thickness units instead of height
units.  Accordingly, the units of one element each in the diagnostics_CS and
wave_speed_CS and a local variable in VarMix_init are also changed to thickness
units.  The unit descriptions of some comments describing diagnostics were also
amended to also describe the non-Boussinesq versions.  Because this is
essentially just changing when the unit conversion occurs, all answers are
bitwise identical, but there are changes to the units of an optional argument in
3 publicly visible routines.

  Added the new runtime parameter BT_RHO_LINEARIZED to specify the density that
is used to convert total water column thicknesses into mass in non-Boussinesq
mode with linearized options in the barotropic solver or when estimating the
stable barotropic timestep without access to the full baroclinic model state.
The default is set to RHO_0 and answers do not change by default.  This new
parameter is used in non-Boussinesq mode with some options in btcalc and
find_face_areas, when LINEARIZED_BT_CORIOLIS = True or BT_NONLIN_STRESS = False,
and in the unit conversion of the ice strength with dynamic pressure.

  Also cancelled out factors of GV%Z_to_H in MOM_barotropic.F90 to simplify the
code and reduce the dependence on the value of GV%Rho_0 in non-Boussinesq mode.
This involved changing the units of 4 variables in the barotropic_CS type, 3
internal variables in btstep and an internal variable in barotropic_init to use
thickness units.  The rescaled internal variable mass_to_Z was also replaced
with the equivalent GV%RZ_to_H.  There are also 4 new debugging messages.  Also
modified the units of the gtot_est argument to match those of pbce.  There is a
new element in barotropic_CS.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode.  Additionally there is a new runtime parameter that will
appear in some MOM_parameter_doc files.

  Use thickness_to_dz to convert layer thicknesses to depths in 4 tracer modules
(DOME_tracer, dye_example, ideal_age_example and nw2_tracers) so that this
conversion is done correctly in non-Boussineq mode, and there is no longer any
dependency on the Boussinesq reference density in that mode.

  This change includes the addition of a thermo_var_ptrs argument to 5 routines
(initialize_DOME_tracer, initialize_dye_tracer, dye_tracer_column_physics
ideal_age_tracer_column_physics and count_BL_layers) and changes to the units of
some internal variables, and the addition of 6 new 2-d or 3-d arrays with the
vertical distance across layers.  An unused param_file_type argument to
initialize_DOME_tracer was also eliminated.

  Comments were also added to describe the units of 5 of the variables in the
ideal age tracer control structure and 7 internal variables in that same module,
and there was some minor cleanup of the formatting cf calls in
tracer_flow_control_init.

  There was some minor refactoring in the ns2_tracers module to use SZK_(GV)
instead of SZK_(G) to declare the vertical extent of some arrays, and the
vertical indexing convention for interfaces in nw2_tracer_dist was revised from
starting at 0 to start at 1 for consistency with all the other code in MOM6.

  Also moved the code to do halo updates for the physical model state variables
and call calc_derived_thermo before calling tracer_flow_control_init, because
some routines there are now using the layer average specific volume to convert
between thicknesses and heights when in non-Boussinesq mode.

  All answers in Boussinesq mode are bitwise identical, but these passive tracer
modules have slightly different answers in non-Boussinesq mode.  There are
changes to the non-optional arguments to 4 public interfaces.

  Changed a recently added OMP directive for plume_flux from private to
firstprivate to reflect how this variable is actually used.  This bug was
introduced with PR mom-ocean#401, but was causing sporadic failures in some of our
pipeline tests with the intel compiler (essentially due to initialized
memory when openMP is used) for subsequent commits.

This patch merges the internal `save_restart` function with the new
`save_MOM6_internal_state` function into a new general MOM restart
function.

It also makes an effort to eliminate `MOM_restart` as a driver
dependency, narrowing the required MOM API for existing and future
drivers.

Also removes the `restart_CSp` argument from `MOM_wave_interface_init`,
since it appeared to be used for nothing.

MOM simulations typically abort of the restart directory (usually
RESTART) are absent.  This patch adds POSIX support for mkdir() and
creates the directory if it is missing.

Using inquire() to check for directory existence is not possible, since
at least one compiler (Intel) does not consider directories to be files.

The inquire call is replaced with a C interface to the POSIX stat()
function.  We do not fully emulate the behavior of stat, but we use its
return value to determine existence of directories.  This provides a
more reliable method for identifying the existence of the directory.

This should resolve many of the observed problems with RESTART creation
in coupled runs.

  Cancelled out factors of GV%Z_to_H in MOM_hor_visc.F90 to simplify the code
and reduce the dependence on the value of GV%Rho_0 in non-Boussinesq mode.  This
involved changing the units of 3 internal variables in horizontal_viscosity and
one element in the hor_visc_CS type to use thickness units or their inverse.
Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode.

  Revised the units of 12 vertvisc_type elements to be based on thicknesses, so
that vertical viscosities (in [H Z T-1 ~> m2 s-1 or Pa s]) are stored as dynamic
viscosites when in non-Boussinesq mode, with analogous changes to the diapycanl
diffusivity (now in [H Z T-1 ~> m2 s-1 or kg m-1 s-1]).  Similarly changed the
units of the 2 Rayleigh drag velocity elements (Ray_u and Ray_v) of the
vertvisc_type from vertical velocity units to thickness flux units and to more
accurately reflect the nature of these fields.  The bottom boundary layer TKE
source element (TKE_BBL) was also revised to [H Z2 T-3 ~> m3 s-3 or W m-2].
This commit also adds required changes to the units of the viscosities or
shear-driven diffusivities returned from KPP_calculate, calculate_CVMix_shear,
calculate_CVMix_conv, Calculate_kappa_shear, Calc_kappa_shear_vertex,
calculate_tidal_mixing and calculate_CVMix_tidal.

 Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.

  Rescaled diapycnal diffusivities passed as arguments in non-Boussinesq mode,
to be equivalent to the thermal conductivity divided by the heat capacity,
analogously to the difference between a kinematic viscosity and a dynamic
viscosity, so that the new units are [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

  This includes changing the units of 4 arguments to set_diffusivity;
3 arguments each to calculate_bkgnd_mixing, add_drag_diffusivity,
add_LOTW_BBL_diffusivity, user_change_diff, calculate_tidal_mixing and
add_int_tide_diffusivity; 2 arguments to KPP_calculate, calculate_CVMix_conv,
compute_ddiff_coeffs, differential_diffuse_T_S, entrainment_diffusive,
double_diffusion, add_MLrad_diffusivity, and calculate_CVMix_tidal; and one
argument to energetic_PBL.

  The units of 36 internal variables were also changed, as were a total of
29 elements in various opaque types, including 8 elements in bkgnd_mixing_cs,
2 in diabatic_CC, 3 in tidal_mixing_diags type, 1 in user_change_diff_CS,
9 in set_diffusivity_CS type, and 6 elements in diffusivity_diags.

  Two new internal variables were added, and several redundant GV%H_to_Z
conversion factors were also cancelled out, some using that
GV%H_to_Z*GV%Rho0 = GV%H_to_RZ.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.

  Changed the units for TKE arguments to [H Z2 T-3 ~> m3 s-3 or W m-2] for
find_TKE_to_Kd, add_drag_diffusivity, calculate_tidal_mixing and
add_int_tide_diffusivity, with similar changes to the units of 21 diagnostics
or internal variables in the same routines and in add_LOTW_BBL_diffusivity and
add_MLrad_diffusivity.  Dozens of unit conversion factors were also cancelled
out with these changes, including using that GV%Z_to_H/GV%Rho_0 = GV%RZ_to_H and
that GV%Rho_0*GV%H_to_Z = GV%H_to_RZ for both Boussinesq or non-Boussinesq
configurations.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.

  Revised the Rossby_front initialization routines to work directly in thickness
units and added completely separate algorithms to initialize the Rossby_front
thicknesses and velocities consistently when the Boussinesq approximation is not
being made.

   To accommodate this change, error handling was added to detect when the
THICKNESS_CONFIG and TS_CONFIG settings are incompatible.  As a part of this
commit the units of the h arguments to Rossby_front_initialize_thickness and
Rossby_front_initialize_temperature_salinity are changed. MAXIMUM_DEPTH is now
read in and rescaled via get_param in the Rossby_front routines rather than
simply being pulled from the ocean grid type.

  There are are also changes to the units of 13 internal variables and 14 new
internal variables in the Rossby_front routines.

  Also pass max_depth as a new argument the internal function Hml to replace the
use of G%max_depth and permit simpler changes to the units being worked with.

  All answers are bitwise identical in Boussinesq mode, but there are
substantial changes (improvements?) in answers in non-Boussinesq mode that are
now independent of the value of the Boussinesq reference density.  There are
changes to the units of arguments to two publicly visible routines.

  Corrected the logic of a recently added test for inconsistently specified
thickness units during the initialization of a new run to only apply to a new
run.  This was causing an incorrect fatal error with some restart tests.  All
answers are bitwise identical in cases that worked previously, but it corrects
the problem with the restarted cases that had been aborted by the test that
was added with the revision to the Rossby_front_2d_initialization code.

  Added the new subroutine find_rho_bottom to return a 1-d slice of the in situ
density averaged over a specified distance from the bottom when in fully non-
Boussinesq mode or an array filled with the Boussinesq reference density.  This
new routine is not yet used with this commit, but it has been tested in another
commit that will follow shortly.  All answers are bitwise identical, but there
is a new public interface.

  Use the in situ near bottom density to calculate the internal tide drag,
energy input and energy loss terms when in non-Boussinesq mode.  This change
includes the addition of an argument containing the near-bottom density to
propagate_int_tide, itidal_lowmode_loss and find_N2_bottom.  The recently added
routine find_rho_bottom is used to calculate this near-bottom density.

  Several instances where the Boussinesq reference density or GV%Z_to_H were
used have been eliminated from use in non-Boussinesq cases by this change,  to
simplify the code and reduce the dependence on the value of GV%Rho_0 in
non-Boussinesq mode.  This involved changing the units of 4 internal variables
in find_N2_bottom to use thickness units or related units.  In some places,
GV%Rho0 was replaced with GV%H_to_RZ.  It also includes the rescaling of a
variable in int_tide_CS, and a new element with the bottom drag in the
int_tide_input_type.

 All answers are bitwise identical in Boussinesq mode, but some solutions will
change in non-Boussinesq mode with this change, and there are new arguments to
publicly visible subroutines and a new element in a transparent type.

  Set a minimum answer date of 20230701 when the model is known to be in
non-Boussinesq mode, and ignore all ANSWERS_2018 flags in that same mode.  This
change only applies to the code called from the src directories, and not to code
within or called directly from the config_src/drivers directories or from
another component model like SIS2 where the information about whether the model
is in non-Boussinesq mode is not available.  A verticalGrid_type arguemnt was
added to porous_barriers_init to support one of these changes.  This commit will
change answers in some non-Boussinesq configurations, but all answers in
existing test cases appear to be bitwise identical.  There are fewer entries
logged in non-Boussinesq MOM_parameter_doc files.

MOM_file_parser's open_param_file() contains explicit inquire() calls
when assessing the correctness of opening such a file.  As written,
these could be called by any rank, and are not thread safe.  In rare
cases (usually related to testing), this would cause a race condition
and raise an error.

Even ignoring the errors, it is probably better if only one rank makes
these calls, rather than all of them.

The following patch modifies the function so that only root PE invokes
inquire().

There is not much to celebrate about this patch; it does not try to
clean up the intrinsic weirdness of the IO handling.  But it does appear
to fix some of the most apparent problems.

  The surface gravity waves invariably work with depths in meters to match the
units of the horizontal wavenumbers, even in non-Boussinesq mode, so several of
the variables that are passed to or used in the MOM_wave_interface module have
been revised to use units of [Z ~> m] and not [H ~> m or kg m-2].  There are
changes to the names or units of a total of 8 arguments to Stokes_PGF,
Update_Stokes_Drift, Get_StokesSL_LiFoxKemper, Get_SL_Average_Prof and
get_Langmuir_Number, and a new argument to StokesMixing.  To accommodate these
changes, there are now calls to thickness_to_dz that precede the
Update_Stokes_Drift or Stokes_PGF calls and a new 3d array in step_MOM.

  Additionally, four hard-coded dimensional constants in Stokes_PGF were given
the needed scaling factors to pass dimensional consistency testing.  This change
also required the addition of a unit_scale_type argument to Stokes_PGF, and
corresponding changes to the calls to this routine from step_MOM_dyn_split_RK2.
Three comments pointing out probable bugs or instances of inaccurate algorithms
were also added.  Incorrect units were also fixed in the comments describing one
internal variable in Update_Stokes_Drift.

  Also added the new runtime parameter RHO_SFC_WAVES to set the surface seawater
density that is used in comparison with the typical density of air set by
RHO_AIR to estimate the 10-meter wind speed from the ocean friction velocity and
the COARE 3.5 stress relationships inside of get_StokesSL_LiFoxKemper() when the
Li and Fox-Kemper (2017) wave model (i.e. LF17) is used.  As a part of this
change, there is a new verticalGrid_type argument to the internal routine
set_LF17_wave_params().  By default, RHO_SFC_WAVES is set to RHO_0.

  Other specific changes include:

 - Changed the units of the publicly visible KvS element of the
   wave_parameters_CS from [Z2 T-1] to [H Z T-1]

 - Replaced the layer thickness argument (h with units of [H ~> m or kg m-2]) to
   Update_Stokes_Drift(), Get_SL_Average_Prof(), get_Langmuir_Number() and
   Stokes_PGF() with a vertical layer extent argument (dz in units of [Z ~> m])

 - Add a vertical layer extent argument (dz) to StokesMixing()

 - Add a unit_scale_type argument to Stokes_PGF()

 - Multiply a hard-coded length scale in Stokes_PGF by the correct unit scaling
   factor and added comments to highlight these hard-coded lengths

 - Corrected the unit description of the DecayScale internal variable

 - Changed the units of 3 internal variables and added 1 new internal variable
   in StokesMixing

 - Added two new arrays to step_MOM for the layer vertical extent and friction
   velocity being set via thickness_to_dz and find_ustar and being passed to
   Update_Stokes_Drift

 - Added new arrays for the layer vertical extent and calls to thickness_to_dz
   in KPP_compute_BLD and energetic_PBL, as well as a new dz argument to
   ePBL_column to provide an altered argument to get_Langmuir_Number

 - Use find_ustar to set the friction velocity passed to Update_Stokes_Drift

 - Update fluxes%ustar or fluxes%tau_mag halos as necessary when the waves are
   in use

  A total of 28 unit conversion factors were eliminated with these changes.

  All answers are bitwise identical in Boussinesq mode, but they are changed in
non-Boussinesq mode by the use of the layer specific volume to convert between
thicknesses and vertical extents.  The units of several arguments to publicly
visible routines are altered as is a publicly visible element in
wave_parameters_CS and there is a new runtime parameter that appears in some
MOM_parameter_doc files.

Two fixes to the ice shelf melt parameterization in MOM_ice_shelf.F90:

(1) the removal of an extra von Kármán constant which differs from the Holland and Jenkins (1999) formulation. My working is detailed in this document.

(2) a modification to the it3 loop in shelf_calc_flux subroutine, which currently does not iterate because wB_flux does not get updated to its new value via the Newton solver method. Instead, the loop either runs 30 times with the same value, or is below the threshold and exits on the first loop. I added a line to update the value of wB_flux.

Specific changes include:

* Remove KV in n_star_term definition

* Add line to ice shelf param it3 to correct iteration

* Reinstate von Karman constant and rearrange to remove double division

* Remove unneeded salt iteration line that overrides option of using false position method

* Change value of ZETA_N to be consistent with Holland & Jenkins 99

* Fix typo, ZETA_N should be 0.13 not 0.013

* Add parameter for ZETA_N

* Make RC and VK runtime parameters

* Add runtime parameters to fix buoyancy iteration bug and salt overwriting false position method bug

* Add runtime parameter for buoyancy iteration Newton's method convergence criteria

  Use SpV_avg in both hybgen_unmix and regridding_main to estimate the nominal
ocean bottom depth in thickness units when in fully non-Boussinesq mode.  This
change eliminates certain dependencies on the Boussinesq reference density via
GV%Z_to_H.  Also added a call to calc_derived_thermo in ALE_regrid_accelerated
that is necessary for the specific volume to be updated, along with tests for
the validity of the specific volume with the 1-point halos that are currently
used with ALE regridding and remapping.

  Also, use the total thickness place of the nominal depth in build_grid_rho and
build_grid_sigma when in fully non-Boussinesq mode.  This is mathematically
equivalent but changes answers at roundoff.

  All answers are bitwise identical in Boussinesq mode, but ALE-based solutions
change in fully non-Boussinesq mode with this commit.

  This commit completes the non-Boussinesq revision of the diabatic_driver code.
The revised code uses thickness_to_dz and works with internal variables in
vertical distances in the denominator of diffusive flux calculations in
diabatic_ALE_legacy, diabatic_ALE and layered_diabatic.  The code now uses
find_ustar to extract the friction velocities passed to KPP_compute.  The (tiny)
boundary layer tracer diffusivity is also rescaled to [H2 T-1].  With this set
of changes, all implicit references to Boussinesq reference density are
eliminated from the calculations in diabatic_driver when in non-Boussinesq mode.

  A total of 14 thickness rescaling factors were cancelled out, and there are 15
new or renamed variables in the diabatic routines.  4 new checksum calls are
also added to diabatic_ALE_legacy to assist in debugging.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode the use of the layer
averaged specific volume to convert thicknesses to vertical distances leads to
changing answers.

  Ignore SURFACE_ANSWER_DATE in non-Boussinesq mode and always allocate
tv%SpV_avg in fully non-Boussinesq mode, setting it to the inverse of the layer
densities in calc_derived_thermo() if there is no equation of state.  Also when
in fully non-Boussinesq mode cancelled out some rescaling factors in
dz_to_thickness_EOS and dz_to_thickness_tv.  Also revised dz_to_thickness_simple
in non-Boussinesq and non-layered mode to use RHO_KV_CONVERT instead of RHO_0 to
rescale vertical distances to thicknesses.  Also set the default value of
CALC_RHO_FOR_SEA_LEVEL to true when fully non-Boussinesq.  All Boussinesq
answers are bitwise identical, but some non-Boussinesq answers do change and
become less dependent on the Boussinesq reference density.  Because
SURFACE_ANSWER_DATE is no longer being used in non-Boussinesq mode, is is no
longer being logged in the MOM_parameter_doc files for these experiments.

  This commit revises full_convection to avoid any dependency on the Boussinesq
reference density when in non-Boussinesq mode.  Specifically, it changes the
units of the Kddt_smooth argument to full_convection and its counterpart in
smoothed_dRdT_dRdS to use units of [H Z ~> m2 or kg s-1], which becomes a mass
based diffusivity when in non-Boussinesq mode.  This change also uses vertical
distances for internal calculations in smoothed_dRdT_dRdS, and includes a call
to thickness_to_dz in the full_convection routine.

  This commit also revises the unit conversion of the (absurdly large) mixing
length in unstable points in full_convection and of the (tiny) added smoothing
distance used in the denominators of some expressions in smoothed_dRdT_dRdS to
avoid division by zero with massless layers so that they are both based on the
density used in rescaling input parameters (RHO_KV_CONVERT), and do not depend
directly on the Boussinesq reference density (RHO_0).  These parameters would
have been set via get_param calls, but there is no control structure for
full_convection parameters.

  All Boussinesq answers are bitwise identical, but non-Boussinesq answers are
altered by the use of the layer specific volumes, rather than the Boussinesq
reference density, to convert layer thicknesses into vertical distances.  This
commit includes a change to the units of an argument (Kddt_smooth) to a publicly
visible interface (full_convection).

  Revised the code in MOM_vertvisc.F90 to work internally with thickness-based
units for the viscous coupling between layers, which eliminates any dependence
on the value of GV%Rho_0 in non-Boussinesq mode.  Various kinematic viscosities
are replaced with dynamic viscosities in non-Boussinesq configurations,
including revising the scaled units of the viscosities to [H Z T-1 ~> m2 s-1 or
Pa s].  This commit also modifies the code to explicitly use vertical distances
rather than thicknesses when calculating the vertical viscous coupling
coefficients in vertvisc_coef and find_coupling_coef.

  This commit changes the units of numerous variables to use thickness, vertical
distance, dynamic viscosity or other related units, including:

 - 14 elements in the vertvisc_CS type
 - 6 arguments to the private routine find_coupling_coef
 - 2 arguments to the private routine find_coupling_coef_gl90
 - 1 arguments ("a") to the public routine write_u_accel
 - 1 arguments ("a") to the public routine write_v_accel
 - 1 internal variable in vertvisc
 - 1 internal variable in vertvisc_remnant
 - 23 internal variables in vertvisc_coef,
 - 7+4+4+3 internal variables in find_coupling_coef,
 - 1 internal variable in find_coupling_coef_gl90

  Local variables that are no longer needed were eliminated in vertvisc and
vertvisc_remnant, while 2 new local variables were added to find_coupling_coef
and 6 new local variables were added to vertvisc_coef. In 6 places the
Boussinesq reference density was replaced with GV%H_to_RZ, which is equivalent
to the reference density in Boussinesq mode, but scales to 1 in non-Boussinesq
mode.

  The previous dimensional rescaling factor for KD_GL90 was incorrect (and
inconsistent with the correct scaling factor used when reading in the analagous
spatially varying kappa_gl90_2d); this has been corrected in this commit.

  A total of 59 GV%H_to_Z or GV%Z_to_H unit conversion factors or references to
GV%Rho_0 were eliminated with these changes, and in non-Boussinesq mode there is
no longer any dependence on the Boussinesq reference density.

  Replaced the forces argument to find_coupling_coef with an array of the
friction velocities and use find_ustar to set them.  When in non-Boussinesq
mode, this has the effect of using forces%tau_mag and tv%SpV_avg instead of
forces%ustar and GV%Rho0 when interpolating the friction velocity and stress
magnitude in find_coupling_coef.

  Revised units used to set the GL90 viscosities and rescale to convert
diagnostics of the vertical viscosities in the MOM_vert_friction module, so that
they do not depend on RHO_0 when in non-Boussinesq mode.

  To accomodate this change in vertvisc, the units of the "a" arguments to
write_u_accel and write_v_accel were also changed to use thickness-based
arguments.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode.  In non-Boussinesq mode, the answers are changed
by the replacement of the Boussinesq reference density by expressions using the
layer-averaged specific volumes.  This commit changes the units of 2 arguments
in public (diagnostic) subroutine interfaces.

  Added the new tau_mag_gustless element to the forcing type, and set this array
in set_derived_forcing_fields if it is associated.  Also added new tau_mag
optional arguments to the public routines allocate_forcing_by_group() and
allocate_mech_forcing_by_group(), with similar mandatory arguments added to the
private subroutines get_forcing_groups() and get_mech_forcing_groups().  Tests
for an associated pointer were added before all calls setting the ustar, tau_mag
or ustar_gustless arrays, preparing for these ustar pointers only to be set when
the model is run in Boussinesq mode when the ustar actually make sense.

  Also use specific volume derivatives to calculate non-Boussinesq mode buoyancy
fluxes in calculateBuoyancy_Flux1d, leaving the Boussinesq buoyancy flux
calculations unchanged.

  All answers are bitwise identical in Boussinesq mode, but there is a new
element in a transparent type and new optional arguments publicly visible
subroutines, and answers will change in non-Boussinesq cases that depend on the
surface buoyancyFlux that is returned by calculateBuoyancy_Flux1d.

  This commit has a set of changes that set the tidal mixing in non-Boussinesq
mode without reference to the Boussinesq reference density, by converting
thicknesses to depths via a call to thickness_to_dz, by calculating an average
near-bottom density explicitly via a call to find_rho_bottom, and by working
with an appropriate mix of layer vertical extents and thicknesses.

  The contents of find_N2, calculate_tidal_mixing, calculate_CVMix_tidal and
add_int_tide_diffusivity are all revised to work directly with vertical
distances, set via thickness_to_dz, rather than rescaled thicknesses.  This
includes replacing the 3-d layer thickness argument to calculate_tidal_mixing,
calculate_CVMix_tidal and add_int_tide_diffusivity with a 2-d slice of vertical
distances across layers.  For find_N2, the 2d-slice of vertical distances is a
new argument.

  The revised code also uses the in situ bottom density when adding certain
contributions to non-Boussinesq diffusivities.  This change includes the
addition of a new bottom density argument to find_N2 and calculate_tidal_mixing,
along with a call to find_rho_bottom to set this density.  There is also a new
runtime parameter, DZ_BBL_AVG_MIN, stored in set_diffusivity_CS, that helps
determine the vertical extent of that average.

  There are a total of 6 new (or renamed) internal variables, and 6 new
subroutine arguments, with revisions to the rank of one internal variable. The
units of 2 internal variables and 1 element of the tidal_mixing_cs type are
altered.  This commit eliminates 22 factors of GV%H_to_Z and two direct uses
of GV%Rho0.

  All answers are bitwise identical in Boussinesq mode, but solutions will
change in non-Boussinesq mode with this commit.  There are several new arguments
to a publicly visible subroutine, and a new runtime parameter.

  Added the name to the do_taumag argument in a call to allocate_forcing_type
in allocate_forcing_by_ref to account for the fact that there are unused
wave-related optional arguments in this interface.  When this was omitted in
the current code, the wrong arrays are being allocated during rotation tests
with resultant segmentation faults in those tests.  This commit corrects a bug
that was recently added with MOM6 dev/gfdl PR mom-ocean#445.

  Add missing dimensional rescaling factors for hard-coded temperature and
salinity offsets in debugging conservation checks in regularize_surface.  All
answers are bitwise identical, but without this change some configurations can
erroneously give fatal errors for some rescaling settings when DEBUG=True.

  This commit revises MOM_CVMix_conv, MOM_CVMix_ddiff and MOM_CVMix_shear to
work in an appropriate mixture of thickness and vertical extent variables to
enable their use in non-Boussinesq mode, using thickness_to_dz to convert
between the two, while retaining the previous answers in Boussinesq mode.  This
includes the use of a layer thicknesses rather than a vertical distance in the
denominator of the calculation of the buoyancy frequency and to the units of the
internal rescaled gravity variable in CVMix_conv.  This commit eliminates any
direct or indirect dependency on the Boussinesq reference density in these 3
modules when in non-Boussinesq mode.

  This set of changes includes changing the units of 7 internal variables and
the addition, renaming, or change in the dimensions of 4 internal variables.
One element in the CVMix_ddiff_cs type is rescaled to use thickness units.
These changes lead to the removal of 4 factors of GV%H_to_Z in these 3 modules.

  Answers will change in non-Boussinesq mode when USE_CVMix_CONVECTION,
USE_CVMIX_DDIFF, USE_PP81 or USE_LMD94 are true, but they are bitwise identical
in all Boussinesq test cases.

  This commit revises MOM_CVMix_KPP to work in an appropriate mixture of
thickness and vertical extent variables to enable their use in non-Boussinesq
mode, using thickness_to_dz to convert between the two, while retaining the
previous answers in Boussinesq mode.  This includes the use of a layer
thicknesses rather than a vertical distance in the denominator of the
calculation of the buoyancy frequency and the replacement of the layer thickness
argument (h) to KPP_smooth_BLD with a layer vertical extent (dz).  When in
non-Boussinesq mode, the buoyancy difference between layers is normalized by the
average of the density of the two layers rather than the Boussinesq reference
density. This commit eliminates any direct or indirect dependency on the
Boussinesq reference density in CVMix_KPP when in non-Boussinesq mode.

  This set of changes includes changing the units of 1 internal variables and
the addition of 2 new internal variables, and a change to the name and units of
one argument to KPP_smooth_BLD.  These changes lead to the removal of 5 factors
of GV%H_to_Z.

  Answers will change in non-Boussinesq mode when USE_KPP is true, but they are
bitwise identical in all Boussinesq test cases.

- Tracer advection diagnostics in symmetric mode had errors at processor boundaries because they were computed over only the tracer indices. The tracer advection diagnostics have been updated to compute over a range of indices identical to the fluxes, which allows it to also work on the u/v point indices in both non-symmetric and symmetric mode.
- These code updates only impact the tracer diagnostics, and in no cases were the model solutions impacted by the bug outside of these select diagnostic fields.
- The corrected diagnostics reproduce the old non-symmetric mode results in both symmetric and non-symmetric mode.
- Some logic has been updated and comments were added to improve the code, but none of these changes impact the model solution.

  Refactored interface_filter when in non-Boussinesq mode to avoid any
dependencies on the Boussinesq reference density, and to translate the volume
streamfunction into the mass streamfunction using an appropriately defined
in-situ density averaged to the interfaces at velocity points.  This form is
similar to the one that is used in thickness_diffuse.

  No public interfaces are changed, and all answers are bitwise identical in
Boussinesq or semiBoussinesq mode, but they will change in non-Boussinesq mode
cases that use the interface filter.

  Revise calc_isoneutral_slopes to eliminate any dependence on the Boussinesq
reference density when in non-Boussinesq mode.  This includes the addition of
two new arrays to hold the rescaled product of the gravitational acceleration
and the specific volume interpolated to the interfaces at velocity points.
The answers change in non-Boussinesq mode when one of several parameterizations
that use the isoneutral slopes are in use, but are bitwise identical in
Boussinesq or semi-Boussinesq mode.

  Rescaled CFC flux diagnostics in the CFC cap via a conversion factor in their
register diag_field calls, rather than by doing array syntax math in their
post_data calls.  All answers are bitwise identical.

  This commit revises set_viscous_BBL and set_viscous_ML to work in fully
non-Boussinesq mode, eliminating all dependencies on the Boussinesq reference
density when in non-Boussinesq mode.  These changes include using the specific
volume interpolated to velocity points for the conversion between thickness and
height changes when tv%SpV_avg is allocated, including in the calculation of the
rescaled ustar in the top and bottom boundary layers and in the Rayleigh drag
calculation.  In set_viscous_ML, the derivatives of specific volume are used to
determine the reduced gravity across interfaces when in non-Boussinesq mode.

  Both set_viscous_BBL and set_viscous_ML now work extensively in height units
rather than thickness units in various places where it is more appropriate to do
so when in non-Boussinesq mode, and both use calls to thickness_to_dz to convert
between thicknesses and vertical distances.  In some places this leads to extra
calculations using separate arrays that are rescaled duplicates of each other in
Boussinesq mode, which will probably slow the model down a little.

  There is one rescaling bug due to a hard-coded unrescaled dimensional constant
that was fixed.  It occurs when setting a tiny floor in a thickness in one
inverse calculation.  However, this appears to be so small that it does not
change any answers in the MOM6-examples test suite, and could only do so if
ANGSTROM is set to be many orders of magnitude less than the typical value of
1e-10 m, so it was decided that no runtime bug-flag is needed in this case.

  A call to find_star is now used to specify the friction velocity used in
set_viscous_ML.  When in non-Boussinesq mode, this has the effect of using
forces%tau_mag and tv%SpV_avg instead of forces%ustar and GV%Rho0 to determine
the friction velocity.  In set_viscous_BBL the layer specific volumes are used
in non-Boussinesq mode to find the average velocities used in the linearization
of the nonlinear bottom drag, with similar changes in set_viscous_ML for the
drag under ice-shelves.

  When in non-Boussinesq mode, the units of Kd_shear, Kv_shear and Kv_shear_Bu
in the restart files are scaled to their natural MKS units of [Pa s] or [kg m-1
s-1] rather than [m2 s-1] to avoid round-off level changes across restarts, with
the units in the files documenting these changes.  These altered units do not
apply to Boussinesq mode restart files.

  In set_visc_init, the rescaling factor used to set CS%Hbbl was changed from
GV%Z_to_H to (US%Z_to_m*GV%m_to_H) so that it does not depend directly on RHO_0.

  These changes are extensive but localized to this file, with one new element
of the opaque set_visc_CS type and 6 elements with altered units, 37 new or
renamed internal variables, and 29 existing internal variables that have revised
units.  In particular one extensively used thickness curvature variable that had
been called just "a" was renamed "crv" to more accurately reflect its purpose
and to make it easier to find when rescaling it, with similar changes to two
other closely related variables.  These changes include the elimination of 29
rescaling factors between thickness and height units, as were two spots that
inappropriately (and unnecessarily, as it turns out) used GV%Rho0 even when in
non-Boussinesq mode.  There are also 4 new debugging checksum output calls.

  In Boussinesq or semiBoussinesq mode, no public interfaces are changed and all
answers are bitwise identical, but in non-Boussinesq mode answers will change
and there are changes to the units of 3 variables in restart files.

  This commit revises the internal routines called by set_diffusivity to work in
fully non-Boussinesq mode, eliminating all dependencies on the Boussinesq
reference density when in non-Boussinesq mode.  The publicly visible interfaces
to this module and the external routines it calls have already been revised, so
only this file needs to be updated.  The specific changes include:

 - Refactored add_LOTW_BBL_diffusivity, add_MLrad_diffusivity, and set_BBL_TKE
   to work in units of layer vertical extent rather than layer thickness to give
   results in non-Boussinesq mode that avoid dependence on the Boussinesq
   reference density.

 - Work with internal variables in vertical distances in the denominator of
   diffusive flux calculations in find_TKE_to_Kd, while other expressions in
   this routine are recast in terms of thicknesses to avoid conversions.

 - Use tv%SpV_avg instead of 1/RHO_0 in find_TKE_to_Kd when in fully
   non-Boussinesq mode.

 - Use layer target density differences in place of g_prime in
   set_density_ratios in mathematically equivalent expressions when
   non-Boussinesq.

 - Use thickness_to_dz in 3 places to convert layer thicknesses into vertical
   distances.

 - The thickness argument to add_MLrad_diffusivity (in [H ~> m or kg m-2]) has
   been replaced with a vertical extent argument (in [Z ~> m]).

 - Use fluxes%tau_mag in place of fluxes%ustar in add_MLrad_diffusivity when in
   non-Boussinesq or semi-Boussinesq mode.  There is a new thermo_var_ptrs type
   argument to the internal routine add_MLrad_diffusivity to permit this
   changes.

 - Use the in situ near-bottom density when adding certain contributions to
   non-Boussinesq diffusivities.  This change includes the addition of a new
   bottom density (rho_bot) argument to add_int_tide_diffusivity,
   add_LOTW_BBL_diffusivity and add_drag_diffusivity.

 - Use GV%dZ_subroundoff in 4 spots in place of GV%H_to_Z*GV%H_subroundoff.

 - A long-standing comment questioning whether there is double-counting of tidal
   mixing has been addressed (there is not) and the comment has been revised
   accordingly.

  These changes involved changing the units of 21 internal variables and 1
element in the set_diffusivity_CS type.  There are 11 new internal variables,
while 9 internal variables were eliminated.  A total of 44 thickness rescaling
factors were eliminated, and there are 2 places where GV%Rho_0 was being used
explicitly that were changed into equivalent rescaling factors.

  All answers are bitwise identical in Boussinesq mode, but some solutions will
change in non-Boussinesq mode with this commit.

  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.

  Added the new runtime parameters RESTORE_FLUX_RHO and TKE_TIDAL_RHO to specify
the densities that are used to convert the piston velocities into restoring heat
or salt fluxes and to translate tidal velocities into tidal TKE input.  Their
defaults are set to RHO_0 to reproduce previous answers.  Also added tau_mag
arguments to 2 calls to allocate_forcing_type() and 4 calls to
allocate_mech_forcing() in the FMS_cap and solo_driver code.  There are new
rho_restore elements in the FMS and solo_driver versions of surface_forcing_CS
and in MESO_surface_forcing_CS and user_surface_forcing_CS.  By default, all
answers are bitwise identical, but there are new entries in some
MOM_parameter_doc files.

  Add tau_mag argument to calls to allocate_forcing_type() in
initialize_ice_shelf_fluxes and the mct_cap and noupc_cap versions of
convert_IOB_to_fluxes and to calls to allocate_mech_forcing() in
initialize_ice_shelf_forces and the mct_cap and noupc_cap versions of
convert_IOB_to_forces.  All answers are bitwise identical.

  Calculate forces%tau_mag in idealized_hurricane_wind_forcing() and
SCM_idealized_hurricane_wind_forcing(), and call allocate_mech_forcing with
the new tau_mag argument so that this array is sure to be allocated.  All
answers are bitwise identical in existing test cases, but this step was
necessary for this code to work in fully non-Boussinesq configurations.

  Use RESTORE_FLUX_RHO in SCM_CVMix_tests and dumbbell_surface_forcing to
specify the density that are used to convert the piston velocities into
restoring heat or salt fluxes.  As with other analogous changes, the default is
set to RHO_0 to reproduce previous answers.  Also set forces%tau_mag in
SCM_CVMix_tests_wind_forcing if it is associated.  There is a new rho_restore
element in the control structures for the SCM_CVMix_tests two module, while the
units of an element in the dumbbell_surface_forcing module are changed. By
default, all existing answers are bitwise identical, but there are new entries
in some MOM_parameter_doc files.

The units of S_REF were first registered in "PSU" and then later by
benchmark in "ppt" which led to a WARNING. Although there is no uniformity
to the units of S_REF throughout the user code (currently PSU, ppt, and 1e-3)
are all used, benchmark is the only one that leads to a warning in the current
suite and changing all to be consistent would 1) unnecessarily update the doc
files, and 2) not be correct for all models. Bottom line, I'm punting on
the best way to handle units of salinity. This commit resolves the only WARNING
we currently get in the MOM6-examples suite.

- Update for missing halo update of field wpup_filtered in MOM_mixed_layer_restrat.F90 during initialization.
- The missing halo update caused the model to diverge when running from restart files.
- With the halo update the model now reproduces from restart files.

  Revised the code in MOM_mixed_layer_restrat.F90 to work internally with
thickness-based units for the restratification timescale calculation and other
internal calculations, which eliminates the dependence in this module on the
value of the Boussinesq reference density in non-Boussinesq mode.  Several other
minor issues (which might not change any answers in production runs) were also
fixed.  The changes with this commit include:

 - When in non-Boussinesq mode, the mixed layer buoyancy gradients are
   determined from the average specific volume referenced to the surface
   pressure, rather than from the average potential density.

 - Use find_ustar to set the friction velocities in the appropriate units in the
   various mixed_layer_restrat routines.

 - A logical branch was added based on the correct mask for land or OBC points
   to avoid potentially ill-defined calculations of the magnitude of the Bodner
   parameterization streamfunction, some which were leading to NaNs.

 - Set a tiny but nonzero default value for MIN_WSTAR2 to avoid NaNs in some
   calculations of the streamfunction magnitude.

 - Within the function mu, the expression for dd was revised in a mathematically
   equivalent way to avoid any possibility of taking a fractional exponential
   power of a tiny negative number due to truncation errors, which was leading
   to NaNs in some cases while developing and debugging the other changes in
   this commit.  This does not appear to change any answers in the existing test
   cases, perhaps because the mixed layer restratification "tail" is not being
   activated by setting TAIL_DH to be larger than 0.

 - The addition of code to both mixedlayer_restrat_Bodner and
   mixedlayer_restrat_OM4 to determine the mixed layer thickness from its
   vertical extent when in non-Boussinesq mode.

  This commit includes changes to the units of the Kv_restrat, ustar_min and
wpup_filtered elements in the mixedlayer_restrat_CS type and the units of four
arguments to the private function growth_time.  CS%wpup_filtered also appears in
the restart files generated with some mixed layer restratification settings, and
it is rescaled to units of vertical distance or mass per unit area in the
restart files depending on whether the model is Boussinesq.

  There are 17 new or renamed internal variables, while the units of 21 internal
variables were changed.  19 rescaling factors were cancelled out or replaced.
There are also comments where variable units were corrected or added.

  The rescaling of several chksum calls for thicknesses was modified to
GV%H_to_mks to avoid any dependence on RHO_0 when non-Boussinesq.

  No public interfaces are changed.  All answers are bitwise identical in
Boussinesq mode (at least when TAIL_DH=0.), but solutions will change in
non-Boussinesq mode when mixed layer restratification is enabled, including
changes to the units of a variable in the restart files.

  Remove dependency on the Boussinesq reference density in MOM_kappa_shear when
in fully non-Boussinesq mode.  This is done by using calls to calculate_density
and calculate_specific_vol_derivs to calculate the thermal expansion and haline
contraction coefficients, and by using thickness_to_dz to convert layer
thicknesses to vertical distances.

  A large part of the changes to the kappa_shear code involved refactoring it to
differentiate between thicknesses and geometric heights by working internally
with dynamic viscosities and diffusivities in non-Boussinesq mode, and to keep
thicknesses in thickness units rather than converting them to distances in
height units.  The internal variable dz was renamed to h_lay in
Calculate_kappa_shear and Calc_kappa_shear_vertex for greater clarity.

  The scaling factor for dz_massless is now set to (US%Z_to_m*GV%m_to_H) in
Calculate_kappa_shear to reduce the dependency on the Boussinesq reference
density; this does not impact answers if RHO_0 = RHO_KV_CONVERT.

  The internal changes related to this are extensive.  This commit includes
changing the units of 63 internal variables, and the addition or renaming of 17
internal variables.  There are 2 new arguments to find_kappa_tke (one is
necessary to differentiate thicknesses from vertical distances and while the
other is used to preserve Boussinesq-mode answers), and 1 new argument to
kappa_shear_column plus the renaming of another.  The units of 13 internal
subroutine arguments and 3 elements of the Kappa_shear_CS type change. A total
of 6 thickness unit conversion factors were eliminated with these changes.

  The answers will change in non-Boussinesq mode when the Jackson et al. shear
mixing parameterization is in use, but they are bitwise identical in all
Boussinesq test cases.  There are no changes to external interfaces.

  Use tv%SpV convert thicknesses to vertical distances in apply_sponge when it
is allocated to replace multiplication by GV%H_to_Z, thereby eliminating the
dependence on GV%RHo_0 in this modue in non-Boussinesq mode.  The new internal
array dz_to_h is used to reduce the code duplication as a result of these
changes.  All answers in Boussinesq test cases are bitwise identical, but
answers change in fully non-Boussinesq mode.  In semi-Boussinesq mode answers
are mathematically equivalent but change at roundoff unless RHO_0 is an integer
power of 2.

  Set DOME inflow properties using the average density rather than RHO_0 when in
non-Boussinesq mode.  The existing Boussinesq solutions are unchanged, but this
changes the non-Boussinesq DOME solutions, which no longer depend on the
Boussinesq reference density.

Calculation for self-attraction and loading (SAL) is separated as a new
module (MOM_self_attr_load) from module MOM_tidal_forcing, as SAL is a
process not limited to tides.

The new module includes both online spherical harmonics method and
scalar approximation. Read-in method (TIDAL_SAL_FROM_FILE) is kept in
the tidal forcing module as it is specific to tides. For the iterative
method (USE_PREV_TIDES), the updating part that is tied to the scalar
approximation is moved to the new SAL model, while the read-in
part remains in the tidal_forcing module.

The tidal forcing module now only contains calculations independent from
 the ocean's state and the only input variables is the current time.

* A new parameter CALCULATE_SAL is added, which controls SAL
calculation in PressureForce independent of whether tides is on or not.
The default of CALCULATE_SAL is TIDES to avoid making changes in old
MOM_input.
* For the unplit mode, runtime parameters calculate_SAL and use_tides
are moved from init subroutines to control structures. This allows safe
deallocations of the corresponding module CSs.
* A new control structure for the SAL module is used by the dynamical
cores and pressure force modules.
* For SAL related parameters, their names still incorrectly contain
TIDE or TIDAL. This will be addressed in the following commits.
* A new diagnostic is added in PressureForce to output calculated SAL
fields. Note that the 'e_tidal' diagnostic is unchanged and
still includes SAL field for backward compatibility.
* Subroutine tidal_forcing_sensitivity, which is used by the barotropic
 solver to take into account the scalar approximation, is renamed to
scalar_SAL_sensitivity.
* Documentations are updated for the cited papers.

The output field eta_tidal from subroutine calc_tidal_forcing originally
 contains both equilibrium tides and tidal SAL that is based on previous
 solutions (TIDAL_SAL_FROM_FILE and USE_PREV_TIDES). This commit
 decompose the two fields for better diagnostics.

This also makes answers using USE_PREV_TIDES is not changed if only one
 tidal constituent is used.

New diagnostics ('e_tidal_eq' and 'e_tidal_sal') are added in pressure
force modules for the decomposed fields. Note that 'e_tidal' still
includes all tidal fields + SAL fields to be backward compatible.

SSH for SAL is modified for grid points with topography above the
reference height z_ref (assumed to be land that can be flooded). Instead
 of eta anomaly referenced to z_ref, eta anomaly referenced to the
 bottom depth is used for these grid points.

Remove "TIDAL" and "TIDE" from the relevant input names of the SAL
module, as SAL is not specific to the tides. This affects both scalar
approximation and the fully online spherical harmonic options.

* For scalar SAL, old parameter names are still acceptable, but a
WARNING is given if these old names appear in MOM_input.
* For read-in SAL, no change is made.
* For iterative method (use_prev_tides), the use of
TIDE_SAL_SCALAR_VALUE is completely deprecated, as this is a feature
that is rarely used.
* For harmonic SAL, a relatively recent feature, a hard obsolete is
applied, i.e. if the old parameters are specified, a FATAL error is
given, unless the new parameters also exist and match the values of the
old parameters.

List of input names changed:
* TIDE_USE_SAL_SCALAR -> USE_SAL_SCALAR
* TIDE_SAL_SCALAR_VALUE -> SAL_SCALAR_VALUE
* TIDAL_SAL_SHT -> USE_SAL_HARMONICS
* TIDAL_SAL_SHT_DEGREE -> SAL_HARMONICS_DEGREE

This commit addresses the issue of bit level answer change with tides
in Boussinesq mode.

* A runtime parameter TIDES_ANSWER_DATE is added to restore old answers
before the SAL module is added. The answer change is due to an reorder
of summations of the SAL and tidal forcing terms.
* The new version flag only applies to the analytical pressure force in
Boussinesq mode, which is virtually the only configuration widely used
and included in the test suite.
* For Montgomery pressure and non-Boussinesq mode, the code is
refactored in a more readable way.

Consistent with the non-Boussinesq option with the new answer option

This commit revises changes of some parameter names from a previous
commit.

* Logical switch TIDE_USE_SAL_SCALAR is obsolete and replaced by
SAL_SCALAR_APPROX.
* TIDE_SAL_SCALAR_VALUE is replaced by SAL_SCALAR_VALUE. The old
parameter name is still accepted but a warning is given.
* Logical switch TIDAL_SAL_SHT is obsolete and replaced by
SAL_HARMONICS.
* TIDAL_SAL_SHT_DEGREE is obsolete and replaced by
SAL_HARMONICS_DEGREE.
* RHO_E is replaced by RHO_SOLID_EARTH.

* An incident of incorrect indent is fixed in SAL_int.
* Parameters read in SAL_int is moved out from if statements. DO_NOT_LOG
is used to prevent logging unused parameters.
* Reading SAL_SCALAR_VALUE in tidal_forcing_init is now consistent with
SAL_init.
* Some unused variables in tidal_forcing_init are removed.

  Use TIDES_ANSWER_DATE with an additional if-block to recover the dev/gfdl
answers in semi-Boussinesq mode runs with tides enabled, similarly to what is
already being done in Boussinesq mode.  This changes some answers for the
nonBous_global and tides_025 test cases back to those that would be obtained
with the code on dev/gfdl.   TIDES_ANSWER_DATE is now logged in all cases with
tides, so the MOM_parameter_doc files change in some non-Boussinesq cases.

  This commit revises the internal code in step_MOM and extract_surface_state
to work in fully non-Boussinesq mode without any dependencies on the Boussinesq
reference density and provide updates needed to reproduce answers across
restarts in non-Boussinesq mode.  The subroutine calls from this module have
already been revised, so only this file needs to be updated.  The specific
changes include:

 - Rescaled the nominal mixed layer depth variables HFREEZE, HMIX_SFC_PROP and
   HMIX_UV_SFC_PROP for determining surface properties to eliminate unit
   conversion factors during the run.  This changes the units of 3 elements
   of the MOM_control_struct and one internal variable and adds 3 new internal
   variables.

 - Call calc_derived_thermo to update SpV_avg at the start of step_MOM if the
   internal ocean state depends on the surface pressure and that surface
   pressure has been passed in from the driver.  There are 3 new internal
   variables and in some cases a new blocking halo update to enable these
   changes.

  All Boussinesq solutions are bitwise identical with this commit.  It changes
answers in some non-Boussinesq fully coupled or ice-ocean configurations.
Other non-Boussinesq mode cases are mathematically equivalent in the case where
RHO_0 = RHO_KV_CONVERT and bitwise identical if RHO_KV_CONVERT is an integer
power of 2.  This commit corrects issues with reproducability across restarts
in the non-Boussinesq cases where the answer changes.

  Revised BFB_set_coord and BFB_buoyancy_forcing to be consistent with other
instances of a linear equation of state, and to set g_prime in non-Boussinesq
mode similarly to how it is set in other places.   Also use RESTORE_FLUX_RHO in
place of RHO_0 when to specify the densities that are used to convert the piston
velocities into restoring heat or salt fluxes with BFB_buoyancy_forcing, like in
other places in the MOM6 code.  This change will change answers by default when
BFB_set_coord is used, but the old Boussinesq answers can be recovered by
setting the reference salinity S_REF to 38.75 ppt or by setting RHO_T0_S0 to
1003.0 kg m-3.

  Use the minimal vertical distance GV%dZ_subroundoff in 8 spots in 2 files in
place of GV%H_to_Z*GV%H_subroundoff.  Four internal variables were renamed for
clarity of purpose in porous_widths_layer and PressureForce_Mont_Bouss.  Also
modified the non-Boussinesq units in a description of two integrated energy
diagnostics from [H L4 T-3 ~> kg m2 s-3] to [H L4 T-3 ~> W]; these units are
equivalent but the latter is more informative.  Answers in all Boussinesq test
cases and the existing non-Boussinesq test case in the MOM6-examples test suite
are bitwise identical.

  Rescale the units diapycnal diffusivity used in the offline tracer mode to
work in [H Z T-1 ~> m2 s-1 or kg m-1 s-1], making it consistent with what has
already been done for other diapycnal diffusivities.  This change involves
changes to the units of arguments to update_offline_from_files and
ALE_offline_inputs and to two elements of offline_transport_CS.   It also
includes a new thermo_vary_type argument to offline_diabatic_ale, as well as a
call to thickness_to_dz in that same routine to set the values of the newly
added 3-d array of the vertical distances across layers and a change in name and
units of the inverse thickness variable used to find the fluxes.  All answers
are bitwise identical in Boussinesq mode, but in offline tracer calculations in
non-Boussinesq mode they will change to become independent of the Boussinesq
reference density.

  Obsoleted 18 2018_ANSWERS parameters in 26 modules, with obsoleting message
hints indicating the corresponding ANSWER_DATE parameters.

  The runtime parameters that have been obsoleted are DEFAULT_2018_ANSWERS,
SURFACE_FORCING_2018_ANSWERS, WIND_GYRES_2018_ANSWERS, BAROTROPIC_2018_ANSWERS,
EPBL_2018_ANSWERS, HOR_REGRID_2018_ANSWERS, HOR_VISC_2018_ANSWERS,
IDL_HURR_2018_ANSWERS, MEKE_GEOMETRIC_2018_ANSWERS, ODA_2018_ANSWERS,
OPTICS_2018_ANSWERS,  REGULARIZE_LAYERS_2018_ANSWERS, REMAPPING_2018_ANSWERS,
SET_DIFF_2018_ANSWERS, SET_VISC_2018_ANSWERS, SURFACE_2018_ANSWERS,
TIDAL_MIXING_2018_ANSWERS and VERT_FRICTION_2018_ANSWERS.

  These changes will cause cases that use these older parameters to fail with a
useful error message, but this change has been discussed and agreed to by the
MOM6 community on July 17, 2023.  All answers are bitwise identical in any cases
that work.

This patch updates the macros used in the `ac/deps` directory.  It
includes a bugfix in the library testing, and also bundles the macro for
testing Fortran links to C libraries.  Both of these macros were already
in `ac/m4` but unused by `ac/deps`.

No problems were observed in any recent FMS builds, although in
principle they could have been affected by these errors.  The issues
were only detected when the macros were used in an unrelated project
(LM3 land model).

These issues are due to code duplication in `ac/m4` and `ac/deps/m4`.
This separation was intended to draw a clear logical separation between
MOM6 and its FMS dependency, and only provide what is required by each
respective build.  But perhaps there are some flaws in this thinking,
and may warrant further discussion.

This commit fixes a combination of bugs related to halo update errors, improperly-defined iteration bounds, and an unexpected change in ice-shelf bed_elev after a restart. These issues caused crashes during the ice-shelf velocity (SSA) solution, and ice shelf dynamics restarts were not bitwise identical. The issues were resolved with the following changes:

-Added bed_elev to the ice shelf restart file (previously, bed_elev was instead calculated incorrectly after a restart). Also fixed the subsequent bed_elev pass_var call, where the halo update was failing due to an error in how an optional argument (CENTER) was passed.

-Added additional pass_var/pass_vector calls upon ice shelf restarts to guarantee that all halos are properly filled.

-Disabled the additional ice shelf flow solve (ice_shelf_solve_outer) after a restart, which was unnecessary and prevented bitwise identical restarts.

-Fixed incorrect bounds for iterations within ice_shelf_solve_inner and for the basis functions in calc_shelf_visc

This commit also sets the mol_wt argument to optional in the dummy function aof_set_coupler_flux (from ice_solo_driver/atmost_ocean_fluxes.F90), as that is how it is called from the FMS coupler. This is required for compilation in ice-shelf-only mode.

* implement restart for internal tides

* add a call to register restart diabatic in MOM

* register restart diabatic call register restart internal tides

* internal tides restart uses the extra_axes optional argument

* support for extra_axes added to MOM_restart/MOM_io

* implement 4d restart/ parallel_restart not working

* add proper documentation for extra axes

* fix read restart for 4d array

* passes global_file optional arg to allow reading
  in parallel restart files

* write one restart variable per vertical mode

  vertical modes >= 6 do not propagate hence there is no need for dynamic array size on this dimension

---------

Co-authored-by: Raphael Dussin <[email protected]>

  This commit revises MOM_MEKE to work in with total depths in thickness units
to avoid any dependence on the Boussinesq reference density in non-Boussinesq
mode, while retaining the previous answers in Boussinesq mode.  It also adds the
new runtime parameter MEKE_TOTAL_DEPTH_RHO to specify the density that is used
to convert the bathymetric depth into a nominal ocean mass in non-Boussinesq mode
when MEKE is enabled and MEKE_FIXED_TOTAL_DEPTH is true.

  There is a new element in the MEKE_CS type, and the scaled units of two other
elements (cdrag and MEKE_min_depth_tot) are altered.  The units of a total of 5
arguments to the private routines MEKE_equilibrium, MEKE_equilibrium_restoring,
MEKE_lengthScales and MEKE_lengthScales_0d where changed consistently with the
broader changes here.  There is also a new vertical_grid type argument to
MEKE_equilibrium_restoring.  The units of 6 internal variables are changed, 2
new internal variable were added and 5 internal variables were eliminated.
There are 4 places where multiplication by a GV%H_to_Z thickness rescaling
factor was eliminated, while 4 instances that use GV%Rho0 directly were
effectively replaced by the use of GV%H_to_RZ.

  All Boussinesq answers are bitwise identical, but answers will change in
non-Boussinesq mode and there is a new runtime parameter the MOM_parameter_doc
files for some non-Boussinesq simulations.

The test.summary rule was causing errors in our Gitlab testing due to
multiple runs (concurrent or otherwise) in the same workspace directory.
This patch removes the WORKSPACE directory variable, and each .testing
run happens in its own directory.

Other minor changes:

- The script to generate the summary was moved out of the Makefile and
  into a separate script.

- Unrelated to these changes, error output was extended from 20 to 40
  lines, to provide more readable backtrace output.

  - This option is defaulted to False to retain previous answers,
    but should be set to True for new experiments in order to make
    the Flather OBC routine consistent with the barotropic solver

  - This option only applies for regional OBC cases with Tides and
    scalar self-attraction and loading

  + Try at fixing issue mom-ocean#476

    - Will change answers for problems with OBCs.

  - Get Matt's patch to compile again

Co-authored-by: Kate Hedstrom <[email protected]>

The current version of makedep uses the existence of program units
(module or program) to determine whether a file contains globally
accessible external functions.  It incorrectly ignores files with
externals it happens to contain both externals and a program unit.

We resolve this by tracking whether we are inside or outside of a
program unit while parsing each line.  We ignore keywords which may
appear in comments, as well as tokens which may contain `function` or
`subroutine` as substrings.  This method is not perfect; there are cases
which will incorrectly report entering a block, but this method appears
sufficient for all of our known codebases.

This issue was detected when attempting to parse the AM2 source.  While
still far from a general solution, it is sufficient to handle a few of
the more challenging cases encountered in AM2.

  Add tests to stop iterating when converged to roundoff in dz_to_thickness_EOS
when in fully non-Boussinesq mode.  Also modified the code to keep track of the
layer thickness directly in non-Boussinesq mode rather than setting the layer
thickness based on a difference between interface pressures.  Boussinesq
answers are bitwise identical, but non-Boussinesq answers change.

  This commit revises the internal routines called by energetic_PBL to work in
fully non-Boussinesq mode, eliminating all dependencies on the Boussinesq
reference density when in non-Boussinesq mode.  The publicly visible interfaces
to this module and the external routines it calls have already been revised, so
only this file needs to be updated.  The specific changes include:

 - Work with both thickness (h) and height change (dz) variables in
   energetic_PBL using thickness_to_dz to translate between the two.  In
   Boussinesq mode, these are a simple rescaling by a constant factor, but in
   non-Boussinesq mode they vary by a factor of the specific volume.  Some
   calculations are unnecessarily duplicated in Boussinesq mode, so this might
   slow the model slightly, but they are not duplicative in non-Boussinesq mode.

 - When in non-Boussinesq mode, use forces%tau_mag and tv%SpV_avg instead of
   forces%ustar and GV%Rho0 to determine the friction velocity and mechanical
   TKE input used in energetic_PBL, and to convert the local TKE to a turbulent
   velocity.  This includes the addition of a new energy input argument and an
   interface specific volume argument to ePBL_column and the removal of an
   unused argument to find_mstar.

 - Use tau_mag_gustless times the in situ specfic volume instead of
   ustar_gustless to calculate the Langmuir number used by ePBL when it is run
   in non-Boussinesq mode.

 - The unused Ustar_mean argument to find_mstar was removed.

  This change involves the addition of 2 new arguments to ePBL_column and
changes to the units of another argument.  There are changes to the units of 2
elements of the energetic_PBL_CS type and 14 internal variables. This change
includes the addition of 13 new internal variables, the removal of 8 internal
variable (most of which were renamed to reflect their new units). With this
change, 9 thickness to height conversion factors have been eliminated, and
GV%Rho0 is only used in Boussinesq mode.

  All answers in Boussinesq mode are bitwise identical, but they change in
non-Boussinesq mode cases that use ePBL, which no longer depend on the value of
the Boussinesq reference density.

  Added the optional argument SpV_avg to extract_optics_slice for use along with
an appropriate value for opacity_scale to convert the units of opacity from [Z-1]
to [H-1] in non-Boussinesq mode without making use of the Boussinesq reference
density.  All Boussinesq answers are bitwise identical, but non-Boussinesq
answers will change and become less dependent on the Boussinesq reference density
when this new argument is used.  There is a new optional argument to a publicly
visible subroutine.

  This commit includes a series of distinct changes that enable the use of the
bulk mixed layer code in non-Boussinesq mode, including an option to do the
non-Boussinesq energetic calculations even when the model itself is in
Boussinesq or semi-Boussinesq mode.  When in fully non-Boussinesq mode, there is
no longer any dependence on the Boussinesq reference density.

  Rescaled the units of turbulent kinetic energy in the bulk mixed layer code to
[H L2 T-2 ~> m3 s-2 or J m-2] to reduce the influence of the Boussinesq
reference density in non-Boussinesq configurations, with similar changes to the
internal units of diagnostics and energy sources.  Also revised how Hmix_min is
set in bulkmixedlayer_init to avoid any dependency on the Boussinesq reference
density.

  Add a U_star argument to find_starting_TKE and use find_ustar to set it. In
some places a thickness-based definition of U_star is used. Also added logic to
the code setting the starting TKE and k_Ustar so that they supports a greater
range of valid combinations of available input variables.

  Added the option of using non-Boussinesq energetic calculations in the bulk
mixed layer code, which avoids any dependence on the Boussinesq reference
density, but do the calculations with the approximation that specific volume is
conserved during mixing, thereby ignoring certain weak thermobaric effects.  The
use of this new option is controlled by the new runtime parameter
BML_NONBOUSSINESQ, which is false by default except in fully non-Boussinesq
mode.  All of the new code is wrapped in logical branches that are selected by
the new logical variable CS%nonBous_energetics in the bulk mixed layer control
structure.  This option changes which equation of state routines are called by
the bulk mixed layer module.

  When in non-Boussinesq mode, use forces%tau_mag and tv%SpV_avg instead of
forces%ustar and GV%Rho0 to determine the surface TKE flux.

  Use SpV_avg to rescale opacity in non-Boussinesq mode via the use of an
optional argument to extract_optics_slice

  Use a call to average_specific_vol to translate the mass of the mixed layer
into the mixed layer thickness.

 As a part of these changes, there is extensive but systematic revision to the
code.  Within the bulkmixedlayer_CS type the units of 10 elements (mostly
diagnostics) are changed, and there is one new logical element. There are 17 new
arguments to internal subroutines in the bulk mixed layer module, while the
units of another 11 are changed. There are 35 new or renamed internal variables,
while the units of another 28 internal variables are changed. A total of 23
thickness conversion factors were eliminated, and the remaining references to
the Boussinesq reference density are only used in Boussinesq mode. Apart from
the new runtime parameter, the external interfaces are unchanged.

  By default, answers are bitwise identical for Boussinesq or semi-Boussinesq
configurations, but there is a new entry (BML_NONBOUSINESQ) in some
MOM_parameter_doc files, and answers do change in non-Boussinesq mode and become
independent of the value of RHO_0.

  This commit changes differential_diffuse_TS, diagnoseMLDbyDensityDifference,
diagnoseMLDbyEnergy and applyBoundaryFluxesInOut to make them appropriate for
use in non-Boussinesq mode, and to eliminate dependencies on the Boussinesq
reference density when in that mode.  It also adds a new optional argument to
extract_optics_slice to enable the use of the layer specific volumes to
translate opacities into thickness-based units.

  The specific set of changes include:

 - Add the optional argument SpV_avg to extract_optics_slice and use it along
   with an appropriate value for opacity_scale to optionally convert the units
   of opacity from [Z-1 ~> m-1] to [H-1 ~> m-1 or m2 kg-1] in non-Boussinesq
   mode without making use of the Boussinesq reference density.

 - Use thickness_to_dz and work with internal variables in vertical distances in
   the denominator of diffusive flux calculations in differential_diffuse_T_S,
   diagnoseMLDbyDensityDifference and diagnoseMLDbyEnergy.

 - Refactored diagnoseMLDbyEnergy for probable efficiencies by calling the
   equation of state with contiguous arguments.

 - Use specific volume derivatives to calculate non-Boussinesq mode buoyancy
   fluxes in calculateBuoyancy_Flux1d and applyBoudaryFluxesInOut.

 - Use the inverse of SpV_avg rather than Rho0 in the calculation of the energy
   input used to drive river mixing when in non-Boussinesq mode.

  There are now separate thickness and depth change internal variables in
several places to avoid any dependency on the Boussinesq reference density when
in non-Boussinesq mode.  A total of 8 rescaling factors were eliminated, and in
one place, GV%Rho0 was replaced with GV%H_to_RZ.

  All Boussinesq answers are bitwise identical, but non-Boussinesq answers will
change and become less dependent on the Boussinesq reference density, and there
is a new optional argument to a publicly visible subroutine.