- Added a base class in MOM_EOS_base_type.F90
- All EOS modules now extend this base class
  - This reduces replicated code between the EOS modules
- All existing APIs in MOM_EOS now avoid branching for the type of
  EOS and ultimately pass through to a low-level elemental function
  implementation of the actual EOS
- Added a new elemental function exposed by MOM_EOS
  (currently not used in the main model)
- There is a speed up over the previous form of EOS due to the
  reduced branching
  - For some functions, a local implementation of the base class member is
    needed to gain performance. I deliberately did not implement this
    optimization for UNESCO or Jackett06 so that the generic implementation
    of the base class is utilized and we have code coverage.

- Added rules to .testing/Makefile to invoke build.timing, run.timing for the
  "target" code checked out for regression tests
- Appended to existing GH "perfmon" workflow

  Added or corrected the units in comments describing about 200 real variables
in MOM_hor_bnd_diffusion and MOM_neutral_diffusion, and corrected spelling
errors or other issues in about another 20 comments.  Only comments are changed
and all answers are bitwise identical.

  Fix unpaired parentheses in the format statement used in an error message
about OBC segment data not being on the supergrid.  All answers are bitwise
identical, but there is one less compile-time warning.

  Added the new public interfaces continuity_fluxes and continuity_adjust_vel to
make use of the continuity code without actually updating the layer thicknesses,
and made the existing routines zonal_mass_flux and meridional_mass_flux in the
continuity_PPM module public.   Continuity_fluxes is overloaded to provide
either the layer thickness fluxes or the vertically summed barotropic fluxes.
As a part of this change, the LB arguments to meridional_mass_flux and
zonal_mass_flux were made optional and moved toward the end of the list of
arguments.  The intent of the ocean_grid_type argument to 11 routines in the
continuity_PPM module was changed from inout to in.

  Missing factors of por_face_area[UV] were also added to merid_face_thickness
and zonal_face_thickness when the visc_rem arguments are not present or where
OBCs are in use.  This could change answers, but it seems very unlikely that any
impacted cases are in use yet.

  Answers are bitwise identical all known cases, but there are 4 new public
interfaces, and some bugs were fixed in cases that are not likely in use yet.

  In recognition of the fact that there is only a single continuity scheme that
is accessed via the MOM_continuity module, this commit refactors MOM_continuity
to use 3 simple pass-through interfaces (continuity, continuity_init and
continuity_stencil) to the corresponding routines from MOM_continuity_PPM, while
the MOM_continuity control structure is now alias for continuity_PPM_CS and is
opaque in the MOM_continuity module.  As a part of these changes, the runtime
parameter CONTINUITY_SCHEME was obsoleted with a warning value of "PPM".  All
answers are bitwise identical, and all public interfaces look the same from the
outside, but there is one fewer entry in the MOM_parameter_doc.all files.

  Refactored MOM_continuity_PPM to create the separate publicly visible routines
zonal_edge_thickness and meridional_edge_thickness, and also renamed the
internal routines zonal_face_thickness and merid_face_thickness to
zonal_flux_thickness, meridional_flux_thickness and made them publicly visible
as well.

  This commit also renames a number of internal edge thickness variables from
h_L and h_R to h_S and h_N or h_W and h_E for greater clarity, since left and
right are not so well defined on the grid as north, south, east and west.

  All answers are bitwise identical, but there are 4 new publicly visible
interfaces.

  Moved the calls to zonal_edge_thickness and meridional_edge_thickness out of
zonal_mass_flux and meridional_mass_flux to facilitate the reuse at some later
date of the PPM thickness reconstructions.  As a part of this, there are new
edge thickness arguments to zonal_mass_flux and meridional_mass_flux.  The
interfaces to zonal_edge_thickness and meridional_edge_thickness are new
publicly visible and are used in MOM_continuity.

  This commits also changes the name of the loop_bounds_type to
cont_loop_bounds_type and makes it public but opaque adds the publicly visible
function set_continuity_loop_bounds to enable the continuity loop bounds to be
set from outside of the continuity_PPM module.

  Reflecting these changes there are new calls to zonal_edge_thickness and
meridional_edge_thickness in the 3 routines in MOM_continuity and in
continuity_PPM, and new arrays for holding the edge thicknesses in these
routines.

  All answers are bitwise identical, but there are new publicly visible
interfaces and types and changes to other publicly visible interfaces.  However,
no changes are required outside of MOM_continuity and MOM_continuity_PPM.

  Added the new publicly visible routines zonal_BT_mass_flux and
meridional_BT_mass_flux to return the vertically summed transports that the
continuity solver would generate.  Also revised the routine continuity_2d_fluxes
in MOM_continuity to make use of these new routines.  Because these new routines
are not yet being used, all answers are bitwise identical, but there are new
public interfaces.

  Added the new publicly visible routines continuity_zonal_convergence and
continuity_meridional_convergence to increment layer thicknesses using the
continuity loop bounds type to specify extents.  Also revised continuity_PPM
to use these new routines.  These changes will allow for the reuse of some
of the reconstructions in calls that replace calls to continuity with the
unwrapped contents.  All answers are bitwise identical, but there are two new
public interfaces.

  Move continuity_fluxes and continuity_adjust_vel from MOM_continuity.F90 to
MOM_continuity_PPM.F90, but with these interfaces also offered via the
MOM_continuity module so that no changes are required outside of these two
files.  In addtion, 11 of the recently added public interfaces from
MOM_continuity_PPM are also made available as pass-through interfaces from
MOM_continuity.  All answers are bitwise identical.

  Added the new optional arguments du_cor and dv_cor to continuity_PPM and
zonal_mass_flux or meridional_mass_flux to return the barotropic velocity
increments that make the summed barotropic transports match uhbt or vhbt.  Also
cleaned up and simplified the logic of some of the flags used to apply specified
open boundary conditions, adding a new 1-d logical array thereby avoiding
working unnecessarily on some loops or repeatedly checking for specified open
boundary condition points.  Two openMP directives were also simplified.  All
answers are bitwise identical, but there are new optional arguments to three
publicly visible routines.

This patch adds support to makedep for handling most #ifdef-like
condition blocks.  The following preprocessing commands are handled:

* #define / #undef
* #ifdef / #else / #endif

A new flag is added to provide defined macros (-D), and is chosen to
match the `cpp` flag, so that flags can be shared across programs.

Macros are tracked in a new internal variable and are used for #ifdef
testing.  Nested condition blocks are supported by using an internal
stack of the exclusion state.

Certain cases are still not handled.

* #if blocks containing logical expressions are not parsed.
* CPP content inside of #include is ignored.

No doubt many other cases are still unconsidered, such as exotic macro
names.

The autoconf builds use this feature by passing the generated $(DEFS)
argument to makedep.  This is suitable for now, but we may need to
consider two cases in the future:

* Macros defined in $CPPFLAGS are currently ignored, and perhaps they
  should be included here.  The risk is that it may contain non-macro
  flags.

* At some point, DEFS could be moved to a config.h file, and DEFS would
  no longer contain these macros.  (Note that this is very unlikely at
  the moment, since this feature only works with C.)

  Standardized the syntax for the units of salinities that are read via
get_param calls to be uniformly 'units="ppt"' or similar units for derivatives
with salinity.  The only exceptions are places where practical salinity is used
specifically, which occurs for several arguments in the MOM_EOS code.  All
answers are bitwise identical, but there are changes to a number of
MOM_parameter_doc files.

 - Without this, if part of your OBC is filled with land mask and
   if that land mask contains a masked out tile, you will generate
   a NaN from the phase speed calculation where h is negative in the
   halo neighbor of that masked tile.

  Renamed the runtime parameter FIX_USTAR_GUSTLESS_BUG to USTAR_GUSTLESS_BUG
(with a switch between the meanings of true and false for the two parameters)
for consistency with the syntax of other bug-fix flags in MOM6 and to partially
address dev/gfdl MOM6 issue #237.  Input parameter files need not be changed
right away because MOM6 will still work if FIX_USTAR_GUSTLESS_BUG is specified
instead of USTAR_GUSTLESS_BUG, but USTAR_GUSTLESS_BUG will be logged, so there
are changes to the MOM_parameter_doc files.   By default or with existing input
parameter files, all answers are bitwise identical, and there is error handling
if inconsistent settings of FIX_USTAR_GUSTLESS_BUG and USTAR_GUSTLESS_BUG are
both specified.

…d/flipped input and grids, and under dimensional rescaling to the +/-140 power. Bitwise identical results were achieved under rotation/rescaling using standard 2D test cases such as MISMIP+. Major changes include the following:

- Specified order of operations throughout, and modified the implementation of FEM integration,  h-point-to-node operations, etc for consistency under rotation (see MOM_ice_shelf_dynamics.F90 subroutines calc_shelf_driving_stress, CG_action, CG_action_subgrid_basal, matrix_diagonal, shelf_advance_front, and interpolate_h_to_b).
- Added an ice-shelf version of 'first_direction' and 'alterate_first_direction', allowing users to control the order in which the x- and y-direction ice-shelf advection calls are made. Required for consistency under rotation.
- Dimensional rescaling fixes throughout MOM_ice_shelf_initialize.F90, and within MOM_ice_shelf_dynamics.F90 (see subroutines initialize_ice_shelf_dyn, ice_shelf_solve_outer, and calc_shelf_taub)
- Rotation/rescaling testing revealed two additional bugs that were subsequently fixed:
(1) An index error in the conjugate gradient algorithm (subroutine CG_action)
(2) Discretization error for the east/north fluxes in subroutines ice_shelf_advect_thickness_x and ice_shelf_advect_thickness_y.

The commit also includes several simple changes for code readability and computational efficiency:

- Saved Phi, PhiC, and Phisub in the ice shelf dynamics control structure at the start of each run, rather than reallocating and recalculating these static fields each time they are used. Reshaped the Phisub array for computational efficiency over loops.
- Modified the sub-cell grounding scheme so that it is only called for partially-grounded cells that contain the grounding line (i.e. where float_cond(i,j)==1). Added float_cond to the ice dynamics structure for output.
- Simplified the option to calculate ice viscosity at 4 quadrature points per cell rather than at cell centers, by adding parameter NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS and reshaping CS%ice_visc accordingly.
- Style changes and removal of unused code (e.g. removed subroutine apply_boundary_values and field thickness_bdry_val)

…ficiency

…ion law

  Renamed the arguments u and v to step_MOM_dyn_split_RK2 as u_inst and v_inst
to more clearly differentiate between the instantaneous velocities (u_inst and
v_inst) and the velocities with a time-averaged phase in the barotropic mode
(u_av and v_av).  A comment is also added at one point where the wrong
velocities are being used to calculate and apply the Orlanski-style radiation
open boundary conditions, with the intention of adding the option to correct
this in a subsequent commit.  This commit only changes the name of a pair of
internal variables in one routine, and all answers are bitwise identical.

  Renamed the runtime parameter FIX_UNSPLIT_DT_VISC_BUG to UNSPLIT_DT_VISC_BUG
(with a switch between the meanings of true and false for the two parameters)
for consistency with the syntax of other bug-fix flags in MOM6 and to partially
address dev/gfdl MOM6 issue #237.  Input parameter files need not be changed
right away because MOM6 will still work if FIX_UNSPLIT_DT_VISC_BUG is specified
instead of UNSPLIT_DT_VISC_BUG, but UNSPLIT_DT_VISC_BUG will be logged, so there
are changes to the MOM_parameter_doc files.   By default or with existing input
parameter files, all answers are bitwise identical, and there is error handling
if inconsistent settings of FIX_UNSPLIT_DT_VISC_BUG and UNSPLIT_DT_VISC_BUG are
both specified.

  Add the new module MOM_dynamics_split_RK2b and calls to the routines in this
module to MOM.F90.  These calls are exercised when the run time parameter
SPLIT_RK2B is true.  For now, all answers are bitwise identical when this new
code is being exercised, but there is a new module with multiple public
interfaces and a new entry in many MOM_parameter_doc files.

  Revised step_MOM_dyn_split_RK2b to use the time-filtered velocities as
arguments and reconstruct the instantaneous velocities with the proper phase
from the barotropic solver from the saved increments du_av_inst and dv_av_inst.
As a part of this change, the continuity solver is used to find the thickness
fluxes used in the predictor step Coriolis terms, and the horizontal viscous
accelerations are calculated for both the predictor and corrector steps, thereby
avoiding the need to vertically remap any 3-d fields apart from a single pair of
velocity components.

  The run-time option STORE_CORIOLIS_ACCEL is no longer used when SPLIT_RK2B =
True.  FPMIX was also disabled in this mode due to unresolved dimensional
inconsistency errors in that code.

  Additionally, the time-filtered velocities are now used instead of the
instantaneous velocities in the second call to radiation_open_bdry_conds(),
which should improve the performance of several of the Orlanski-type open
boundary conditions.

  The 2-d fields du_av_inst and dv_av_inst were added to the restart file, while
the 3-d fields u2, v2, diffu and diffv and either CAu and CAv or uh, vh and h2
are all removed from the restart files.  Remap_dyn_split_RK2b_aux_vars() now has
nothing to do, and it should probably be eliminated altogether in a subsequent
commit.

  All answers are changed when SPLIT_RK2B = True, and there are changes to the
contents of the restart files.

  This commit includes further revisions to MOM_dynamics_split_RK2b that avoid
some unnecessary calculations and group some of the halo updates into fewer
group passes.  All answers are bitwise identical.

  Corrected the description of the runtime parameter SPLIT_RK2B that will appear
in the MOM_parameter_doc files and in the doxygen descriptions of the MOM module
to better reflect what was ultimately being done with this new scheme.  All
answers are bitwise identical, but there are changes to the (newly added)
contents of some MOM_parameter_doc files.

  Added the new functions cuberoot and intrinsic_functions_unit_tests to the
MOM_intrinsic_functions module, and call intrinsic_functions_unit_tests from
unit_tests to confirm that this new function works as intended.  Separately,
cuberoot was tested by replacing expressions like A**(1./3.) with cuberoot(A) in
MOM_energetic_PBL and verifying that the answers only change at roundoff, but
that it can give bitwise identical results when the argument is scaled by an
integer power of 8 and then unscaled by the corresponding integer power of 2,
but that change will occur in a subsequent commit as it can change answers
depending on an ANSWER_DATE flag.  With this commit, cuberoot is not yet being
used so all answers are bitwise identical, although there are new publicly
visible routines.

JOB_DIR in the Gaea-specific .gitlab-ci.yml configuration file is
updated to use its F5 filesystem.

  Modified the cuberoot function to do 3 iterations with Halley's method
starting with a first guess that balances the errors at the two ends of the
range of the iterations, before a final iteration with Newton's method that
polishes the root and gives a solution that is accurate to machine precision.
Following on performance testing of the previous version, all convergence
testing has been removed and the same number of iterations are applied
regardless of the input value.  This changes answers at roundoff for code that
uses the cuberoot function, so ideally this PR would be dealt with before the
cuberoot becomes widely used.

Applying the first iteration explicitly appears to speed up the cuberoot
function by a bit over 20%:

 Before:
 Halley Final:  0.14174999999999999

 After:
 Halley Final:  0.11080000000000001

There is an assumption that compilers will precompute the constants like
`0.7 * (0.7)**3`, and that all will do so in the same manner.

* Update forcing.rst

* add more detail to forcing.rst removing ambiguity about yaml format.

This patch adds `do_not_log` to `openParameterBlock`, to prevent logging
of `BLOCK%` `%BLOCK` entry and exit calls.

The argument was not added to `closeParameterBlock`, since this state
can be tracked inside the `block` with a new `log_access` field in
`parameter_block`.

This flag does not extend to parameters within the block, since (as far
as I know) there is no way for a `get_param` to know if it is within a
block or not.  Even if it could know this, there would need to be some
careful handling of nested blocks.

The potential block/parameter inconsistency should be supported at some
point, but for now it is the user's responsibility to consistently apply
`do_not_log` to blocks and its contents.

Fix the error caused by readinging `MLE%BODNER23`, and instead
explicitly opening and closing the parameter blocks.

  Added a missing unit conversion factor to a hard-coded 10 m distance in
oil_tracer_column_physics.  This will not change answers in Boussinesq cases
without any dimensional rescaling, but it will correct answers in a hypothetical
non-Boussinesq case.  Also made some white space in expressions near this fix
more closely match the MOM6 style guide. No answers are affected in any known
existing regression test cases or other runs.

  Use GV%H_to_MKS instead of GV%H_to_m when undoing the dimensional rescaling of
thicknesses when taking weighted averages in horizontally_average_diag_field,
global_layer_mean and global_volume_mean.  In Boussinesq mode, these are
identical, but in non-Boussinesq mode using GV%H_to_m introduced a
multiplication and then division by the Boussinesq reference density, whereas
GV%H_to_MKS avoids this by rescaling to a volume or mass-based coordinate
depending on the mode.  Several comments were also updated to reflect these
conditional changes in the units of some internal variables.  All expressions
are mathematically equivalent, and this does not impact any solutions, but there
can be changes in the last bits in some non-Boussinesq averaged diagnostics.

This patch replaces the intrinsic-based exponent rescaling with explicit
bit manipulation of the floating point number.

This appears to produce a ~2.5x speedup of the solver, reducing its time
from embarassingly slow to disappointingly slow.  It is slightly faster
than the GNU cbrt function, but still about 3x slower than the Intel
SVML cbrt function.

Timings (s) (16M array, -O3 -mavx -mfma)

| Solver              |  -O2  |  -O3  |
|---------------------|-------|-------|
| GNU x**1/3          | 0.225 | 0.198 |
| GNU cuberoot before | 0.418 | 0.412 |
| GNU cuberoot after  | 0.208 | 0.187 |
| Intel x**1/3        | 0.068 | 0.067 |
| Intel before        | 0.514 | 0.507 |
| Intel after         | 0.213 | 0.189 |

At least one issue here is that Intel SVML is using fast vectorized
logic operators whereas the Fortran intrinsics are replaced with slower
legacy scalar versions.  Not sure there is much we could even do about
that without complaining to vendors.

Also, I'm sure there's magic in their solvers which we are not
capturing.  Regardless, I think this is a major improvement.

I do not believe it will change answers, but probably a good idea to
verify this and get it in before committing any solutions using
cuberoot().

Some modifications were made to the cuberoot rescale and descale
functions:

* The machine parameters were moved from function to module parameters.
  This could dangerously expose them to other functions, but it prevents
  multiple definitions of the same numbers.

* The exponent is now cube-rooted in rescale rather than descale.

* The exponent expressions are broken into more explicit steps, rather
  than combining multiple steps and assumptions into a single
  expression.

* The bias is no longer assumed to be a multiple of three.  This is true
  for double precision but not single precision.

A new test of quasi-random number was also added to the cuberoot test
suite.  These numbers were able to detect the differences in GNU and
Intel compiler output.  A potential error in the return value of the
test was also fixed.

The volatile test of 1 - 0.5*ULP has been added.  The cube root of this
value rounds to 1, and needs to be handled carefully.

The unit test function `cuberoot(v**3)` was reversed to `cuberoot(v)**`,
to include testing of this value.  (Cubing would wipe out the anomaly.)

In separate testing, we observed that Intel would use the `pow()`
function to evaluate the cubes of some numbers, causing different
answers with GNU.

In this patch, I replace the cubic x**3 operations with explicit x*x*x
multiplication, which appears to avoid this substitution.

Well, for the moment, at least.

  Use the new cuberoot() function in place of **(1./3.) to calculate the
turbulent velocity vstar in ePBL_column when EPBL_ANSWER_DATE is 20240101 or
higher.  This is mathematically equivalent to the previous version, but it does
change and answers at roundoff and it allows several dimensional scaling factors
that had previously been required to be eliminated.  All answers are
mathematically equivalant, but answers do change if EPBL_ANSWER_DATE is 20240101
or higher and the description of EPBL_ANSWER_DATE changes in some
MOM_parameter_doc files.

  Fixed two bugs on a single line of convert_temp_salt_for_TEOS10.  The first
bug was a reversal in the order of the temperature and salinity arguments to
poTemp_to_consTemp, resulting in temperatures that closely approximate the
salinities.  The second bug that was fixed on this line was temperatures being
rescaled with a factor that is appropriate for salinities. This bug-fix will
change answers dramatically for any cases that use the ROQUET_RHO, ROQUET_SPV
and TEOS10 equations of state and initialize the model with
INIT_LAYERS_FROM_Z_FILE = True.

  This commit add the new optional arguments h_neglect and h_neglect_edge to
ALE_remap_scalar to allow for the thicknesses used in this routine to be
provided in any self-consistent units, including [Z ~> m], instead of just
[H ~> m or kg m-2].  To help make use of this new capability, this commit also
adds the new functions set_h_neglect and set_dz_neglect to the MOM_regridding
module.  build_grid_rho and build_grid_HyCOM1 have been refactored to use
set_h_neglect in place of the corresponding duplicated code blocks.

  This commit also adds the new optional argument h_in_Z_units to
MOM_initialize_tracer_from_Z, which in turn uses this new capability for
ALE_remap_scalar to use vertical layer extents (in Z units) rather than
thicknesses (in H units).

  Although there are new optional arguments to public interfaces, they are not
yet being exercised with this commit so no answers are changed.  Moreover, even
if they were being exercised, all Boussinesq solutions would give identical
answers.

  Revise MOM_temp_salt_initialize_from_Z in cases when Z_INIT_REMAP_GENERAL is
False to call ALE_remap_scalar with vertical layer extents (in Z units) rather
than layer thicknesses (in H units).  When in fully non-Boussinesq mode, this
same routine uses dz_to_thickness (using the full equation of state) rather than
dz_to_thickness_simple to initialize the layer thicknesses.  Boussinesq answers
are bitwise identical, but answers can change in some fully non-Boussinesq
cases.

  Revised initialize_MOM_generic_tracer to use thickness_to_dz to get the layer
vertical extents and then provide these to MOM_initialize_tracer_from_Z to read
in initial generic tracer concentrations from Z-space files.  The previous
approach inappropriately used an simple multiplicative rescaling (via a call to
dz_to_thickness_simple in MOM_initialize_tracer_from_Z) by a factor that
includes the Boussinesq reference density when in non-Boussinesq mode.  A new
thermo_vars_type arguments was added to initialize_MOM_generic_tracer to allow
for this change.

  Also revised MOM_generic_tracer_column_physics to use thickness_to_dz instead
of a simple multiplicative rescaling to get the layer vertical extents (in m)
that are used in calls to generic_tracer_source.  The multiplicative factor that
was used previously (GV%H_to_m) includes the Boussinesq reference density and
hence is inappropriate in non-Boussinesq mode; using thickness_to_dz avoids
this.

  Also added comments documenting the meaning and units of about 30 real
variables in the MOM_generic_tracer routines.

  There is a new mandatory argument to initialize_MOM_generic_tracer.  All
Boussinseq mode answers are bitwise identical, but in non-Boussinesq mode mode
generic tracer answers are changed by avoiding the use of the Boussinesq
reference density in several places.

Moved call for filename_appendix outside of loop.

Fixes issue where restart filename was
./MOM.res.ens_01.nc
./MOM.res.ens_01.ens_01_1.nc
./MOM.res.ens_01.ens_01.ens_01_2.nc
./MOM.res.ens_01.ens_01.ens_01.ens_01_3.nc

instead of
./MOM.res.ens_01.nc
./MOM.res.ens_01_1.nc
./MOM.res.ens_01_2.nc
./MOM.res.ens_01_3.nc

  Use the new cuberoot function in the Bodner estimate of u'w' when the new
runtime parameter ML_RESTRAT_ANSWER_DATE is 20240201 or higher, which avoids the
need to undo and redo the dimensional scaling in this calculation.  In addition,
the refactoring associated with this change explicitly revealed that as it was
implemented it introduced an undesirable dependency on the value of the
Boussinesq reference density, RHO_0, when in non-Boussinesq mode.  To avoid
this, a new version of the Bodner u'w' calculation was introduced in fully
non-Boussinesq mode, which does change answers with this combination; because
there is not yet a known case that used this combination, we have chosen not to
add a runtime parameter to preserve the old answers when the Bodner
parameterization is used in fully-Boussinesq mode.  This change will modify the
contents of some MOM_parameter_doc files with USE_BODNER=True, and it changes
answers in cases that are also fully non-Boussinesq.  The new runtime parameter
ML_RESTRAT_ANSWER_DATE might need to be set below 20240201 to retain some
existing Boussinesq answers.

  Eliminated some unused calculations in wave_speed and thickness_diffuse_full
that would use the Boussinesq reference density when in non-Boussinesq mode,
either by adding logical flags around the calculations so that they only occur
in Boussinesq mode.  All answers are bitwise identical.

  Added or amended comments describing about 90 variables in the MOM_regridding,
coord_adapt, coord_rho, and coord_sigma modules to document the purpose and
units of real variables using the MOM6 standard syntax.  A missing doxygen
comment was also added to describe inflate_vanished_layers_old().  A few unused
internal variables were also eliminated, including the max_depth_index_scale
element of the regridding_CS type.  All answers are bitwise identical and only
comments are changed.

  Modified comments to document the units of 17 variables in 5 files in the core
directory and another 5 variables in MOM_EOS_linear.F90.  Only comments are
changed and all answers are bitwise identical.

  Revised the ALE_sponge code to convert the model's thicknesses (in thickness
units of [H ~> m or kg m-2]) into layer vertical extents (in height units of
[Z ~> m]) to match how the input datasets vertical grids are specified and to
avoid using the Boussinesq reference density to do the opposite conversion when
in non-Boussinesq mode.  The internal conversion in apply_ALE_sponge between
layer thicknesses and layer vertical extents matches those in thickness_to_dz,
and uses the layer-averaged specific volumes when in non-Boussinesq mode.

  As a part of these changes, 12 internal variables in MOM_ALE_sponge were
renamed (e.g., from data_h to data_dz) to more clearly reflect that they are the
layer vertical extents instead of thicknesses, and there are new 3-d variables
dz_model, that is created from h when the sponges are applied to the velocities,
and data_dz to handle the conversion of the input thicknesses or vertical
extents in initialize_ALE_sponge_fixed.

  This commit includes adding a new thermo_var_type argument to
apply_ALE_sponge, a new unit_scale_type argument to rotate_ALE_sponge and
initialize_ALE_sponge_varying and the new optional argument data_h_is_Z to
initialize_ALE_sponge_fixed to indicate that the in thicknesses arguments are
vertical distances (in [Z ~> m]) and do not need to be rescaled inside of the
ALE_sponge code, and similarly for the new optional argument data_h_in_Z to
get_ALE_sponge_thicknesses.

  This commit also adds halo size and Omit_Corners arguments to the pass_var
calls in the ALE_sponge code.

  There are new optional arguments to 2 publicly visible routines, and new
mandatory arguments to 3 others, but by default the meaning and units of all
previous arguments are the same, and answers will be bitwise identical in any
Boussinesq cases.  However, answers will change in any non-Boussinesq cases that
use the ALE_sponges.

  Use thickness_to_dz to convert the layer thicknesses into the geometric
depths used to set the Bryan-Lewis background diffusivities instead of
multiplication by GV%H_to_m, thereby avoiding division by the Boussinesq
reference density in non-Boussinesq mode. No answers are changed in any
Boussinsesq configurations, but answers will change in non-Boussinesq
configurations that have BRYAN_LEWIS_DIFFUSIVITY = True.

  Added the new optional scale argument to tracer_Z_init to allow the units of
the tracers that it initializes from a z-space file to be rescaled, similarly to
what is already done in tracer_z_init_array.  This also required the addition of
the new missing_scale argument to the private routine read_Z_edges to specify
how the missing values are rescaled. Comments were also added or modified
describing 28 real variables or their units in this module.  By default all
answers are bitwise identical, but when this new option is exercised, answers
could change at roundoff for some non-Boussinesq cases.

* Update MOM_mixed_layer_restrat.F90

Adding and calculating diagnostic front length scale (lf_bodner) to
mixed layer restratification code

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

* Bodner diag: Fix ustar and dims

Fix two issues which emerged from the Bodner length scale diagnostic
after the cuberoot ANSWER_DATE was introduced.

The split of the loop due to the cuberoot caused ustar to be assigned in
a separate loop, leaving it uninitialized for the length scale
diagnostics loop.  It also had its dimensions removed.

However, the whole point of the new cuberoot() function was to preserve
its dimensions, so there was no need for the intermediate dimensional
manipulation, and all dimensionality should be shifted to the
registration.

This patch removes the intermediate dimensions and also explicitly
defines it as Z2/H rather than L.

* Bodner diagnostic: Restore scaling to L

The scaling of the Bodner lengthscale diagnostic is correctly restored
to L in this patch.  The Z2 H-1 constant is factored into the
calculation of the lengthscale; I believe this represents a conversion
from a [Z2 T-2]-like momentum flux to a [LH T-2]-like momentum flux.

A generic `uflux_rescale` term was introduced to be used in both
ld_bodner and wpup, to avoid a double division by SpV.

The if-block for computation of ld_bodner was also moved outside of the
loop, in order to facilitate vectorization.

---------

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

Inside the `MLE%` ... `%MLE` parameter block, the default answer date
was being read while inside the block.  Since the namespaces inside
blocks are not shared or inherited, this was creating a new parameter,
`MLE%DEFAULT_ANSWER_DATE`, whose value differed from the global
`DEFAULT_ANSWER_DATE`.  This was causing unintuitive behavior of
`MLE%ML_RESTRAT_ANSWER_DATE`, which ignored this global answer date.

This patch reads `DEFAULT_ANSWER_DATE` outside of the block, so that
`ML_RESTRAT_ANSWER_DATE` now defaults to the global value.

  Added explicit names to the h_neglect and h_neglect_edge arguments in a call
to build_rho_column in diag_remap_update.  These were problematic because these
arguments were previously in the wrong positions and hence were being used for
the wrong variables inside of build_rho_column.  In Boussinesq mode, the results
will be the same as long as there is no dimensional rescaling being applied to
thicknesses, and it this were being tested it would have resulted in a
dimensional consistency test failure.  In non-Boussinesq mode, the values of
H_subroundoff as used for this coordinate will be off by a factor of
RHO_KV_CONVERT (usually 1035), but these thicknesses are so small that they do
not matter unless ANGSTROM is explicitly set to zero (when no dimensional
rescaling is being used).  To summarize, this fixes a bug that has been around
since 2017, and while this fix could slightly change some diagnostics that are
remapped to density, it only does so in circumstances that are unlikely to have
been realized.  No answer changes were detected in the MOM6-examples regression
suite, and almost certainly none ever would be in any existing use of MOM6.

  Revised the update_OBC_segment_data code to keep the z-space input data in
height units ([Z ~> m]) rather than rescaling them to thickness units.  This
change means that the non-Boussinesq open boundary condition calculations avoid
using the Boussinesq reference density.  Associated with this change, 5 internal
variables in update_OBC_segment_data were renamed to reflect that they are layer
vertical extents instead of thicknesses.  Also added or amended comments
describing the purpose and units of 25 real variables in this module and
corrected a handful of lines in this module that were not adhering to the
guidance in the MOM6 style guide.  Answers are bitwise identical in any
Boussinesq cases.  However, answers will change in any non-Boussinesq cases that
use open boundary conditions that use time_interp_extrnal to read in Z-space
data.

  Removed the duplicative multiplication by the Boussinesq density and related
unit conversion factors in the calculation of the diagnosed mole concentration
of the "cfc11" and "cfc12" diagnostics.  This unit rescaling is duplicative of
the conversion factors in the register_diag_field calls for these diagnostics
(at about line 263).  This will change these two diagnostics by a factor of
RHO_0.  This bug arose from separate and independent changes coming from GFDL
and NCAR to add the same missing rescaling factor, and it might have gone
undetected because of the vast range of valid values for CFC concentrations or
because all of our regression testing uses calendar dates that precede the
invention and first release of CFCs in the 1930s.  Although this commit will
correct the code that is on the dev/gfdl branch, care will have to be taken when
merging in modified versions of this file from other branches that might have
also corrected this bug to avoid reintroducing it (perhaps in reverse).

  Also replaced the expression for the CFC flux rescaling factor for use with
KPP_NonLocalTransport with GV%RZ_to_H, which is equivalent to the previous
expression in Boussinesq mode but avoids a multiplication and division by the
Boussinesq reference density in non-Boussinesq mode.  All Boussinesq answers are
identical, but some non-Boussinesq CFCs could be altered in the last bits.

  This PR will change diagnostics of CFC concentrations by about 3 orders of
magnitude. The internal representation of the CFCs is bitwise identical in all
Boussinesq cases, but they may change at roundoff in non-Boussinesq cases that
use KPP.

Minor version updates to multiple packages used in the generation of the
documentation introduce dependencies needing Sphinx >= 5.0, which breaks
the sphinx extensions we use in documenting the MOM6 APIs. I have added
versions for all the packages needed to keep things working with Sphinx4
for now, but we really do need to find a way to work with the newer
versions.

More pinningaof of requirements for readthedocs

More pinning of requirements for readthedocs

This patch introduces two new macros, BUILD and WORK, to permit
relocation of the build/ and work/ directories.  It also makes the
following smaller changes:

* deps/ is now defined by the DEPS macro.  If unset, deps/ is placed in
  the BUILD directory.

* results/ is moved into WORK.

* Compiler flags which track directories now use $(abspath ...) to
  allow for arbitrary paths.

* GitHub CI paths were adjusted to support these new settings.

* DO_* flags are now used as on/off with ifdef testing, rather than
  checking for `true` values.

* mkmf macros have been removed from the coupled test config.

* The default FMS infra has been changed to FMS2 in all components,
  including the configure.ac outside of .testing.

This work will enable testing of multiple FMS libraries in our CI.

  Refactored set_viscous_BBL to separate out the routines setting the open
interface lengths used for the channel drag, shortening a 1070 line long routine
to 915 lines and reducing the scope of a number of temporary variables.  A
number of logical branch points have been moved outside of the innermost do
loops.  This refactoring will also make it easier to provide alternatives to
some of the solvers that do not use the trigonometric functions to solve for the
roots of a cubic expression and avoiding the issues noted at
/issues/483.  All answers are bitwise identical and public
interfaces are unchanged.

  Added the new routine find_L_open_concave_iterative to use iterative Newton's
method approaches with appropriate limits to solve the cubic equation for the
fractional open face lengths at interfaces that are used by the CHANNEL_DRAG
code.  These solutions are analogous to those given by the previous expressions
that are now in find_L_open_concave_trigonometric, and the two differ at close
to roundoff, but the new method is completely independent of the transcendental
function library, thereby addressing dev/gfdl MOM6 issue #483.  This new routine
is called when the new runtime parameter TRIG_CHANNEL_DRAG_WIDTHS is set to
false, but by default the previous answers are recovered.  By default all
answers are bitwise identical, but there is a new runtime parameter in some
MOM_parameter_doc files.

  Added the new debugging or testing subroutine test_L_open_concave along with
extra calls when DEBUG = True that can be used to demonstrate that the iterative
solver in find_L_open_concave_iterative is substantially more accurate but
mathematically equivalent to the solver in find_L_open_concave_trigonometric.
This extra code is only called in debugging mode, and it probably should be
deleted in a separate commit after find_L_open_concave_iterative has been
accepted onto the dev/gfdl branch of MOM6.  All answers are bitwise identical
and no output or input is changed.

  Carried out minor refactoring in set_viscous_BBL as suggested by the reviews
of this PR, including the elimination of some unnecessary error handling and the
replacement of C2pi_3 as an argument to find_L_open_concave_trigonometric with
an internal parameter.  All answers are bitwise identical and there are no
changes to publicly visible interfaces.

It was blowing up with "forrtl: error (65): floating invalid" when
accessing dz in the halo at the boundary, but just sometimes.  My
default layout is trouble while my testing layout of 48 cores is not.

  Revised the interfaces to the myStats routine in the horizontal_regridding
module to avoid segmentation faults due to inconsistent horizontal indices and
array extents in global indexing mode.  Rather than passing in absolute array
extents to work on, an ocean grid type is now passed as an argument to myStats,
with the new optional full_halo argument used to capture the case where the
tracer statistics are being taken over the full data domain.  The most
frequently encountered problems occurred when the hard-coded debug variable in
the horiz_interp_and_extrap_tracer routines are changed from false to true.
When global indexing is not used, this revised work exactly as before, but when
it is used with global indexing, it avoids segmentation faults that were
preventing the model from running in some cases with all debugging enabled.

  Corrected bugs in the horizontal indexing in apply_topography_edits_from_file
that led to differing answers depending on the value of GLOBAL_INDEXING.  This
change gives identical results when GLOBAL_INDEXING is used, as can be seen by noting that G%idg_offset = G%isd_global + G%isd and that without global
indexing G%isd = 1, but with global indexing the new expressions give the same
answers as without them.  Because global indexing is not typically used,
answers are not changed for any cases in the MOM6-examples test suite.

  Moved the post-initialization halo updates for thicknesses and temperatures
and salinities in initialize_MOM to occur immediately after the last point where
they are modified and before the remapped diagnostic grids are set up.  This
does not change any existing answers but it will enable the future use of the
thermo_var_ptr type in calls to thickness_to_dz when setting up remapped
diagnostics, and perhaps elsewhere in the initialization of other auxiliary
variables.  All answers are bitwise identical.

  Refactored MOM_diag_remap to work with global indices and to move logical
branches outside of do loops.  This was done by adding internal routines that
set the loop indices consistently with the NE c-grid convention used throughout
the MOM6 code and converting the optionally associated mask pointer into an
optional argument.  It also simplifies the logic of many of the expressions
within the remapping code.  There is also a new element, Z_based_coord, in the
diag_remap_ctrl type, to indicate whether the remapping is working in thickness
or height units, but for now it is always set to false.

  The function set_h_neglect or set_dz_neglect is used to set the negligible
thicknesses used for remapping in diag_remap_update and diag_remap_do_remap,
depending on whether the remapping is being done in thickness or vertical height
coordinates.  Diag_remap_init has a new vertical_grid_type argument, and
diag_remap_do_remap has a new unit_scale_type argument.

  For REMAPPING_ANSWER_DATES later than 20240201, diag_remap_updated does an
explicit sum to determine the total water column thickness, rather than using
sum function, which is indeterminate of the order of the sums.  For some
compilers, this could change the vertical grids used for remapping diagnostics
at roundoff, but no such change was detected for any of the compilers used with
the MOM6 regression test suite.

  All answers and diagnostics in cases that worked before are bitwise identical,
but there are new arguments to two publicly visible interfaces.

  Do remapping in Z-space for some remapped diagnostics, depending on which
coordinate is used.  The subroutine thickness_to_dz is used with the
thermo_vars_type to do the rescaling properly in non-Boussinesq mode.

  A new thermo_var_ptrs argument was added to diag_mediator_init, replacing
three other arguments that had the same information, and its call from MOM.F90
was modified accordingly.

  The various calls to the diag_remap routines from post_data_3d and
diag_update_remap_grids were modified depending on whether a z-unit or h-unit
vertical grid is being remapped to.

  All answers and diagnostics are identical in Boussinesq mode, but some
remapped diagnostics are changed (by not using expressions that depend on the
Boussinesq reference density) in the non-Boussinesq mode.  There are altered
or augmented public arguments to two publicly visible routines.

  Corrected indexing problems in downsample_mask that would cause masked
reduced-resolution diagnostics to work improperly when the model is in symmetric
memory mode or when global indexing is used.  This involves passing the starting
index in memory of the native-grid field being downsampled in all calls to
downsample_mask.  In global-indexing mode, this change avoids a series of
segmentation faults that stopped the model runs when compiled for debugging.
All solutions are bitwise identical in all cases but some down-scaled
diagnostics will change in some memory modes, hopefully becoming consistent
across all memory modes (although this has yet to be tested in all modes).

  Replaced the IIm1 and JJm1 variables in Update_Stokes_Drift with 'I-1' and
'J-1' in Update_Stokes_Drift.  The previous expressions could give incorrect
solutions at the southern and western edges of the global domain with global
indexing and serve no purpose when global indexing is not used.  In addition,
the i-, j- and k- index variables in Update_Surface_Waves, Update_Stokes_Drift,
get_Langmuir_Number and Get_SL_Average_Prof were changed from ii to i, jj to j,
and kk to k to follow the patterns used elsewhere in the MOM_wave_interface
module and throughout the rest of the MOM6 code.  All answers are bitwise
identical in cases that do not use global indexing.

  Added or amended comments to document the mostly arbitrary units of about 360
variables in 13 low-level remapping modules.  Only comments are changed and
all answers are bitwise identical.

This patch contains several bugfixes associated with the rotational grid
testing.

The following checksums are declared as scalar pairs, rather than
vectors:

* eta_[uv] (in porous barrier)
* por_face_area[UV]
* por_layer_width[UV]
* Ray_[uv] (Rayleigh drag velocity)

Fluxes and surface fields are now permitted to contain tracer fluxes
(tr_fluxes) when rotation is enabled.  The fields are retained in their
unrotated form, since these are accessed and handled outside of MOM6.

The rotated p_surf_SSH pointer in forces now correctly points to either
p_surf or p_surf_full.

read_netcdf_nc() now correctly uses the unrotated horizontal index
struct HI, used to access the contents of the file.  Previously, it was
using the model HI, which may be rotated.

Reading chlorophyll with time_interp_external now uses rotation to
correctly fetch its output.

NOTE: This could be cleaned up so that the rotation details are hidden
from users, but there are some unresolved issues around how to approach
this.

  Added the option to do neutral tracer diffusion with expressions that satisfy
rotational symmetry.  This option is enabled by setting the new runtime
parameter NDIFF_ANSWER_DATE to be greater than 20240330.  By default, this
parameter is set to use the previous expressions, although the default may be
changed later to follow DEFAULT_ANSWER_DATE.  By default all answers are bitwise
identical, but there are changes to some MOM_parameter_doc files due to the
introduction of a new runtime parameter.

Runtime parameter USE_CONT_THICKNESS is added in hor_visc to let it use
velocity-point thickness consistent with the continuity solver.
Thicknesses are borrowed from BT_cont so only split mode is supported.

  When ROTATE_INDEX is true, the model keeps both its rotated internal grid and
an unrotated grid for setting initialization and forcing.  When the model is run
in symmetric memory mode, there is an auxiliary non-symmetric domain that is
needed for reading in fields at velocity points.  The auxiliary domain in the
unrotated grid (G_in) was not being set previously, causing the model to fail to
run some cases when ROTATE_INDEX was true.  That auxiliary domain in the
unrotated grid is now being set up properly.  All answers and output are bitwise
identical for any cases that worked previously.

  Added code to tht two versions of clone_MD_to_MD in the FMS1 and FMS2 versions
of MOM_domain_infra.F90 to properly change the flags when a tripolar grid is
rotated, so that it does not lead to misleadingly incorrect answers.  With the
current versions of FMS, doing grid rotation testing with a tripolar grid will
lead to error messages about the incomplete implementation of FMS, but these
failures are preferable to the model silently working incorrectly.  All answers
are bitwise identical in cases that worked correctly before.

   Add get_netcdf_filename and use it to add useful details (the field and
filenames in question) to a fatal error message in get_field_nc.  All answers
are bitwise identical, but there is a new public interface and some output is
changed in cases where get_field_nc is failing.

  Added the option to do epipycnal tracer diffusion between a bulk mixed layer
and the interior ocean with expressions that satisfy rotational symmetry.  This
option is enabled by setting the new runtime parameter HOR_DIFF_ANSWER_DATE to
be greater than 20240330.  By default, this parameter is set to use the previous
expressions, although the default may be changed later to follow
DEFAULT_ANSWER_DATE.  Also corrected two bugs with the tracer limits used to
repartition the fluxes between layers within tracer_epipycnal_ML_diff; this
correction is enables by setting the new HOR_DIFF_LIMIT_BUG to .false., but to
retain previous answers by default it is set to true. By default all answers are
bitwise identical, but there are changes to some MOM_parameter_doc files due to
the introduction of two new runtime parameters.

  Pass vertical extents (dz in [Z ~> m]) instead of thicknesses (h in
[H ~> m or kg m-2] and use data_h_is_Z flag in calls to initialize_ALE_sponge
calls from 5 user modules and avoid extra calls to dz_to_thickness in these
routines.  All Boussinesq solutions are bitwise identical, but in any
non-Boussinesq configurations using ALE sponges, the previous conversion from dz
to thickness and back again can change dz in the last bits, so some
non-Boussinesq answers could change.

  When in non-Boussinesq mode, write out mass-based tracer inventories from
calls to MOM_tracer_chkinv and mean thicknesses in mass per unit area from calls
to MOM_state_stats, rather than the approximate volumes that depend on the
Boussinesq reference density.  Similarly the thicknesses and transports output
by write_u_accel and write_v_accel are in mass per unit area or mass flux per
unit face length in non-Boussinesq mode.  This is done specifically by using
GV%H_to_MKS in place of GV%H_to_m; these are identical in Boussinesq mode, but
GV%H_to_MKS avoids rescaling by the Boussinesq reference density in
non-Boussinesq mode.  Several comments were also updated to reflect these
changes.  All solutions are identical, but there are changes in the units of
several diagnostic output fields that the model can write out when debugging
non-Boussinesq mode runs.

  Added units in comments describing about 141 real variables in 7 framework
modules.  All answers are bitwise identical, and only comments are changed.

  Added or corrected the units in comments describing about 115 real variables
scattered across 14 tracer modules and 10 other modules.  All answers are
bitwise identical, and only comments are changed.

This patch removes the code coverage upload requirement.  Constraints
around codecov.io upload rules have made it impossible to keep this as a
requirement.  However, we will still attempt an upload, which should be
more successful for accounts with a stored URL token, such as NOAA-GFDL.

This patch adds a relatively robust parser for C preprocessor
expressions inside of an #if statement.

The following are supported:

* Nearly all operators, including arithmetic, logical, and bitwise,

* Parentheses within expressions,

* defined() evaluations.

The following are explicitly not supported:

* Function macros,

* Multiline preprocessors.

No doubt there are other lingering issues, but this is comprehensive
enough to handle both MOM6 as well as current and legacy FMS source
codes.

Existing Makefile.dep output files appear to be mostly unchanged.  One
rule (data_override.o) had its arguments reordered but is otherwise
unchanged.  mpp_data.o had its rule corrected to use mpp_util_mpi.inc
rather than mpp_util_nocomm.inc.

Some fixes and adjustments were made to the overall makedep source:

* Input macros (-D) are now stored as key-value dicts, rather than
  simply a list of macro names.

* Input macros are now passed to all scan_fortran_file() calls, rather
  than just the Fortran source.

* Input macros are now correctly passed to FMS makedep.  Previously,
  these were omitted from the Makefile generation.

* Previously, #if blocks were always set to True, even though the
  comments indicated that they were always set to False.  Given that
  neither of these was ever correct, it's amazing that we were able to
  survive this long without prior incident.

The motivation for this PR comes from issues with Makefile generation in
FMS.  Older versions of FMS were unable to correctly resolve their
dependencies in fft.f90 on certain systems (perhaps caused
by filesystem peculiarities).  Newer versions of FMS were unable to
handle the #if block default from True to False.  Inevitably, we threw
up our hands and solved the underlying problem.