Replaces LU_pred to L_diag, since now this logical only controls
if diagnostics should be posted.

This commit adds the latest updates to the vertFPmix subroutine
after Bill Large did some cleaning. We have highlight places in
the code where work must be done.

omega_w2x is the counter-clockwise angle of the wind stress with
respect to the horizontal abscissa (x-coordinate) at tracer
points [rad]. This variable is needed in the vertPFmix subroutine.

This line of code was lost during the last merge.

* Rename MOM_lateral_boundary_diffusion.F90 to
  MOM_hor_bnd_diffusion.F90.
* Following the suggestion from a reviewer, MOM_lateral_boundary_diffusion
  has been renamed to MOM_hor_bnd_diffusion. Many submodules related to the
  'old; lateral diffusion have been renamed throughout the code.
  LBD has been replaced to HBD.
* Tested that answers for GMOM do not change.

Build and store the HBD grid outside the tracer loop since the
same grid is used in all tracers. This makes this module
more computationaly efficient. A GMOM case run for 10 days and with
3 tracer is ~ 7.5 % faster.

add desc argument to log_param calls in MOM_CFC_cap

CFC_BC_FILE must be specified if USE_CFC_CAP=.true.

use hemispheric averages poleward of 10 degrees latitude
linearly interpolate between 10S and 10N

correct bug that atm cfc12 was used in cfc11 flux computation

CFC_cap forcing update

* Add subroutine hbd_grid_test, which mimics subroutine hbd_grid
  but it is only used in the unit tests;

* Add unit tests for hbd_grid_test and fix existing tests for
  fluxes_layer_method;

* Delete unused code and fix the format of doxygen throughout
  the module.

…OM_CFC_cap

refs moved out of nuopc cap code, MOM_forcing_type, MOM_variables

call CFC_cap_set_forcing in call_tracer_set_forcing

add call to call_tracer_set_forcing in nuopc cap

add arguments to call_tracer_set_forcing

increase width in MOM_CFC_cap unit test output

correct typo in oil_tracer

…of_cap

migrate nearly all refs to CFC_cap into MOM_tracer_flow_control and MOM_CFC_cap

Improvements in the MOM_lateral_boundary_diffusion module

* changes in nuopc cap, infra, and MOM.F90 to receive ensembe id from the coupler (alternative to FMS ensemble mngr)

* multi-instance logfile name correction in nuopc cap

* append ensemble suffix to _doc files

* changes in rpointer and restart file name handling to accommodate multi-instance CESM runs

* remove fms2_io_mod usage in FMS1/MOM_ensemble_manager_infra.F90

* rm whitespace in mom_cap

…ix_draft_18April2023

This commit adds the option to apply a linear decay in the
neutral diffusion fluxes within a transition zone defined
by the boundary layer depths of adjacent columns. This option
is controlled by a new parameter NDIFF_TAPERING, which is only
available when NDIFF_INTERIOR_ONLY=True. By default
NDIFF_TAPERING=False and answers are bitwise identical.

…0406

Merge changes from GFDL to main (2023-04-06)

Writes useful fields when the diffusivity of viscosity is less than zero.
The should help understanding the root cause of such cases and facilitate
the necessary adjustments.

…ngmuir_kpp

Output relevant fields when diff or visc < 0

Simplifies and reduces the code by adding hbd to the neutral
diffusion contril structure. This avoid the need to "extract"
hbl multiple times. Answers are bitwise indenticals.

This commit adds a vertical dimension to the tracer diffusivities
(Kh_u and Kh_v) and associated coefficiets (coef_x and coef_y).
The following diagnostics were changed from 2D (lat/lon) to 3D
(lat/lon/depth): KhTr_u, KhTr_v, and KhTr_h.

To preserve old answers, the values of all modified arrays are
depth independent by default. The option to apply the equivalent
barotropic structure as the vertical structure of the tracer
diffusivity is also introduced and this can be controlled
via a new parameter: KHTR_USE_EBT_STRUCT (default is false).

Following up on the previous commit, where a vertical dimension
was added to the tracer diffusivities, this commit modifies the
HBD module to work with this change. To do so, parameter khtr_u
(diffusivity times the time step) is calculated at cell centers
and then remapped onto the HBD vertical grid. All unit tests
in this module were updated to conform with this change.

This commit also makes the default value of HBD_DEBUG equal to
the value set for DEBUG.

This commit modifies the neutral diffusion module to work with
3D diffusivities. When the diffusivities are depth dependent
(KHTR_USE_EBT_STRUCT=True), a new array (Coef_h, with values at
tracer points and at vertical interfaces) with a four-point average
between Coef_x and Coef_y is introduced. This array is then used to
calculate zonal and meridional neutral fluxes via optional arguments
and using an existing four-point average (vertical interfaces of two
tracer cells) inside subroutine neutral_surface_flux. The same approach
is already used when tapering the neutral diffusive fluxes. In this case,
however, the unit of the output from neutral_surface_flux (Flx) is
modified because the flux of the tracer between pairs of neutral layers is
multiplied by the average of Coef_h. To avoid double counting Coef_h,
the code block for updating the tracer concentration from divergence
of neutral diffusive flux components also had to be modified
for when KHTR_USE_EBT_STRUCT=True. Similar for
diagnostics trans_x_2d and trans_y_2d.

This commit also makes the default value of NDIFF_DEBUG equal to
the value set for DEBUG.

* Document and unit test for mu(z) in MLE parameterization

- Renamed function from psi(z) to mu(sigma)
- Added comments and units in function mu(sigma)
- Added [numerical] unit tests for mu(z), including special limits,
  special values, and one test value (checked against a python script).

* Adds the Bodner et al. 2023 version of MLE

Changes:

- Allow MLE parameterization to see surface buoyancy flux return from
  PBL scheme (affects MOM.F90, MOM_variables.F90:vertvisc_type,
  MOM_diabatic_driver.F90, MOM_set_viscosity.F90)

- Adds the Bodner et al., 2023, parameterization of restratification by
  mixed-layer eddies to MOM_mixed_layer_restrat.F90

  - This is a new subroutine rather than embedded inside the previous
    "OM4" version. It uses different inputs, different parameters,
    filters the BLD differently,

- Renamed mixedlayer_restrat_general to mxiedlayer_restrat_OM4 to better
  distinguish the two versions.

- Added function rmean2ts to extend the resetting running-mean time
  filter used in OM4 to use different time scales when growing or
  decaying. While mathematically the same in the limit of a zero
  "growing" time-scale, the implementation differs in the use of a
  reciprocal instead of division so was not added to the OM4 version.

- Updated module documentation

Co-authored-by: Abigail Bodner <[email protected]>

* Add Bodner MLE testing

This patch adds the Bodner MLE testing parameters to the tc2.a test.

* +Add Pa_to_RL2_T2 and Pa_to_RLZ_T2 to US type

  Add the combined unit scaling factors Pa_to_RL2_T2 and Pa_to_RLZ_T2 to the
unit_scale_type to rescale pressures and wind stresses.  All answers are bitwise
identical, but there are two new elements in a public type.

* Use US%Pa_to_RL2_T2 to rescale pressures

  Use the new combined unit scaling factor US%Pa_to_RL2_T2 to rescale input
pressure fields and US%Pa_to_RLZ_T2 to rescale input wind stresses in various
places in the MOM6 code, including in the solo_driver and FMS_cap drivers.
Analogous changes could also be made to the mct and nuopc surface forcing files,
but have been omitted for now.  All answers are bitwise identical.

* +Add runtime parameter TAUX_MAGNITUDE

  Added the new runtime parameter TAUX_MAGNITUDE to set the strength of the
zonal wind stresses when WIND_CONFIG = "2gyre", "1gyre" or "Neverworld", with a
default that matches the previous hard-coded dimensional parameters that were
used to specify the wind stresses in these cases.  Also use US%Pa_to_RLZ_T2 to
rescale wind stresses throughout solo_driver/MOM_surface_forcing.F90.  By
default, all answers are bitwise identical, but there is a new runtime parameter
in the MOM_parameter_doc files for some test cases.

* Correct MLD_EN_VALS rescaling

  Correct inconsistent dimensional rescaling of the input values of MLD_EN_VALS,
setting them all to [R Z3 T-2 ~> J m-2] to reflect that these are energies
associated with vertical turbulent mixing.  This fixes a rescaling bug when
these energies are set to non-default values at runtime, but all answers and
output are bitwise identical when no rescaling is used.

* Add better error handling to read_var_sizes

  Add better error handling to read_var_sizes when a missing file or missing
variable is provided as an argument.  Without this change the model fails with a
segmentation fault on line 768 of MOM_io.F90 if a bad file or variable name is
provided.  With this change, a useful error message is returned.  All answers
are bitwise identical in all cases that worked previously.

* Checksum unrescaled non-Boussinesq thicknesses

  Redid the scaling of 52 checksum or check_redundant calls for thickness or
transports to use the MKS counterparts of the thickness units (i.e., m and m3/s
or kg/m2 and kg/s, depending on the Boussinesq approximation), rather than
always rescaling them to m or m3/s.  In Boussinesq mode, everything remains the
same, but in non-Boussinesq mode, this means that the model's actual variable
are being checksummed and not a version that is rescaled by division by the
(meaningless?) Boussinesq reference density.  All solutions are bitwise
identical, but some debugging output will change in non-Boussinesq mode.

* (*)Use conversion factor for masscello diagnostic

  Use a conversion factor to rescale the units of masscello, just like every
other diagnostic.  This does not change the diagnostic itself, but it changes
the order of the rescaling and the vertical remapping of this diagnostic onto
other coordinates (like z) or spatial averaging of this diagnostic, which can
change values in the last bits for this diagnostic for Boussinesq models (but
not for non-Boussinesq models, for which the conversion factor is an integer
power of 2).  As a result some of the diagnostics derived from masscello can
differ and this commit nominally fails the TC testing for reproducibility across
code versions.  All solutions and primary diagnostics, however, are bitwise
identical, and even the derived diagnostic calculations are mathematically
equivalent.

* +Remove rescaling factors from restart files

  Remove the code to account for unit rescaling within the restart files.  This
rescaling within the restart files has not been used in the code since March,
2022, and the model will work with older restart files provided that they did
not use dimensional rescaling, and even if they did they can be converted not to
use rescaling with a short run with the older code that created them.  Also
removed the publicly visible routines fix_restart_scaling and eliminated the
m_to_H_restart element of the verticalGrid_type; in any cases of non-standard
code using this element, it should be replaced with 1.0.  The various
US%..._restart elements and fix_restart_unit_scaling are being retained for now
because they are still being used in the SIS2 code.

  These changes significantly simplify the code, and they lead to a handful
of constants that are always 1 not being included in the MOM6 restart files.
All answers are bitwise identical, but a publicly visible interface has been
eliminated, as has been an element (GV%m_to_H_restart) of a transparent type.

* +Add MOM_EOS_Wright_Full

  Added the new module MOM_EOS_Wright_full to enable the use of the version of
the Wright equation of state that has been fit over the larger range of
temperatures (-2 degC to 40 degC), salinities (0 psu to 40 psu) and pressures (0
dbar to 10000 dbar), than the does the restricted range fit in MOM_EOS_Wright,
which had been fit over the range of (-2 degC to 30 degC), (28 psu to 38 psu)
and (0 to 5000 dbar).  Comments have been added to both modules to clearly
document the range of properties over which they have been fitted.  The new
equation of state is enabled by setting EQN_OF_STATE = "WRIGHT_FULL".  In
addition, the default values for TFREEZE_FORM and EOS_QUADRATURE were changed
depending on the equation of state to avoid having defaults that lead to fatal
errors.  All answers are bitwise identical in any cases that currently work, but
there are new entries in the MOM_parameter_doc files.

  For now, only the coefficients have been changed between MOM_EOS_Wright and
MOM_EOS_Wright_full, but this means that it does not yet have all of the
parentheses that it should, as github.com/mom-ocean/MOM6/issues/1331 discusses.
A follow up PR should add appropriate self-consistency and reference value
checks (with a tolerance) for the various EOS routines, and then add enough
parentheses to specify the order of arithmetic and hopefully enhance the
accuracy.  Ideally this can be done with the new equation of state before it
starts to be widely used, so that we can avoid needing a extra code to reproduce
the older answers.

* Fix and tidy Wright_EOS API documentation

  Cleaned up the comments describing the routines and added a proper doxygen
namespace block at the end of the MOM_EOS_Wright and MOM_EOS_Wright_full
modules, based on changes that A. Adcroft had on a detached branch of MOM6.
Only comments are changed, and all answers are bitwise identical.

* (*)Rearranged parentheses in MOM_EOS_Wright_full

  Added parentheses to all expressions with three or more additions or
multiplications in the MOM_EOS_Wright_full code, so that different compilers
and compiler settings will reproduce the same answers in more cases.  In doing
this, an effort was made to add the smallest terms first to reduce the impact
of roundoff.  In some cases, the code was deliberately rearranged to cancel
out the leading order terms more completely.  In addition, two bugs had been
identified in calculate_density_second_derivs_wright_full.  These were corrected
and the entire routine substantially refactored with renamed variables to make
the derivation easier to follow and verify.  Apart from the bug corrections in
the calculation of drho_dt_dt and drho_dt_dp, the changes in the expressions
are mathematically equivalent, but they might make the model less noisy in some
cases by reducing contributions from round-off errors.

  Also added comments highlighting two bugs in the drho_dt_dt and drho_dt_dp
calculations in calculate_density_second_derivs_wright in the original
MOM_EOS_Wright code, but did not correct them to preserve the previous answers.

* +Created the new module MOM_EOS_Wright_red

  Created a new module, MOM_EOS_Wright_red, that uses the reduced range fit
coefficients from the Wright EOS paper, but uses the parentheses, expressions
and bug fixes that are now in MOM_EOS_Wright_full.  To use this new module, set
EQN_OF_STATE="WRIGHT_RED". This new form is mathematically equivalent using
EQN_OF_STATE="WRIGHT" (apart from correcting the bugs in the calculations of
drho_dt_dt and drho_dt_dp), but the order of arithmetic is different, so the
answers will differ.  This change is probably as close as we can come to
addressing the issues discussed at github.com/mom-ocean/MOM6/issues/1331, so
that issue should be closed once this commit is merged onto the main branch.
Also corrected some misleading error messages in MOM_EOS and modified the code
to properly handle the case for equations of state (like NEMO and UNESCO) that
do not have a scalar form of calculate_density_derivs, but do have an array
form.  By default, all answers are bitwise identical.

* *Fix bug in calculate_spec_vol_linear with spv_ref

  Corrected a sign error in calculate_spec_vol_array_linear and
calculate_spec_vol_scalar_linear when a reference specific volume is provided.
This bug will cause any configurations with EQN_OF_STATE="LINEAR" and
BOUSSINESQ=False (neither of which is the default value) to have the wrong sign
of the pressure gradients and other serious problems, like implausible sea
surface and internal interface heights.  This combination of parameters would
never be used in a realistic ocean model.  There are no impacted cases in any of
the MOM6-examples tests cases, nor those used in the ESMG or dev/NCAR test
suites, and it is very unlikely that any such case would work at all.  This bug
was present in the original version of the calculate_spec_vol_linear routines,
but it was only discovered after the implementation of the comprehensive
equation of state unit testing.  This will change answers in configurations that
could not have worked as viable ocean models, but answers are not impacted in
any known configuration, and all solutions in test cases are bitwise identical.

* +Add EOS_unit_tests

  Added the new publicly visible function EOS_unit_tests, along with a call to
it from inside of unit_tests.  These tests evaluate check values for density and
assess the consistency of expressions for variables that can be derived from
density with finite-difference estimates of the same variables.  These tests
reveal inconsistencies or omissions with several of the options for the equation
of state.  The EOS self-consistency tests that are failing are commented out for
now, so that this redacted unit test passes.  All answers are bitwise identical,
but there can be new diagnostic messages written out.

* Fix doxygen labels in EOS_Wright_full and _red

  Changed recently added doxygen labels in the two newly added EOS_Wright_red and
EOS_Wright_full modules to avoid reusing names that were already being used
by EOS_Wright.  All answers are bitwise identical, but the doxygen testing that
had been failing for the previous 5 commits is working again.

* *+NEMO equation of state self-consistency

  Corrected numerous issues with the NEMO equation of state so that it is now
self consistent:

- Modified how coefficients are set in MOM_EOS_NEMO so that they are guaranteed
to be internally self-consistent, as verified by the EOS unit tests confirming
that the first derivatives of density with temperature and salinity are now
consistent with the equation of state.  Previously these had only been
consistent to about 7 decimal places, and hence the EOS unit tests were failing
for the NEMO equation of state.

 - Added new public interfaces to calculate_density_second_derivs_NEMO, which
had previously been missing.

 - Added code for calculate_compress_nemo that is explicitly derived from the
NEMO EOS.  The previous version of calculate_compress_nemo  had worked only
approximately via a call to the gsw package

  With these changes, the NEMO EOS routines are now passing the consistency
testing in the EOS unit tests.  Answers will change for configurations that use
the NEMO EOS to calculate any derivatives, and there are new public interfaces,
but it does not appear that the NEMO equation of state is in use yet, at least
it is not being used at EMC, FSU, GFDL, NASA GSFC, NCAR or in the ESMG
configurations.

This commit addresses the issue raised at github.com/mom-ocean/MOM6/issues/405.

* +Add calculate_density_second_derivs_UNESCO

  Added the new public interface calculate_density_second_derivs_UNESCO, which
is an overload for both scalar and array versions, to calculate the second
derivatives of density with various combinations of temperature, salinity and
pressure.  Also added a doxygen block at the end of MOM_EOS_UNESCO.F90 to
describe this module and the papers it draws upon.  Also replaced fatal
errors in MOM_EOS with calls to these new routines.  All answers are bitwise
identical, but there are newly permitted combinations of options that previously
failed.

* (*)+Added calc_density_second_derivs_wright_buggy

  Added the new public interface calc_density_second_derivs_wright_buggy to
reproduce the existing answers and corrected bugs in the calculation of the
second derivatives of density with temperature and with temperature and pressure
in in calculate_density_second_derivs_wright.  Also added the new runtime
parameter USE_WRIGHT_2ND_DERIV_BUG to indicate that the older (buggy) version of
calculate_density_second_derivs_wright is to be used.  Most configurations will
not be impacted, but by default answers will change with configurations that use
the Wright equation of state and one of the Stanley or similar nonlinear EOS
parameterizations, unless USE_WRIGHT_2ND_DERIV_BUG is explicitly set to True.

  This commit also activates the self-consistency unit testing with the Wright
equation of state (now that it passes) and limited unit testing of the TEOS-10
equation of state, omitting the second derivative calculations, one of which is
failing (the second derivative of density with salinity and pressure) due to a
bug in the TEOS10/gsw code.  Also added a unit test for consistency of the
density and specific volume when an offset reference value is used.

* *Refactor MOM_EOS_UNESCO.F90

  Refactored the expressions in MOM_EOS_UNESCO.F90, adding parentheses to
specify the order of arithmetic, starting with the highest-order terms first for
less sensitivity to round-off.  Also added comments to better describe the
references for these algorithms.  Although the revised expressions are all
mathematically equivalent, this commit will change answers for any cases that
use EQN_OF_STATE = "UNESCO".  However, it is believed based on a survey of the
MOM6 community that there are no active configurations that use this equation of
state.

* *Refactor MOM_EOS_NEMO.F90

  Refactored the expressions in MOM_EOS_NEMO.F90, adding parentheses to specify
the order of arithmetic, starting with the highest-order terms first for less
sensitivity to round-off.  A number of internal variables were also renamed for
greater clarity, and a number of comments were revised to better describe the
references for these algorithms..  Although the revised expressions are all
mathematically equivalent, this commit will change answers for any cases that
use EQN_OF_STATE = "NEMO".  However, there is another recent commit to this file
that also changes answers (specifically the density derivatives) with this
equation of state, and it is believed based on a survey of the MOM6 community
that there are no active configurations that use this equation of state.

* +Add MOM_EOS_Roquet_SpV.F90

  Added the new equation of state module MOM_EOS_Roquet_SpV with the polynomial
specific volume fit equation of state from Roquet et al. (2015).  This equation
of state has also been added to MOM_EOS, where it is enabled by setting
EQN_OF_STATE="ROQUET_SPV".  Two other new valid settings have been added to
EQN_OF_STATE, "ROQUET_RHO" and "JACKETT_MCD", which synonymous with "NEMO" and
"UNESCO" respectively, but more accurately reflect the publications that
describe these fits to the equation of state.  The EoS unit tests are being
called for the new equation of state (it passes).  By default, all answers are
bitwise identical, but there are numerous new publicly visible interfaces.

* +Add MOM_EOS_Jackett06.F90

  Added the new equation of state module MOM_EOS_Jackett06 with the rational
function equation of state from Jackett et al. (2006).  This uses potential
temperature and practical salinity as state variables, but with a fit to more
up-to-date observational data than Wright (1997) or UNESCO / Jackett and
McDougall (1995).  This equation of state has also been added to MOM_EOS, where
it is enabled by setting EQN_OF_STATE="JACKETT_06".  The EoS unit tests are
being called for the new equation of state (it passes).  This commit also adds
slightly more output from successful EoS unit tests when run with typical levels
of verbosity.  By default, all answers are bitwise identical, but there are
numerous new publicly visible interfaces.

* *+Add calculate_specvol_derivs_UNESCO

  Added the routine calculate_specvol_derivs_UNESCO to calculate the derivatives
of specific volume with temperature and salinity to the MOM_EOS_UNESCO module.
Also added some missing parentheses elsewhere in this module so that the answers
will be invariant to complier version and optimization levels.  Also revised the
internal nomenclature of the parameters in this module to follow the conventions
of the other EOS modules.  Although the revised expressions are mathematically
equivalent, this commit will change answers for any cases that use EQN_OF_STATE
= "UNESCO".  However, it is believed based on a survey of the MOM6 community
that there are no active configurations that use this equation of state.  There
is a new publicly visible routine.

* +Add EOS_fit_range and analogs for each EoS

  Added the new publicly visible subroutine EOS_fit_range and equivalent
routines for each of the specific equation of state modules to return the range
of temperatures, salinities, and pressures over which the observed data have
been fitted.  This is also tested for in test_EOS_consistency to indicate
whether a test value is outside of the fit range, but the real purpose will be
to flag and then figure out how to deal with the case when the ocean model is
called with properties for which the equation of state is not valid.  Note that
as with all polynomial or other functional fits, extrapolating far outside of
the fit range is likely to lead to bad values, but things may not be so bad for
values that are only slightly outside of this range.  However the question of
how far out of the fit range these EoS expressions become inappropriate for each
of temperature, salinity and pressure is as yet unresolved.  All answers and
output are bitwise identical, but there are 10 new public interfaces.

* Do not include MOM_memory.h in EoS modules

  Removed unused and unnecessary #include <MOM_memory.h> statements from 5
equation of state modules.  All answers are bitwise identical.

* *Refactor calculate_specific_vol_wright_full

  Refactored the specific volume calculations for the WRIGHT_FULL and WRIGHT_RED
equations of states for simplicity or to reduce the impacts of roundoff when
removing a reference value.  Also added code to multiply by the reciprocal of
the denominator rather than dividing in several places in the int_spec_vol_dp
routines for these same two equations of state, both for efficiency and greater
consistency across optimization levels.  These changes are mathematically
equivalent but will change answers at roundoff with these two equations of state,
but they are so new that they can not have been used yet.

* +Renamed MOM_EOS_NEMO to MOM_EOS_Roquet_rho

  Renamed the module MOM_EOS_NEMO to MOM_EOS_Roquet_rho to more accurately
reflect its provenance, although setting either EQN_OF_STATE = NEMO or
EQN_OF_STATE = ROQUET_RHO will still work for using this code.  All answers
are bitwise identical, and previous input files will still work, but there are
some minor changes in the MOM_parameter_doc files.

* *Avoid re-rescaling T and p in MOM_EOS_Roquet_rho

  Refactored MOM_EOS_Roquet_rho and MOM_EOS_Roquet_SpV to work directly with
conservative temperatures in [degC] and pressures in [Pa] rather than
normalizing them as in the original Roquet publication.  However, the
coefficients are still set using the values directly copied from that paper, but
rescaled where they are declared as parameters, enabling (or requiring)
compilers to precalculate them during compilation.  These changes are
mathematically equivalent but will change answers at roundoff with these two
equations of state, but they are not believed to be in use yet.

* +Add calculate_TFreeze_TEOS_poly

  Added the overloaded interface calculate_TFreeze_TEOS_poly to MOM_TFreeze to
use the 23-term polynomial expression from TEOS-10 for the freezing point in
conservative temperature as a function of pressure and absolute salinity.  This
gives results that agrees to within about 5e-4 degC with the algorithm used by
calculate_TFreeze_TEOS10, which calls the gsw TEOS10 code that does an iterative
inversion of a balance of chemical potentials to find the freezing point (see
the TEOS10 documentation for more details).  Also added testing for the freezing
point calculations to the EOS_unit tests via the new internal subroutine
test_TFr_consistency.  This new freezing point calculation is invoked by setting
TFREEZE_FORM = TEOS_POLY.  By default, all answers are bitwise identical, but
there are some minor changes in the comments in some MOM_parameter_doc files,
and there are several new interfaces.

* +*Add MOM_temperature_convert.F90

  Added the new module MOM_temperature_convert, which contains the elemental
functions poTemp_to_consTemp and consTemp_to_poTemp to convert potential
temperature to conservative temperature and the reverse.  These routines are
mathematically equivalent to the TEOS-10 functions gsw_ct_from_pt and
gsw_pt_from_ct, but with some refactoring and added parentheses to help ensure
identical answers across compilers or levels of optimization.  Also added the
new subroutines pot_temp_to_cons_temp and prac_saln_to_abs_saln, and added the
new optional argument use_TEOS to convert_temp_salt_for_TEOS10, and
cons_temp_to_pot_temp and abs_saln_to_prac_saln.

  The equivalency between the new code and their gsw_ counterparts is
demonstrated in new tests in the new function test_TS_conversion_consistency,
which in turn is called from EOS_unit_tests.

  All answers are mathematically equivalent, but because of the choice to use
the new code by default there could be changes at the level of roundoff in some
cases that use conservative temperature as their state variable but initialize
it from potential temperature.  There are not any such cases yet in the
MOM6-examples test suite, nor are there believed to be any such MOM6
configurations that are widely used.  This commit introduces a new module and
several new functions or subroutines with public interfaces.

* Update _Equation_of_State.dox

  Updated _Equation_of_State.dox to reflect the new options for the equation of
state and freezing point calculations.

* +Eliminate use_TEOS arg to cons_temp_to_pot_temp

  Eliminate use_TEOS optional arguments that were recently added to
cons_temp_to_pot_temp and 4 other thermodynamic variable conversion functions,
along with calls to gsw_pt_to_ct and similar conversion functions.  All answers
in the MOM6-examples test suite are bitwise identical.

* +Make calculate_density_array private

  Removed calculate_density_array from the overloaded public calculate_density
interface, and similarly for the other EOS calculate_..._array routines, to help
standardize how they are called.  Calculate_density_derivs_array  is the one
exception is because it is being called from SIS2 and has to stay publicly
visible for now.  Additionally, the scalar and 1-d versions of the
calculate_stanley_density routines were refactored to just use calculate_density
and calculate_density_second_derivs call and avoid any EoS-specific logic, while
the unused routine calculate_stanley_density_array is eliminated altogether.
All answers are bitwise identical, including in extra tests that use the
stanley_density routines.

* +Rename WRIGHT_RED to WRIGHT_REDUCED

  Revised the setting EQN_OF_STATE to select the Wright equation of state with
the reduced-range fit to "WRIGHT_REDUCED" (instead of "WRIGHT_RED") for greater
clarity, in response to a comment in the review of the pull request with this
sequence of code revisions.  All answers are bitwise identical, but this changes
the text for a recently added input parameter and it leads to changes in some
comments in the MOM_parameter_doc files.

* Removal of FMS1 I/O from FMS2 I/O infra

This patch removes the calls to FMS1 I/O (fms_io_mod, mpp_io_mod) from
the FMS2 infra layer, and now exclusively uses FMS2 for those
operations.

FMS2 I/O is currently restricted to files which use domains; files which
do not use them are delegated to the native netCDF layer.  The reasoning
for this is that FMS is required to define the formatting of
domain-decomposed I/O; for single-file I/O, this is not necessary.

This does not remove all references to FMS1 I/O from MOM6, only those in
the I/O layer.

Several minor changes are included to accommodate the change:

* MOM restart I/O now always reports its MOM domain.  Previously, the
  domian was omitted when PARALLEL_RESTARTFILES was false, in order to
  trick FMS into handling this as a single file.  We now generate a new
  domain with an IO layout of [1,1] when single-file restarts are
  requested.

* The interface acceleration (g') was incorrectly set to the layer grid
  (Nk) rather than the interface grid (Nk+1).  This did not appear to
  change any answers, but when Vertical_coordinate.nc was moved to the
  netCDF layer, it detected this error.  This is fixed in this patch.

* Remove FMS1 calls from MOM_domains_infra

* Add .nc extension to ALE Vertical_coordinate.

The `Vertical_coordinate.nc` files has two points of creation,
MOM_coord_initialization and MOM_ALE.  Having moved the file from the
infra to netCDF I/O layer, the .nc extension is no longer automatically
applied.

The extension was explicitly added to `Vertical_coordinate` in
MOM_coord_initialization, but not to MOM_ALE.  This patch adds the
extension.

Thanks to Kate Hedstrom for detecting this and Keith Lindsay for the
proposed fix.

* +Remove optional argument eta_to_m from find_eta

  Eliminate the unused optional argument eta_to_m from the two find_eta routines
for simplicity and code clarity.  These were used during the transition of the
units of the interface height variables, but they are now using [Z ~> m] units
everywhere, with the unscaling occurring via conversion factors in the
register_diag calls.  All answers are bitwise identical, but there is al
optional argument that is removed from a public interface.

* +Initialize thicknesses in height units

  Pass arguments in height units rather than thickness units to most of the
routines that initialize thickness or temperatures and salinities.  These
routines are already undoing this scaling and working in height units, and it is
not possible to convert thicknesses to thickness units in non-Boussinesq mode
until the temperatures and salinities are also known.  The routines whose
argument units are altered include:

- initialize_thickness_uniform
- initialize_thickness_list
- DOME_initialize_thickness
- ISOMIP_initialize_thickness
- benchmark_initialize_thickness
- Neverworld_initialize_thickness
- circle_obcs_initialize_thickness
- lock_exchange_initialize_thickness
- external_gwave_initialize_thickness
- DOME2d_initialize_thickness
- adjustment_initialize_thickness
- sloshing_initialize_thickness
- seamount_initialize_thickness
- dumbbell_initialize_thickness
- soliton_initialize_thickness
- Phillips_initialize_thickness
- Rossby_front_initialize_thickness
- user_initialize_thickness
- DOME2d_initialize_temperature_salinity
- ISOMIP_initialize_temperature_salinity
- adjustment_initialize_temperature_salinity
- baroclinic_zone_init_temperature_salinity
- sloshing_initialize_temperature_salinity
- seamount_initialize_temperature_salinity
- dumbbell_initialize_temperature_salinity
- Rossby_front_initialize_temperature_salinity
- SCM_CVMix_tests_TS_init
- dense_water_initialize_TS
- adjustEtaToFitBathymetry

Similar changes were made internally to MOM_temp_salt_initialize_from_Z to defer
the transition to working in thickness units, although the appropriate call to
convert_thickness does still occur within MOM_temp_salt_initialize_from_Z and
the units of its arguments are not changed.

  The routine convert thickness was modified to work with a new input depth
space input thickness argument and return a thickness in thickness units, and it
is now being called after all of the routines to initialize thicknesses and
temperatures and salinities, except in the few cases where the thickness are
being specified directly in mass-based thickness units, as might happen when
they are read from an input file.

  The new option "mass_file" is now a recognized option for the THICKNESS_CONFIG
runtime parameter, and this information is passed in the new mass_file argument
to initialize_thickness_from_file.  The description of the runtime parameter
THICKNESS_IC_RESCALE was updated to reflect this change.

  The unused thickness (h) argument to soliton_initialize_velocity was
eliminated.

  The unused thickness (h) argument to determine_temperature was eliminated, as
was the unused optional h_massless argument to the same function.

  This commit also rearranges the calls to do adjustments to the thicknesses to
account for the presence of an ice shelf or to iteratively apply the ALE
remapping to occur before the velocities are initialized, so that there is a
clearer separation of the phases of the initialization.

  Also added optional height_units argument to ALE_initThicknessToCoord to
specify that the coordinate are to be returned in height_units.  If it is
omitted or false, the previous thickness units are returned, but when called
from MOM_initialize_state the new argument is being used.

  The runtime parameter CONVERT_THICKNESS_UNITS is no longer meaningful, so it
has been obsoleted.

  All answers are bitwise identical, but there are multiple changes to the
arguments to publicly visible subroutines or their units, and there are changes
to the contents of the MOM_parameter_doc files.

* +Add the new overloaded interface dz_to_thickness

  Renamed convert_thickness from MOM_state_initialization to dz_to_thickness_tv
in MOM_density_integrals, so that it can be called from other lower-level
modules. This new version also takes the tv%p_surf field into account and it has
an optional halo_size argument, analogous to that in the other routines in the
MOM_density_integrals module.  The dz_to_thickness interface is overloaded so
that it can also be used directly with temperature, salinity, and the equation
of state type if the thermo_var_ptrs is not available.  There is also a new and
separate variant of this routine, dz_to_thickness_simple, that can be used in
pure layered mode when temperature and salinity are not state variables, or
(more dangerously) if it is not clear whether or not there is an equation of
state.  This simpler version is being kept separate from the main overloaded
interface because its use may need to be revisited later in some cases.  All
answers are bitwise identical, but there are two new public interfaces,
dz_to_thickness and dz_to_thickness_simple.

* (*)Improve non-Boussinesq initialization

  This commit includes three distinct sets of changes inside of
MOM_state_initialization.F90 to better handle the initialization of
non-Boussinesq models, none of which change any answers in Boussinesq models.
These include:

 - Refactored trim_for_ice to have a separate, simpler form appropriate for use
   in non-Boussinesq mode.  The units of the min_thickness argument to
   cut_off_column top were also changed to thickness units.

 - Initialize_sponges_file was refactored to work in depth-space variables
   before using dz_to_thickness to convert to thicknesses, but also to properly
   handle the case where the input file has a different number of vertical
   layers than the model is using, in which case the previous version could have
   had a segmentation fault.

 - Code in MOM_temp_salt_initialize_from_Z was reordered to more clearly group
   it into distinct phases.  It also uses the new dz_to_thickness routine to
   convert input depths into thicknesses.

  All answers are bitwise identical in all Boussinesq test cases and all test
cases in the MOM6-examples regression suite, but answers could
be changed and improved in some non-Boussinesq cases.

* (*)Use dz_to_thickness in 4 user modules

  Use dz_to_thickness to convert vertical distances to layer thicknesses in the
sponge initialization routines in the DOME2d_initialization,
ISOMIP_initialization, dumbbell_initialization and dense_water_initialization
modules, and also in MOM_initialize_tracer_from_Z.  For the user modules,
the presence or absence of an equation of state is known and handled properly,
but MOM_initialize_tracer_from_Z works with the generic tracer code and it
it outside of the scope of MOM6 code to provide any information about the
equation of state or the state variables that would be needed to initialize
a non-Boussinesq model properly from a depth-space input file.  For now we
are doing the best we can, but this should be revisited.  All examples in
existing test cases are bitwise identical, but answers could change (and be
improved) in any non-Boussinesq variants of the relevant test cases.

* Update the Gitlab .testing modules for c5

In preparation for the migration to C5, this patch updates the modules
required to run the .testing suite.

* POSIX: generic wrappers for all setjmp.h symbols

This patch extends the generic wrappers of sigsetjmp to all of the *jmp
wrapper functions in <setjmp.h>

The C standard allows these to be defined as macros, rather than
explicit functions, which cannot be referenced by Fortran C bindings, so
we cannot assume that these functions exist, even when using a compliant
libc.

As with sigsetjmp, these functions are now disabled on default, and
raise a runtime error if called by the program.  Realistically, they
will only be defined by an autoconf-configured build.

This is required for older Linux distributions where libc does not
define longjmp.

* Autoconf: External FMS build configuration

This patch modifies the `ac/deps` Makefile used to build the FMS
depedency.  The autoconf compilation is now done entirely outside of the
`ac/deps/fms/src` directory.  This keeps the FMS checkout unchanged and
allows us to better track any development changes in that library during
development.

The .testing/Makefile was also modified to use existing rules in
deps/Makefile rather than duplicating them.

Dependency of the m4 directory is also now more explicit (albeit still
somewhat incomplete).

* Autoconf: Explicit MOM_memory.h configuration

MOM6 requires an explicit MOM_memory.h header to define its numerical
field memory layout.  Previously, autoconf provided a flag to configure
this with `--enable-*`, but was prone to two issues:

* The binary choice of symmetric/nonsymmetric prevented use of static
  headers.

* It was an incorrect use of `--enable-*`, which is intended to enable
  additional internal features; it is not used to select a mode.

To address these issues, we drop the flag and replace it with an
AC_ARG_VAR variable, MOM_MEMORY, which is a path to the file.  This
variable will default to dynamic symmetric mode,

    config_src/memory/dynamic_symmetric/MOM_memory.h

so there should be no change for existing users.

To the best of my knowledge, no one used the `--enable-*` flag, nor was
it used in any automated systems (outside of .testing), so there should
be no issue with dropping it.

.testing/Makefile was updated to use MOM_MEMORY.

* Profiling: subparameter parser support

The very crude MOM_input parser in the automatic profiler did not
support subparameters (e.g. MLE% ... %MLE), which caused an error when
trying to read the FMS clock output.

This patch adds the support, or at least enough support to avoid errors.

* +*Redefine GV%Angstrom_H in non-Boussinesq mode

  Redefined GV%Angstrom_H in non-Boussinesq mode so that it is equal to
GV%H_to_Z*GV%Angstrom_Z, just as it is in Boussinesq mode.  This will change
answers (slightly) in all cases with BOUSSINESQ = False.  In addition, this
commit adds the elements semi_Boussinesq, dZ_subroundoff, m2_s_to_HZ_T,
HZ_T_to_m2_s and HZ_T_to_MKS to the verticalGrid_type.   The first 3 new
elements are used in rescaling vertical viscosities and diffusivities.  The last
two elements are set using the new runtime parameters SEMI_BOUSSINESQ and
RHO_KV_CONVERT, which are only used or logged when BOUSSINESQ = False.  All
answers and output are identical in Boussinesq cases, but answers change and
there are new runtime parameters in non-Boussinesq cases.

* +Set_interp_answer_date and REGRIDDING_ANSWER_DATE

  Add the ability to set the answer date for the regridding code, including the
addition of the new subroutine set_interp_answer_date and the new runtime
parameter REGRIDDING_ANSWER_DATE to specify the code vintage to use with state-
dependent vertical coordinates.  There is also new optional argument to
set_regrid_params.  By default, all answers are bitwise identical, but there are
new or modified public interfaces and there is a new entry in some
MOM_parameter_doc files.

* *+Revise non-Boussinesq find_coupling_coef calcs

  Restructure one of the find_coupling_coef calculations to draw out the
stress-magnitude terms, in preparation for future steps to reduce the dependency
on the Boussinesq reference density.  Using a value of VERT_FRICTION_ANSWER_DATE
that is below 20230601 recovers the previous answers with non-Boussinesq test
cases, but this is irrelevant for Boussinesq test cases.  This updated code is
mathematically equivalent to the previous expressions but it does change answers
at roundoff in non-Boussinesq cases for recent answer dates.  There are
modifications to some comments in MOM_parameter_doc files.

* +Code to calculate layer averaged specific volumes

  Add routines to calculate and store the layer-averaged specific volume, along
with code to do the unit testing of this new capability. The new public
interfaces include avg_specific_vol, average_specific_vol, avg_spec_vol_Wright,
avg_spec_vol_Wright_full, avg_spec_vol_Wright_red and avg_spec_vol_linear.
There is also a new optional argument to test_EOS_consistency to control whether
these new capabilties are tested for a particular equation of state.  All
answers are bitwise identical, and the new capabilities pass the unit testing
for self consistency.

* +Add thickness_to_dz and calc_derived_thermo

  Added the new overloaded interface thickness_to_dz to convert the layer
thicknesses in thickness units [H ~> m or kg m-2] into vertical distances in
[Z ~> m], with variants that set full 3-d arrays or an i-/k- slice.  Also added
a field (SpV_avg) for the layer-averaged specific volume to the thermo_vars_ptr
type and the new subroutine calc_derived_thermo to set it.  This new subroutine
is being called after halo updates to the temperatures and salinities.  The new
runtime parameter SEMI_BOUSSINESQ was added to determine whether tv%SpV_avg is
allocated and used; it is stored in GV%semi_Boussinesq.  Also added the new
element GV%dZ_subroundoff to the verticalGrid_type as a counterpart to
GV%H_subroundoff but in height units.  All answers are bitwise identical, but
there is a new runtime parameter in some MOM_parameter_doc files, new elements
in a transparent type and a new public interface.

* wave structure computation into wave_speeds

wave_speeds now computes the wave structures
(eigenvectors) for each mode speed (eigenvalue)
similarly to the wave_speed (singular) function.
This is a replacement for the MOM_wave_structure
function, which could be removed in a subsequent PR.

Additional arrays for mode strucures and integral
quantities are passed as output hence this is a
breaking change for the call to wave_speeds.
However it is only called once in diabatic_driver
and is used exclusively for internal tides ray tracing.

The dimensional solutions for the wave structures are now
computed inside MOM_internal_tides, and new diagnostics
are added.

An out-of-bounds bug is also corrected for the computation
of an averaged coriolis parameter.

* remove wave_structure broken code

* Autoconf: Better Unicode Python support in makedep

The `open()` commands in `makedep` for reading Fortran source now
includes an `errors=` argument for catching bytes outside of the file
character set.  Unknown characters are replaced with the "unknown"
character (usually �) rather than raising an error.

This avoids problems with Unicode characters and older Pythons which do not
support them, as well as characters from legacy encodings which can
cause errors in Unicode.

Substitution does not break any behavior, since Unicode is only
permitted inside of comment blocks and strings.

This fixes several errors which were silent in `.testing` but were
observed by some users which using autoconf to build their own
executables.

* Autoconf: Fix Python test and allow configuration

The AC_PATH_PROGS macros used in Python testing were incorrectly using
AC_MSG_ERROR in places where a missing value for PYTHON should be if the
executable was not found.

It also did not permit for a configurable PYTHON variable, since the
autodetect was always run, even if PYTHON were set.

This has been updated so that Python autodetection only runs if
PYTHON is unset.  It also correctly reports a failed configuration if
PYTHON is not found.

(It does not, however, test of PYTHON is actually a Python interpreter,
but we can deal with that at a later date.)

* Fix PGI runtime issue with class(*)

- Some tests such as global_ALE_z crash under PGI (ncrc4.pgi20 or
ncrc5.pgi227) with
FATAL from PE    27: unsupported attribute type: get_variable_attribute_0d: file:INPUT/tideamp.nc- variable:GRID_X_T attribute: axis
- PGI in general has issues with class(*) construct and in this case
cannot recognize the axis argument to be a string.
- This mod helps PGI recognize that the argument is a string.

* Use fileset rather than threading for decompositon

MOM IO was using the `threading` flag rather than `fileset` to determine
whether a file should be forced as single file rather than
domain-decomposed.  This patch applies the correct flag.

* FMS2 interpolation ID replaced with derived type

All instances of an FMS ID to the internal interpolation content is
replaced with a derived type containing additional metadata recording
the field's origin filename and fieldname.

This additional information is required in order to replicate the axis
data from the field, which is no longer provided by FMS2.

The abstraction of this type also allows us to either extend it or
redefine it in other frameworks as needed in the future.

This primarily affects the usage of the following functions:

- init_external_field
- time_interp_external
- horiz_interp_and_extrap_tracer

The following solvers are updated:

- MOM_open_boundary
- MOM_ice_shelf
- MOM_oda_driver
- MOM_MEKE
- MOM_ALE_sponge
- MOM_diabatic_aux

Of these, OBC was the most significant.  The integer handle (fid) was
previously used to determine if each segment field was constant or (if
negative) read from a file.  After being replaced by the derived type, a
new flag was added to make this determination.

All of the coupled drivers have been modified, since they support time
interpolation of T and S fields.

- FMS
- MCT
- NUOPC

The NUOPC driver also includes modifications to its CFC11 and CFC12
fields.  Changes to the MOM CFC modules replaces an `id == -1`-like
test, which is not used by the derived type.  This check has been
removed, and we now solely rely on the `present(cfc_handle)` test.

While this could change behavior, there does not seem to be any scenario
where init_external_field would return -1 but would be passed to the
function.  (But I may eat these words.)

* FMS2: Remove MPP-based axis data access

With removal of axis-based operations in FMS2 I/O, this patch removes
references to these calls and replaces them with MOM `axes_info` types.
References to FMS1 read into an `axistype`, but the contents are
transferred to an `axis_info`.  FMS2 directly populates the `axis_info`
content.

The `get_external_field_info` calls are modified to return `axis_info`
rather than `axistype`.

The redundant `get_axis_data` function is also removed from
`MOM_interp_infra`, since `get_axis_info` provides an equivalent
operation.

Generally speaking, this is not an improvement of the codebase.  The
FMS1 layer does a redundant copy of data from `axistype` to `axis_info`.
The FMS2 layer is significantly worse, and re-opens the file to read the
axis data for each field!  But if the intention is to leverage the
existing API, then I don't think we have any choice at the moment.

Assuming this is a relatively infrequent operation, this should not
cause any measureable issues, but it needs to be watched carefully.

* FMS2: Update time_interp_external functions

This patch shifts all remaining time_interp_external functions from
time_interp_external to equivalent ones in time_interp_external2.

Internally, time-interpolated fields are initialized with `ongrid` set
to `.true.`, and such fields are assumed to be on-grid.

This seems to hold for all existing instances of `time_interp_external`,
but needs to be monitored in the future somehow.

* FMS2: Case-insensitive init_external_field

The FMS1 implementation of init_external_field is case-insensitive, but
the FMS2 implementation is case-sensitive, which can cause errors in
older established input files.

This patch sweeps through the fields of the input files and checks for a
case-insensitive match (using lowercase()).  This requires an additional
open/close of the file.

* Implementation of ZB sheme

* Filters for ZB. Regression changed (FGR changed to amplitude)

* Rotate test is passed. Regression changed (order of operatrions)

* ZB submitted via PR

* ZB: Response to the code review

* Update icebergs source path in nolibs build

The icebergs project now includes drivers and tests which can interfere
with the coupled nolibs build, so we only pass its src directory to
mkmf.

* +Make units argument mandatory for get_param_real

  This commit includes changes to the get_param_real and log_param_real
interfaces to make the units arguments mandatory.  It also adds an optional
unscale argument to the log_param_real interfaces.  Without other changes in the
previous commits, it will cause the MOM6 code to fail to compile.  However, by
itself this commit does not change any answers or output.

* github workflows: update to use actions/checkout@v3

- Update actions/checkout from v2 to v3 (suggested at
  https://github.com/NCAR/MOM6/pull/231#issuecomment-1347224581 thanks
  to @jedwards4b)

* FMS2: Safe inspection of unlimited dim name

The FMS2 function `get_unlimited_dimension_name` raises a netCDF error
if no unlimited dimension is found.  This is problematic for legacy or
externally created input files which may have not identifed their time
axis as unlimited.

This patch adds a new function, `find_unlimited_dimension_name` which
mirrors the FMS2 function but returns an empty string if none are found.

This is an internal function, not intended for use outside of the
module.

* +Refactor internal_tides interface

  Refactors the internal tide code in MOM_internal_tides and MOM_diabatic_driver
to consolidate it in the MOM_internal_tides module and allow the control
structure for that module to be made opaque.  This includes moving the internal
wave speed diagnostics and the call to wave_speeds or other code setting the
internal wave speeds into propagate_int_tide.  The get_param calls for
INTERNAL_WAVE_CG1_THRESH and UNIFORM_TEST_CG were moved from the diabatic module
to the MOM_internal_tides module. The wave_speed_CS and uniform_test_cg were
removed from diabatic_CS and added to int_tide_CS.  The Nb argument to
propagate_int_tide has been made intent inout, as it is now usually set via the
call to wave_speeds in that routine, but for certain tests it could use the
value passed in from diabatic_driver.

  All answers are bitwise identical, but there are changes to public interfaces
and types, and the order of some entries in the MOM_parameter_doc files and the
available_diags files is changed for some cases.

* +Add fluxes%tau_mag and forces%tau_mag

  Add new allocatable tau_mag arrays to the forcing and mech_forcing types to
hold the magnitude of the wind stresses including gustiness contributions.
There is also a new tau_mag diagnostic.  This same information in tau_mag is
being transformed into ustar, but these changes avoid division by the Boussinesq
reference density (GV%Rho0), and allow for a more accurate calculation of
derived fields when in non-Boussinesq mode, without having to multiply and
divide by GV%Rho0.  There is also a new optional tau_mag argument to
extract_IOB_stresses to support these changes.  These new arrays are not being
used yet in the MOM6 solutions, but they are being allocated and populated in
the routines that set the ustar fields, and they have been tested in changes to
the modules that use ustar that will come in a subsequent commit. This commit
also adds the new RLZ_T2_to_Pa element to the unit_scale_type to undo the
scaling of wind stresses and it makes use of it in some of the new code.  All
answers are bitwise identical, but there are new arrays or elements in three
transparent public types.

* *+Fix problems in mixedlayer_restrat_Bodner

  Fixed several problems with the recently added Bodner mixed layer
restratification parameterization code.

- Corrected the dimensional rescaling in the expressions for psi_mag by adding a
  missing factor of US%L_to_Z.

- 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.

- Revised the expression for dd within the mu function 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 that are
  not included 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.

- Corrected or added variable units in comments in the mixedlayer_restrat
  control structure.

  These could change answers (and avoid NaNs) in some cases with
USE_BODNER23=True, MLE_TAIL_DH > 0 or MLE%TAIL_DH > 0, and there will be changes
to the MOM_parameter_doc files for some cases, but given how recently this code
was added, it is expected that all answers are bitwise identical in the existing
test cases.

* FMS2: New interface to set/nullify_domain

This patch adds wrappers to the set_domain and nullify_domain functions
used in FMS1 for internal FMS IO operations.  These are not used in
FMS2, so the wrapper functions are empty.

This is required to eliminate FMS1 IO dependencies in SIS2.

* *Correct nuopc_cap tau_mag bug

  Correct a recently added bug in the expression for tau_mag in the nuopc_cap
version of convert_IOB_to_forces, where CS%gust(i,j) was used in place of
CS%gust_const, even though the 2-d array was not being set.  This commit changes
answers in some recent versions of the code back to what they had been
previously, and it addresses concerns that had been raised with the first
version of gfdl-candidate-2023-07-03 and its PR to the main version of MOM6.

* Fms2 io read3d slice (#399)

* 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

* Update MOM_variables.F90

---------

Co-authored-by: Alistair Adcroft <[email protected]>
Co-authored-by: Abigail Bodner <[email protected]>
Co-authored-by: Marshall Ward <[email protected]>
Co-authored-by: Robert Hallberg <[email protected]>
Co-authored-by: Marshall Ward <[email protected]>
Co-authored-by: Raphael Dussin <[email protected]>
Co-authored-by: Niki Zadeh <[email protected]>
Co-authored-by: Pavel Perezhogin <[email protected]>
Co-authored-by: Matthew Harrison <[email protected]>
Co-authored-by: Matthew Thompson <[email protected]>

Follow Marshall Ward suggestion and rearrange the code
to be closer to what the compilers will do (or we hope
they would do).
This commit aims to potentially enhance performance.
Answers are bit-wise identical.

…on_tapering

Option to taper neutral diffusion

Changes needed for introducing 3D tracer diffusivities

* Add Leith+E

This commit adds the 2D Leith+E closure, which uses a modified 2D Leith biharmonic viscosity paired with a harmonic backscatter. ('Modified' here is not used in the same sense as 'modified 2D Leith'; it just means that the biharmonic coefficient is modified to account for enstrophy backscatter.) Variables are often named 'leithy' to refer to Leith+E.

The parameterization is controlled by three main entries in user_nl_mom:
1. USE_LEITHY = True
2. LEITH_CK = 1.0
3. LEITH_BI_CONST = 8.0

To use Leith+E you should have LAPLACIAN=True and BIHARMONIC=True. (It doesn't hurt to be explicit and also set LEITH_AH=False, along with any other viscous closures, but this is not required. If USE_LEITHY=True it will not use any of the other schemes. It does use the background value of the biharmonic coefficient as a minimum, but ignores the
background harmonic value.) LEITH_CK is the fraction of energy dissipated by the biharmonic term that gets backscattered by the harmonic term (it's a target; the backscatter rate is not exact.) Recommended values between 0 and 1. LEITH_BI_CONST is Upsilon^6 where Upsilon is the ratio between the grid scale and the dissipation scale for enstrophy.
Values should be greater than or equal to 1; 8 is a good place to start.

The code is sensitive to the background value of Ah; specifically, if Ah is too large, the code is unstable. This is because the backscatter coefficient is proportional to Ah, and if Ah is large then you get large backscatter. If your code is unstable, consider reducing, e.g., `AH_VEL_SCALE`.

* Background Ah

This commit updates the code so that it uses the background Ah as
a minimum. Previously, if `SMAGORINSKY_AH = True`, Leith+E would
use the Smag value of Ah as the minimum, which is incorrect.

* Improve logging

Removed `do_not_log` condition on `USE_LEITHY`

* Fix Leithy Logic

Added one line to fix the fact that the code would only work as
intended if either (i) writing out Ah_h, or (ii) in debug mode.
Also swapped .le. and .lt. for <= and <.

…ugust2023

Add option to apply non-local momentum flux

* Makes set_u_at_v and set_v_at_u public

* First draft for fpmix

* Change name of logical

Replaces LU_pred to L_diag, since now this logical only controls
if diagnostics should be posted.

* Updates to vertFPmix

This commit adds the latest updates to the vertFPmix subroutine
after Bill Large did some cleaning. We have highlight places in
the code where work must be done.

* Add missing use for vertFPmix

* Add omega_w2x to fluxes and forces

omega_w2x is the counter-clockwise angle of the wind stress with
respect to the horizontal abscissa (x-coordinate) at tracer
points [rad]. This variable is needed in the vertPFmix subroutine.

* Add mssing call to get_param for FPMIX

This line of code was lost during the last merge.

* Pass wavebands from coupler to wave_parameters_CS

This commit passes the waveband information recieved
from the coupler to wave_parameters_CS. This information
is set to public so that it can be used elsewhere. To
exercise this code the following must be set:

SURFBAND = COUPLER
WAVE_METHOD = SURFACE_BANDS

No answer changes.

* Describe local variables and make code consistent

* Removed L_diag and moved variables in vertFPmix

* Revert order of variables in vertFPmix

* Move mct_cap/ to STALE_mct_cap/. mct cap is no longer supported and will soon be removed for good.

* remove mct from CI testing

* Remove mct test from github workflows

fix multi-instance PR

Remove fms_io_mod import