Added parentheses to expressions taking the squares of the wind stress
components in 70 lines in 7 files so that these expressions will be rotationally
invariant when fused-multiply-adds are enabled.  All answers are bitwise
identical in cases without FMAs, but answers could change with FMAs.

  Added parentheses to 9 expressions like `curv_3 = h_W(i) + h_E(i) - 2.0*h(i)`
in PPM_limit_pos, zonal_flux_layer, zonal_flux_thickness, merid_flux_layer and
merid_flux_thickness, changing them to `curv_3 = (h_W(i) + h_E(i)) - 2.0*h(i)`.
This change enforces the order of arithmetic that is required to give rotational
symmetry, but it also is the order that the Intel, GNU, and Nvidia compliers
were all already using in these expressions.  Moreover, had the order of
arithmetic ever been anything else, this would have led to failures in our
rotational consistency and redundant point consistency testing, and almost
certainly would have been detected before. However, by adding these parentheses,
there is a remote chance that the addition of these parentheses could change
answers for some compiler or compiler settings we have never tested before.
This change should not impact any FMA-enabled calculations.  All answers are
bitwise identical in the MOM6-examples regression suite as run on Gaea.

  Added parentheses to 16 expressions setting the grad_gradient arrays with
oblique_grad open boundary conditions and setting cff_new with all kinds of
oblique boundary conditions so that they will be rotationally invariant when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers with certain types of open boundary conditions could
change with FMAs.

  Added parentheses to 150 lines in the 5 generic density integral routines
(int_density_dz_generic_pcm, int_density_dz_generic_plm,
int_density_dz_generic_ppm, int_spec_vol_dp_generic_pcm and
int_spec_vol_dp_generic_plm) in the MOM_density_integrals module so that they
will be rotationally invariant when fused-multiply-adds are enabled.  All
answers are bitwise identical in cases without FMAs, but answers could change
with FMAs.

  Added parentheses to 140 lines in 8 int_density_dz and int_spec_vol_dp
routines for the linear, Wright, Wright_full and Wright_red equations of state
so that they will be rotationally invariant when fused-multiply-adds are
enabled.  All answers are bitwise identical in cases without FMAs, but answers
could change with FMAs.

  Removed recently added parentheses around expressions like '+ (hL*hR)' in 110
lines in MOM_density_integrals and 4 equation of state module to reflect that
these parentheses are not necessary for rotational symmetry in
FMAs.  All answers are bitwise identical in cases without FMAs, but
answers could change with FMAs.

  Added parentheses to 4 lines in PressureForce_Mont_nonBouss and
PressureForce_Mont_Bouss so that they will be rotationally invariant when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change with FMAs in cases that use the
Montgomery potential form of the pressure gradient accelerations.

  Added parentheses to 4 lines in calc_isoneutral_slopes so that they will be rotationally invariant when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change with FMAs.

  Added parentheses to 2 lines in MOM_calc_varT so that they will be
rotationally invariant when fused-multiply-adds are enabled.  In this case, FMAs
can still be applied to the impacted lines, exploiting that the masks are always
0 or 1.  Also added parentheses to 2 other lines used to generate the stochastic
pattern for rotational symmetry with FMAs.  All answers are bitwise identical in
cases without FMAs, but answers could change with FMAs.

  Added parentheses to the calculation of the diffusive temperature changes in
tracer_hordiff so that it will be rotationally invariant when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change with FMAs.

  Added parentheses to the calculation of the iceberg contribution to the
fractional area of ice shelves in iceberg_forces so that it will be rotationally
invariant when fused-multiply-adds are enabled.  All answers are bitwise
identical in cases without FMAs, but answers could change with FMAs enabled in
cases with tabular icebergs.

  Added parentheses to the calculation of the Stokes-drift Coriolis velocity
increments in CoriolisStokes so that it will be rotationally invariant when
fused-multiply-adds are enabled.  All answers are bitwise identical because
CoriolisStokes is still under development and is never called, with a fatal
error occurring if anyone tries to use it.  Also added parentheses to two
expressions calculating the magnitude of the Stokes velocity in
get_Langmuir_Number.  Answers could change for some cases that use Langmuir
turbulence parameterizations with FMAs enabled.

  Added the new element Coriolis2Bu to the ocean_grid_type and the
dyn_horgrid_type to hold the square of the Coriolis parameter, and use this
array in 10 routines (including btstep, set_dtbt, calculate_diagnostic_fields,
VarMix_init, propagate_int_tide, Calculate_kappa_shear, Calc_kappa_shear_vertex
and add_MLrad_diffusivity) that had been calculating and averaging the square of
the Coriolis parameter.  This could change some answers with FMAs enabled
because the compilers were previously free to split up some of the squares
when averaging the squared Coriolis parameter, but without FMAs all answers are
bitwise identical.  This commit does add a new element to two transparent
types.

  Added parentheses to 17 expressions in thickness_diffuse_full,
thickness_diffuse and thickness_diffuse_init to give rotationally consistent
solutions when fused-multiply-adds are enabled.  One comment was also added to
note that the calculation of PE_release_h is does not exhibit rotational
symmetry when MEKE_GM_SRC_ALT is set to true.  All answers are bitwise identical
in cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 2 expressions in the Zanna_Bolton code and rearranged
another line so that the u- and v-discretizations introduce terms in the same
order so that the Zanna_Bolton code will exhibit rotationally consistent
solutions when fused-multiply-adds are enabled.  All answers are bitwise
identical in cases without FMAs, but answers could change with FMAs enabled in
cases that use the Zanna-Bolton parameterization.

  Added parentheses to 19 expressions in the MOM_internal_tides propagation code
to exhibit rotationally consistent solutions when fused-multiply-adds are
enabled.  All answers are bitwise identical in cases without FMAs, but answers
could change when FMAs are enabled in models that use the ray-tracing based
internal tides code.

  Added parentheses to 10 expressions in find_uv_at_h to exhibit rotationally
consistent solutions and treat the velocities at both edges of a tracer cell
equivalently when fused-multiply-adds are enabled. All answers are bitwise
identical in cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 19 expressions in set_viscous_ML, set_u_at_v and
set_v_at_u to treat the velocities at both edges of a tracer cell equivalently
when fused-multiply-adds are enabled, and thereby to exhibit exhibit
rotationally consistent solutions.  Also swapped the order of the u- and
v-components in the u-point calculation of Uh2 to mirror the order of the
corresponging v-point calculation for the same purpose.  All answers are bitwise
identical in cases without FMAs, but answers could change when FMAs are enabled.

  Added mathematically equivalent rearrangements of the code in
calc_kappa_shear_vertex that interpolates velocities, temperatures and
salinities to the vertices to expose the mask variables while ensuring that the
other multiplications occur within parentheses so that they will exhibit
rotational symmetry when fused-multiply-adds are enabled.  FMAs can still occur,
but it will be multiplication by the 0-or-1 masks that are fused with an
addition.  Also added parentheses to 3 expressions calculating the squared shear
in calculate_projected_state for rotational symmetry with FMAs.  All answers are
bitwise identical in cases without FMAs, but answers could change when FMAs are
enabled.

  Added parentheses to 4 expressions in add_drag_diffusivity, set_BBL_TKE and
add_LOTW_BBL_diffusivity setting the bottom-drag contributions to TKE and
friction velocity so that they will exhibit rotationally consistent solutions
when fused-multiply-adds are enabled.  All answers are bitwise identical in
cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 20 expressions in CorAdCalc and one in gradKE to exhibit
rotationally consistent solutions when fused-multiply-adds are enabled.  All
answers are bitwise identical in cases without FMAs, but answers could change
when FMAs are enabled.

  Added parentheses to 18 expressions in btstep, and one more each in set_dtbt
and barotropic_init to exhibit rotationally consistent solutions when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 19 expressions in 5 routines (calc_Visbeck_coeffs_old,
calc_Eady_growth_rate_2D, calc_slope_functions_using_just_e,
calc_QG_Leith_viscosity VarMix_init) in MOM_lateral_mixing_coeffs.F90 to give
rotationally consistent solutions when fused-multiply-adds are enabled.  Also
reordered terms in a sum in the calculation of beta_dx2_u to mirror that of
beta_dx2_v, also for rotational symmetry with FMAs.  All answers are bitwise
identical in cases without FMAs, but answers could change for some parameter
settings when FMAs are enabled.

  Added parentheses to 40 expressions horizontal_viscosity and another 14
expressions in in hor_visc_init and 3 more in align_aniso_tensor_to_grid to give
rotationally consistent solutions when fused-multiply-adds are enabled.   Also
swapped the order of two terms in the expression for Del2u to mirror the order
of the corresponding terms in Del2v for rotational symmetry with FMAs.  All
answers are bitwise identical in cases without FMAs, but answers could change
when FMAs are enabled.

  Added parentheses to 20 sums of squares of x- and y- distances or velocity
components used for initialization in 8 modules to give rotationally consistent
solutions when fused-multiply-adds are enabled.  All answers are bitwise
identical in cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 29 sums of squares of velocity or other vector components
used in parameterizations in 9 modules to give rotationally consistent solutions
when fused-multiply-adds are enabled.  All answers are bitwise identical in
cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 9 diagnostics of Coriolis accelerations or expressions
used in the kinetic energy budgets to give rotationally consistent solutions
when fused-multiply-adds are enabled.  All answers are bitwise identical in
cases without FMAs, but answers could change when FMAs are enabled.

  Added parentheses to 4 tracer edge value calculations used with PPM tracer
advection to give rotationally consistent solutions when fused-multiply-adds are
enabled.  Although these lines may not appear to need parentheses, some
compliers appear to be putting these expressions directly into others, where the
direction of the flow seems to determine which multiplications are incorporated
into FMAs.  All answers are bitwise identical in cases without FMAs, but answers
could change when FMAs are enabled.

  Added parentheses around the full expressions for intz and intp in 10 more
lines in 4 generic density integral routines (int_density_dz_generic_pcm,
int_density_dz_generic_ppm, int_spec_vol_dp_generic_pcm and
int_spec_vol_dp_generic_plm) in the MOM_density_integrals module so that
non-Boussinesq cases will be rotationally invariant when fused-multiply-adds are
enabled.  Although this might not seem to do anything, these parentheses do
matter if these expressions are in-lined in the sums where they are used a few
lines later.  The analogous parentheses had previously been added to
int_density_dz_generic_plm.  All answers are bitwise identical in cases without
FMAs, but answers could change with FMAs.

  Added parentheses to prevent FMAs in 4 expressions in end_value_h4 that rely
on exact vertical symmetry in order to get the cancellations that are necessary
to pass the vertical remapping unit testing.  Before this change, the MOM6
unit-testing was failing when FMAs are enabled, but with it the unit-testing is
passing even when FMAs are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change with FMAs.

…ymmetry_main

(*)Add parentheses for FMA rotational symmetry

This patch makes some changes to makedep output:

* Compile commands now end with the source code, rather than the include
  files.  This makes it easier to see which file is being compiled.

* Blank lines were added inbetween object file rules.

* Redundant spaces from absent include flags is now removed.

Also some minor PEP8-like cleanups:

* Some (but not all) 79+ character lines were reduced or split.

* Some (but not all) single-letter variables were renamed.

* A bad error handler for missing macros was fixed.

Nothing particularly major here, but build output should be more
readable.

  This merge updates dev/gfdl following PR mom-ocean#1634 to main that added parentheses
to give rotational symmetry when fused-multiply-adds are enabled.

Set default HAready = False

Also made additional adjustment to ensure harmonic analysis is not
activated if tide is not used in the model.

* (*)bugfixes for internal tides

- TKE_itidal_input values are set to zero where it is not applied to
  the model so that diagnostics are consistent

- removed useless halo updates

- diagnose energy loss terms as dE/dt instead of equivalent loss rate
  to properly close energy budget

- tot_vel_btTide2 was already a velocity squared, no need to square again

- Froude loss term was not properly initialized

- add missing halo updates for internal tide speed

- compute residual/critical slopes loss only where it has a meaning,
  allows to clean up noise in field

- uh/vh in energy flux routines do not need to be inout

* Refactor debugging internal tides

- put test for negative energt behind debug flag

- do energy budget in debug mode

- add flags to skip refraction, propagation for debugging

- add option to input TKE only at first step for debugging

- add checksums for unit scaling testing

- rewrite terms in refraction to avoid substractions of velocities

- rewrite gradient of coriolis in rotationally invariant form

* Add internal tides diffusivities

- MOM_internal_tides gets new routine that converts
  TKE losses into diffusivities using the TKE_to_Kd array

- MOM_set_diffusivities handles the diagnostics

* extend find_rho_bottom to get bbl thickness/index

* new calculation of BBL, profiles in H units

In the MOM_streaming_filter module, a system of two coupled ODEs is
configured as a streaming band-pass filter that takes the broadband
model output as the input and returns the filtered signal at a user
specified tidal frequency as the output. It is capable of detecting
the instantaneous tidal signals while the model is running, and can
be used to impose frequency-dependent wave drag or to de-tide model
output. An example of filtering the M2 and K1 barotropic velocities
is provided in the MOM_barotropic module.

Added runtime flags for the M2 and K1 filters. Simplified the code
by removing the control structure of the linked list, and now each
filter has its own control structure.

Using hor_index_type argument to avoid calculation in halo regions.

  Added the option of including the atmospheric or ice pressure and sea surface
height displacements from the global reference height in the pressures used in
the density calculations for Boussinesq pressure gradient calculations.  Note
that the full pressures were already being used everywhere apart from the
calculation of the equation of state.  This capability is controlled by the new
runtime parameter SSH_IN_EOS_PRESSURE_FOR_PGF.  This commit changes the Z_0p
argument to int_density_dz and 8 other routines (int_density_dz_generic_pcm,
int_density_dz_generic_plm, int_density_dz_generic_ppm, analytic_int_density_dz,
int_density_dz_wright, int_density_dz_wright_full, int_density_dz_wright_red
and int_density_dz_linear) from scalars into 2-d arrays, as were the internal
z0pres arrays in most of these routines.  By default, all answers are bitwise
identical, but there is a new runtime parameter in the MOM_parameter_doc files
for Boussinesq cases.

  Rearranged the calculation of the layer average pressure used in the density
integrals to have one less multiply and be clearer to read.  This could change
answers at roundoff if REFERENCE_HEIGHT is nonzero (it is 0 by default and in
all known cases) or if SSH_IN_EOS_PRESSURE_FOR_PGF is set to true.  Because the
former is mostly used for approximate self-consistency testing and the latter
option has just been added, it is unlikely that answers will change for any
production runs.

  Added the ability to rotate the ice shelf forcing and initialization.  The
specific changes include correcting the grid passed to a call to
initialize_ice_shelf in initialize_MOM, using a temporary mech_forcing type in
add_shelf_forces when the grid is rotated, uncommenting and correcting the draft
rotate_index code in initialize_ice_shelf, and correcting some error-checking
and the grid passed to add_shelf_fluxes in initialize_ice_shelf_forces.  In
addition, unused hor_index_type elements were removed from the ice_shelf_CS, and
some callTree calls were added to save_MOM_restart and initialize_ice_shelf.
Without rotation, all answers are bitwise identical, but some cases that would
have failed previously will now work.

  Add MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP option for the top boundary, using
top interface height, analogously to what is done near the bathymetry when
MASS_WEIGHT_IN_PRESSURE_GRADIENT.  The information from the is parameter is
encoded in the MassWghtInterp value passed to the various int_density_... and
int_spec_vol_... routines, so the routine interfaces do not change.  For now
this new bit of information about whether to do mass weighting near the surface
under ice shelves is only used in int_density_dz_generic_plm.  By default this
new option is set to false and no answers are changed.

  Implemented the MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP related capabilities to
the 13 density or specific volume integral routines that did not have them yet,
mirroring what was already done in int_density_dz_generic_plm.  This includes
both the density integral routines that are used primarily for Boussinesq cases
and the specific volume integral routines that are primarily used with
non-Boussinesq simulations.

  This change includes the addition of an optional SSH argument to
int_density_dz and 6 other related routines (int_density_dz_generic_pcm,
analytic_int_density_dz, int_density_dz_wright, int_density_dz_wright_full,
int_density_dz_wright_red and int_density_dz_linear).  It also includes the
addition of an optional P_surf argument to int_specific_vol_dp and 7 other
related routines (int_spec_vol_dp_generic_pcm, int_spec_vol_dp_generic_plm,
analytic_int_specific_vol_dp int_spec_vol_dp_wright,
int_spec_vol_dp_wright_full, int_spec_vol_dp_wright_red and
int_spec_vol_dp_linear).  Both of these new optional arguments are required
when the new near-surface mass weighting is activated, and there are test that
would issue a fatal error if they are not provided in such cases.  (Note that
these routines can be called from other places than the pressure gradient force
calculation, in which case the calculations that need these fields are not
done.)

  The diagnose_mass_weight_Z and diagnose_mass_weight_p diagnostic routines were
similarly revised to mirror the expanded range of capabilities in the integral
routines.

  This commit can change answers when MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP is
true, but all answers are bitwise identical otherwise.  There are new arguments
to 15 publicly visible routines.

  Add CORRECTION_INTXPA which makes a quadratic correction to surface pressure
integrals (intx_pa and inty_pa) under ice based on the horizontal gradients of
the in-situ density  anomaly along the surface.  The non-Boussinesq version
corrects the topmost value of intx_za and inty_za and uses in-situ specific
volume gradients along the ocean surface.  By default new the runtime parameter
CORRECTION_INTXPA is false, and answers are unchanged.

  Add CORRECTION_INTXPA_5PT which uses 5 point quadrature to calculate surface
pressure integral and therefore could work with a nonlinear EOS, or (if added)
different subgrid distributions.  This option requires CORRECTION_INTXPA = True.
By default CORRECTION_INTXPA_5PT is false and no answers are changed.

  Add RESET_INTXPA_INTEGRAL which resets intxpa and intypa at a trusted cell in
the interior (non-vanished and non-MWIPG-affected), and then integrates both up
and down to update intxpa and intypa for the interfaces above and below.  This
also adds the new runtime parameter RESET_INTXPA_H_NONVANISHED and to determine
when a cell is trusted.  This option is recommended with
MASS_WEIGHT_IN_PRESSURE_GRADIENT_IS for a quiet zstar ice shelf.  By default,
this option is not on and no answers change, but there are new parameters in
some MOM_parameter_doc files.

  Refactored and revised the recently added code in MOM_PressureForce_FV.F90 to
reduce the number of calls to the equation of state routines, and corrected a
number of minor bugs in the original implementation.  Answers are bitwise
identical unless the new options to reset the pressure gradient calculations are
actively selected.

  Refactored the CORRECTION_INTX_PA calculations to avoid multiple intermediate
steps, while also adding comments documenting the derivation of the final
expression.  Also calculate the pressures used in the equation of state calls
with the Boussinesq CORRECTION_INTX_PA and RESET_INTXPA_INTEGRAL options
consistently with the other terms.  These changes will change the answers when
either of those options are in use, but are bitwise identical when they are not.

  Added non-Boussinesq versions of the code that is exercised when
CORRECTION_INTXPA_5PT and RESET_INTXPA_INTEGRAL are set to true.  These options
use the same names as in their Boussinesq forms, even though the arrays that are
being adjusted are actually intx_za and inty_za.  As a part of the testing of
this commit, several checksums were added to the PressureForce_FV routines;
these are enabled when DEBUG = True.  The changes in this commit will change the
answers in non-Boussinesq cases when either of those options are in use, but are
bitwise identical when they are not, and all Boussinesq answers are bitwise
identical.

  Added code to do 5-point quadrature integrals when RESET_INTXPA_INTEGRAL and
CORRECTION_INTXPA_5PT are true.  Also extended the ranged of interfaces that
can be used for the corrections to include the bottom interface and use the
surface interface for the nonlinear pressure gradient force corrections either
when the ocean surface interface is within the pressure or height range of the
top cell, or when no appropriate interior interface has been found (for lack of
a better idea).  This latter case should probably be reconsidered later.  This
commit will change answers if RESET_INTXPA_INTEGRAL is true and one description
of a runtime parameter in the MOM_parameter_doc files has been revised.

  Eliminated the runtime parameter CORRECTION_INTXPA_5PT by always selecting the
code that would have been used if it were true.  This decision was based on the
relative performance of the options with this set to true or false, and on the
desire to maintain consistency with the 5-point Boole's rule quadrature used
elsewhere in the pressure gradient force calculations.  CORRECTION_INTXPA_5PT
was only added in the same PR as this commit, so it can been simply removed and
there is no need to obsolete it.  This will change answers if CORRECTION_INTXPA
or RESET_INTXPA_INTEGRAL are true and CORRECTION_INTXPA_5PT was set to false,
but otherwise answers are bitwise identical.  There are fewer entries in the
MOM_parameter_doc files as a result of this commit.

  Set MASS_WEIGHT_IN_PRESSURE_GRADIENT = True in the MOM_input files for the
tc1, tc2 and tc4 test cases in .testing, to permit the code to pass this testing
while dealing with inconsistent settings for two diagnostics related to the
mass-weighting in the pressure gradient forces.  This commit redefines these 3
test cases.

Adding null read operations so that compilers will not warn about unused
input variables in dummy functions.

This is a minor patch which allows the code to be compiled under F2023
standardization.

* Arguments are declared in the order that they are used.  In this case,
  it means that array lengths are declared before the arrays using them.

* The syntax of various IO statements has been cleaned up.

Some functions in the FMS cap used arguments which depended on other
arguments, which were declared out of order.

* ocean_model_data2D_get
* ocean_model_get_UV_surf

This patch moves the array index size definitions before the array
definitions.

This is required for language standard compliance.

This patch enables Fortran 2018 standard compliance testing.  We could
do 2023, but our current CI of choice doesn't yet have a compiler which
can do this.

The actual content of this PR is a decoupling of FCFLAGS_DEBUG and
FCFLAGS_FMS.  There is now a default FCFLAGS macro which is used by the
other two macros.  One can now optionally configure FCFLAGS_DEBUG
without worrying about the impact on FCFLAGS_FMS.

The motivation here is that we don't want to test for F2018/2023
compliance in FMS.

Temporarily cut off the control of VISC_REM_CONT_HVEL_FIX from
VISC_REM_BUG and change the default of VISC_REM_CONT_HVEL_FIX to False.

…l ice-sheet velocities (NOAA-GFDL#714)

* Updated with dev/gfdl. Then, added parameters for min ice thickness, min basal traction, max surf slope, and min ice viscosity to use for ice dynamics

* noted how the MIN_BASAL_TRACTION parameter input is in units [Pa m-1 yr], but converts automatically to [Pa m-1 s] in the code

* Added separate scalar ice-shelf diagnostics for Antarctica and Greenland

* Added sx_shelf and sy_shelf as ice shelf diagnostics to save the surface slope fields used in the shallow shelf approximation (SSA). This is particularly helpful for determining whether unrealistic velocities are caused by unrealistically steep surface slopes, which can sometimes arise for example, on coarse grid cells that cover both a steep mountainous region and a realively flat ice shelf. Then, the MAX_SURFACE_SLOPE parameter can be tuned to set an upper bound on the SSA surface slope to avoid these steep-slope-induced problematic velocities.

* Fix doxygen errors for ice-sheet process_and_post_scalar_data routine

* FMA fix associated with enforcement of max allowed ice-shelf surface slope