* Change dumbbell initialization
* Change in Dumbbell Layer Mode
* Fix sponge diagnostics
* Fix ALE Sponge Diagnostics
* A minor style change removing spaces around = in optional. function arguments.
* Fix ALE sponge diagnostics
* character declaration fix

  - Ice shelf needs to be initialized prior to
    ocean state initialization. This fixes any cases with
    an ice shelf using the FMScap.

  Added the new runtime parameter CHANNEL_DRAG_MAX_BBL_THICK, to allow the
CHANNEL_DRAG parameterization to use the full dynamically determined depth of
the bottom boundary layer, or to use a fixed maximum thickness, regardless of
the settings of other parameters, although for now the default value is set
based on the settings of USE_JACKSON_PARAM and DRAG_AS_BODY_FORCE in order to
reproduce existing solutions by default.  There are new entries in some
MOM_parameter_doc files.  All answers in the MOM6-examples test suite are
bitwise identical, and this test suite has been tested and works with revised
values of this parameter.

* Porous barrier variables, which are grouped in porous_barrier_ptrs,
are changed from pointers to allocatables.
* Copies (aliases) of these variables in MOM_CS are removed.
* The name porous_barrier_ptrs is yet to be changed to minimize the
number of files needed to be edited.

* The interface porous_widths is deleted and subroutine por_widths is
renamed as porous_widths.
* porous_barrier_CS is added to control input and diagnostics of the
 module
* Diagnostics for both interface and layer averages weights are added in
subroutine porous_widths.
* An _init subroutine is added to facilitate reading parameters and
registering diagnostics
* checksum debugs are added within subroutine porous_widths.

This commit primarily fixes indexing bugs in subroutine porous_widths.

* Loop range is subroutine porous_widths is changed from data domain to
computation domain.
* find_eta call is now with a proper halo to cover all velocity points.
* Halo update for porous barrier fields is added in step_MOM_*.

Other changes:
* mask_depth, a component of porous_barrier_CS is now used to as
the criterion for masking. This removes the limit that the cell edge
depth has to be below sea surface.
* Output variable eta_cor was unused and is now changed to a local.
* Some unused indexing variables are removed.

parameters

* subroutine set_subgrid_topo_at_vel_from_file is added to read max, min
and avg depth at the edges from file.
* The subroutine is only called when a new runtime parameter
SUBGRID_TOPO_AT_VEL is True. Default is false.

* id_clock is added (as a private variable) to time porous barrier
calculations.
* Some small format fixes

* A runtime parameter PORBAR_ETA_INTERP is added to control different
methods for calculating interface height at the velocity points from
adjacent tracer cells.
* Two small thickness variables are added and scaled the unit of eta.
This is to assit the harmonic mean calculation, but eventually the
if-test eta_s - eta_prev>0 should be relaced by Angstrom.
* Runtime parameter POROUS_BARRIER_MASKING_DEPTH is renamed to make it
a liitlle bit shorter.
* log_version call is added to separate out porous_barriers module in
parameter file.

The main loops in porous_widths are simpified and cleaned up.
* calc_por_layer iss slightly modified to reduced the if-blocks in
the main loops.
* k-loops are moved out.
* To assist this new structure, two 2-D arrays are added to the stack.
* A new function eta_at_uv is added to treat different interpolations.
This is currently done at every point within the loop, and it is
probably adding too much an overhead. It is better pre-calculate
 eta_[uv] all together, which would increase the stack size.

The averaged weight over a layer (por_face_area[UV]) and the weight
at the interface (por_layer_width[UV]) are now calculated
separately, as the only two updates in step_mom are in fact
using only one of them. This reduces
some unnecessary calculations.

* Two new subroutines replaced the original `porous_widths`. They could
be combined with interface if we remove porous_barrier_ptrs type, and
use variables (of different nz) as output.
* There is a place holder switch do_next in the two calc_por subs. This
is used to further reduce calculations for all layers above D_max.
Will be tested and implemented in the next commit.
* A new subroutine calc_eta_at_uv is added to calculate eta at uv.
This simplifies the code, but also increases stack and some performance
overhead.

* A new runtime parameter USE_POROUS_BARRIER is added to control this
feature. The default is True at the moment to assist potential test. It
should be changed to false in the future.
* A boolean array is added to avoid unnecessary porous barrier
 calculations above the shallowest points.
* The layer-averaged weights are bounded by 1.0 explicitly, to avoid
the cases it goes beyond due to roundoff errors.
* For very thin layers (Angstrom), layer-averaged weights are set to
zeros for simplicity.
  * Perhaps it is more reasonable to use the inferace weight for these
  cases.
  * It might be useful to make the thin layer definition a runtime
  parameter.
* A bug related the size of eta_[uv] is fixed.

This commit is potentially answer-changing when the porouss barrier
module is used.

* The recently introduced ANSWER_DATE functionality is used to maintain a
version that should be bit-for-bit reproducible for previous porous barrier
related runs.
* Rename porous_barrier_ptrs to porous_barrier_type
* Some small documentation edits

* Fixed a bug that eta_[uv] was not properly initialized, which caused issues
with global chksums with DEBUG=True.
* Documents added to comply with Doxygen tests

* Fix nudged OBCs for tracers

Failed codecov.io uploads now report the stderr output.

This will allow us to start diagnosing the intermittent failures and,
possibly, find a solution.

Report error logs from codecov upload fail

  Added the new module MOM_interface_filter to allow for a biharmonic or
other order of filtering of the interface heights.  This new capability
is enabled with the new runtime parameter APPLY_INTERFACE_FILTER, and with
rates of filtering determined by the new runtime parameters
INTERFACE_FILTER_TIME, INTERFACE_FILTER_MAX_CFL, INTERFACE_FILTER_ORDER and
INTERFACE_FILTER_ISOTROPIC.  Set APPLY_INTERFACE_FILTER to True and provide a
positive timescale in INTERFACE_FILTER_TIME to enable the filtering.

  The comments describing THICKNESSDIFFUSE and THICKNESSDIFFUSE_FIRST were
also updated to clarify the distinctions with the new capabilities or to
clearly identify which parameters work together.

  By default, all answers are bitwise identical, but there are new entries or
modified comments in the MOM_parameter_doc files.

  Corrected the doxygen comments for the interface filtering module, to avoid
a section-name conflict with the GM module, and renamed some internal variables
for greater clarity when reading the code.  All answers are bitwise identical.

The sigsetjmp function is part of the POSIX, but is not required to be
defined as a symbol, and may be implemented as a macro.  Since Fortran
C bindings require a symbol, we cannot bind to macro implementations.

The prior implementation assumed a Linux glibc binding of __sigsetjmp
(accessed by a `sigsetjmp` macro), but this did not work on BSD and
MacOS builds, which have a dedicated `sigsetjmp` symbol.

Although the autoconf build included a macro to test and assign the
symbol to `SIGSETJMP_NAME`, this did not resolve builds based on mkmf or
similar build systems, and would fail to compile.

To resolve this, the SIGSETJMP_NAME points to a placeholder function,
`sigsetjmp_missing` which permits compilation but raises an error if
called.

Since this function is only used in our unit testing, and even then only
for tests which would otherwise raise FATAL, this change will not
disrupt any simulations.

However, it does mean that only "power" users who build with either
autoconf or `-DSIGSETJMP_NAME=\"...\"` will be able to run the unit
tests.  In practice, it should be sufficient to direct users to the
autoconf builds, and no actual disruptions are expected.

Disable sigsetjmp for default compilation

  Added a new argument to set_diffusivity_init and bkgnd_mixing_init to specify
when a physically based ocean surface boundary layer mixing scheme, like KPP,
ePBL or a bulk mixed layer is in use, and use this to restrict when a constant
diffusivity is used as a crude parameterization of the surface boundary layer.
This constant mixed layer diffusivity was formerly set with KDML, but is now set
with KD_ML_TOT, but KDML will still work for input but is no longer logged and
triggers a warning message if it is used.  Also removed KDML from the
description of the ADIABATIC option.  This leads to changes in the entries of
several MOM_parameter_doc files, and it could lead to answer changes that
set a non-default value of KDML and use ePBL or KPP, but it seems unlikely that
any such cases would exist.

  Replaced the runtime parameter KVBBL, which sets the total viscosity in the
bottom boundary layer when BOTTOMDRAGLAW is false, with KV_EXTRA_BBL, which sets
the viscosity anomaly relative to KV or other turbulent viscosities, and has a
default value of 0 m2/s, instead of the previous nonzero default for KVBBL.
Specifying KVBBL will still give equivalent results, but is no longer logged and
triggers a warning message if it is used.  This leads to changes in the entries
of several MOM_parameter_doc files, and while all existing solutions are bitwise
identical, by changing the definition of KVBBL there could be some changes at
the level of roundoff for cases that do not use BOTTOMDRAGLAW = True.

  Corrected a spelling error, from HARMOINIC to HARMONIC, in the get_param
description of one of the valid options for PORBAR_ETA_INTERP, and a few other
spelling errors in comments in MOM_porous_barriers.F90.  All answers are bitwise
identical, but an entry in some MOM_parameter_doc files will change.

  Renamed the runtime parameter KVML with KV_ML_INVZ2 and modified the comments
describing it to more clearly indicate what it does and to explicitly discourage
its use. Specifying KVML will still give identical results, but is no longer
logged and triggers a warning message if it is used.  The two descriptions of
the HBBL runtime parameter in MOM_vert_friction and MOM_set_viscosity were also
made consistent.  Also added units to several comments describing real variables
and corrected some spelling errors in comments in MOM_vert_friction.F90.  This
commit leads to changes in the entries of several MOM_parameter_doc files, but
all existing solutions are bitwise identical.

  Refactored the way that find_coupling_coef sets the viscous coupling in the
surface boundary layer, in preparation for adding new options.  All answers are
bitwise identical.

  Added the new runtime parameters FIXED_DEPTH_LOTW_ML and LOTW_VISCOUS_ML_FLOOR
to turn on a law-of-the-wall based scheme to specify the viscosity in a boundary
layer of thickness HMIX_FIXED, and to apply a law-of-the-wall based lower bound
on the viscous coupling strengths across interfaces that are within the surface
boundary layer thickness, as determined in one of several different ways.
Because these are false by default, all answers are bitwise identical, but there
are new entries in some MOM_parameter_doc files.

  Set the friction velocity, ustar, when the wind stresses are specified via
wind_forcing_2gyre and wind_forcing_2gyre.  This was done by adding the new
subroutine stresses_to_ustar to start to standardize the calculation of the
friction velocities from the wind stresses.  Neverworld_wind_forcing and
wind_forcing_gyres were already using the same expressions for ustar as are in
stresses_to_ustar, so the new routine is now used in these cases.  In other
routines that set ustar, the answers differ at roundoff due to the difference
between division by a reference density and multiplication by its reciprocal, so
using the new subroutine in those cases should be wrapped in an answer_date
flag, which is being left for a later commit.  This is being added now as a
result of a new commit that uses ustar as a part of the specification of
viscosities even in idealized cases.  All answers in the MOM6-examples test
suite are bitwise identical.

  Made the von Karman constant into a runtime parameter, specified with
VON_KARMAN_CONST or VON_KARMAN_BBL, replacing a hard-coded value in multiple
places.  By default, all answers are bitwise identical, but there are one or
two new entries in the MOM_parameter_doc files.

Obsolete the HENYEY_IGW_BACKGROUND_NEW option.

This patch modifies the calculation of barotropic weight functions
wt_[uv] to a more vectorizable form.  The if-blocks used to evaluate the
piecewise function was replaced with a sequence of min/max evaluations,
which are vectorizable.

A nondimensional subroundoff term is also introducted to address the
division by zero in the evaluation of 1 - 0.5 * Instep / visc_rem.
When visc_rem is zero, this term will never be selected by the max()
operator, since it is negatively asymptotic.  When visc_rem is
sufficiently large to accept this term, the value of subroundoff is too
small to modify the significand of the floating point value.

In our runs, the overall time of btstep was reduced by 6%.

Experiment: benchmark
Platform: AMD Ryzen 5 2600
Compiler: gfortran 12.2
Flags: -O3 -mavx

Timer: (Ocean BT pre-calcs only)

Before:
- 1.781952
- 1.787916
- 1.790447

After:
- 1.544491
- 1.553111
- 1.544491

Timer: (Ocean barotropic mode stepping)

Before:
- 4.042281
- 4.050225
- 4.051413

After:
- 3.785838
- 3.822393
- 3.835511

Speedup:
~13% of "(Ocean BT pre-calc step)"
~6% speedup of btstep

Barotropic weight vectorization

  Add calls to query_averaging_enabled to avoid unnecessary extra calls to
post_data for diagnostics that will never be used in vertvisc and vertvisc_coef.
All answers are bitwise identical, but the C.I. testing appears to fail because
there are fewer calls to post_data.

  Updated the default values for KV_ML_INVZ2, BAROTROPIC_TIDAL_SAL_BUG,
LAYER_Z_INIT_IC_EXTRAP_BUG, and KAPPA_SHEAR_VERTEX_PSURF_BUG. Also completed
obsoleting of HENYEY_IGW_BACKGROUND_NEW.  This PR will change answers in some
cases, unless the MOM_input files already included lines with the following
settings, or otherwise explicitly set these parameters:

  KV_ML_INVZ2 = 0.0
  BAROTROPIC_TIDAL_SAL_BUG = False
  LAYER_Z_INIT_IC_EXTRAP_BUG = False
  KAPPA_SHEAR_VERTEX_PSURF_BUG = False

For each of these parameters, there is another parameter (such as BULKMIXEDLAYER
or TIDES) whose value determines whether they might be used in a particular
case.  To determine whether or why answers for a particular configuration might
change with this PR, compare the MOM_parameter_doc.all files from equivalent
runs before and after these code changes to determine whether any of these
parameters are taking the default value.  This commit could change answers in
some cases that use default values for these parameters, and the entries in some
MOM_parameter_doc files will change.  The answers in the MOM6-examples test
suite are bitwise identical, due to some entries that were recently added to the
MOM_input files for these cases.

* Description of wt_eta related  variables
* Longname of PFu_BT

Logging BT halo size parameters should happen after the modification in
clone_MOM_domain, where they are lower-bounded by the domain halo.
This commit fixes this and therefore may result a change of parameters
in MOM_parameter_doc.layout in past runs.

This commit fixes the conflicting timestamps for diagnostics at the BT
step level (do_hifreq_output=True).
Instead of the actual dtbt, a nominal time interval is used to "squeeze"
output at all nstep+nfilter steps into a time interval of a single
baroclinic timestep.

This patch re-implements the FMS2 implementations of `open_ASCII_file`
and `open_namelist_file` to remove their dependency on FMS1 functions
which have been staged for deletion.

Note that if a file is opened with `mpp_open` but closed with
`close_file_unit`, then it will raise an error in `fms_io_exit`.
This will no longer be an issue after all references to `mpp_open` have
been removed.  But in the meantime, we will need to ensure that all
unit-based `close_file` calls were not opened with `mpp_open`.

There is also a minor patch to `.testing/Makefile` which selects the
framework ("infra") source dependency, rather than hard-set to FMS1.

FMS2: open_ASCII_file and open_namelist_file

  Added the new overloaded interface only_read_from_restarts to allow for
variables that are not registered to be stored in restart files to be read from
the restart files.  Versions of this routine that work on 2d, 3d or 4d arrays
and pairs of 3d arrays have been added for now, although the pattern should be
clear enough if other sizes of arrays are needed.  They use the new private
routine find_var_in_restart_files to determine whether a named variable is in
the restart files and the full path to the appropriate file. These routines
allow for the contents of the restart files to be gracefully altered, while
preserving the ability to use restart files with the previous format.  All
answers are bitwise identical, but there is a new public interface.

  Added the option to store the Coriolis and advective accelerations from the
end of one split RK2 dynamics time step for use in the predictor phase of the
next time step, in order to allow for a future commit that will properly remap
these accelerations in ALE mode.  This change is controlled by the new runtime
parameter STORE_CORIOLIS_ACCEL, and it leads to changes in the fields that are
saved in the restart files.   This commit also includes the necessary changes to
allow the restart files written in one mode to be read in the other and give
bitwise identical restarts.  Because the next step's predictor CAu and CAv
accelerations now have to coexist with the previous step's CAu and CAv
accelerations, which are used for derived diagnostics (like energy budget
terms), the predictor CAu and CAv terms have to be stored in new variables,
CS%CAu_pred and CS%CAv_pred, and there is a new accel_diag_ptrs type in the
dyn_split_RK2 control structure to diagnose truncation errors arising in the
predictor step.  Several comments have also been modified to document the
dimensions of variables.  All answers and diagnostics are bitwise identical.

Changes made to the main equations documentation files for MOM6. Changes
were made on top of Kate's original documentation.

The latest homebrew `netcdf` package appears to no longer include
netCDF Fortran headers and tools.  This patch installs the new
`netcdf-fortran` package which adds them back.

This is an initial step to add the functionality of calculating
self-attraction and loading (SAL) with the spherical harmonic transform
(SHT) method. The code in this commit is a direction translation and
adoption of the "parallel SAL" from MPAS-Ocean model developed at the
Department of Energy's Los Alamos National Laboratory
(Barton et al. (2022) and Brus et al. (2022)).

This commit serves as an "MVP" so that the MPAS-Ocean algorithm is
functional in MOM6. Further refactor, optimization and documentation
should be put in place.

* A new module MOM_spherical_harmonics is introduced to house
precomputed coefficients used for SHT.
* New subroutines are added in module MOM_tidal_forcing for calculating
SAL with the SHT method.

The actual calculation of SHT is stripped out of the tidal SAL
subroutine to make module MOM_spherical_harmonics self-contained.

* Forward and inverse spherical harmonic transform calculations are
imported from MOM_tidal_forcing module, as two new subroutines are in
module MOM_spherical_harmonics. This generalization allows the
spherical harmonic transforms to applications other than SSH based SAL.
* The main loops in the two new subroutine are also simplified.
* Child variables in the control structure for SHT are now private.
* The associated Legendre polynomials in the diagonal (n=m) is now
precomputed and stored, rather than recalculating at every step. This
does not have the memory cost as fully precomputing all the polynomials,
but still reduces some repeat calculations. The performance appears to
be improved with this change.
* SAL calculation is simplified as forward transform, scaling and
inverse transform.
* A few variables names inherited from MPAS-O are changed.

* A new module is introduced with the sole purpose of storing Love
numbers. This makes module MOM_tidal_forcing less chunkier and easier
to read.

* A new function in module MOM_spherical_harmonics is introduced to
calculate the starting index of the first element for each order m,
replacing a previous function of a similar purpose but a bit more
complicated.

* Option to do reproducing sums for forward spherical harmonic transform
is added, which is controlled by a runtime parameter.

* New timers are added to monitor the performance of the more expensive
global sum.

* Associated Legendre polynomials calculation is vectorized.

* Some redundant variables are removed.

* Fix filename suffix for module MOM_spherical_harmonics

* In MOM_tidal_forcing module, spherical harmonic coefficients (for SAL)
are now parts of tidal_forcing_CS to avoid repeated allocations. The
same applies to the Love number scaling factors.

* Allocations for arrays used for reproducing sums are moved to
subroutine spherical_harmonics_init in the MOM_spherical_harmonics
module.

* Documentations are added for all modules related to SHT SAL.
References from the MPAS-O group will need to be updated once they are
published.

* Variables names are shortened and changed to follow MOM6 style.
(from camel to snake)

* Change RhoE and RhoW in Love number scaling to runtime parameters

* Correct a bug in a_recur size

* Local arrays in SHT subroutines are properly initialized.

* Love number scaling coefficients multiplication is re-ordered to avoid
ambiguity and to follow MOM6 style.
* Coefficients for Pmm is simplified.

Both changes introduce bit-wise level answer change to previous SHT SAL
related unpublished commits.

* Patches for Doxygen

  Started the refactoring of code related to ALE_main, including:

 - Split remap_all_state_vars into remap_tracers and remap_velocities
 - Created the new public subroutine pre_ALE_diagnostics, and call it from
   step_MOM
 - Eliminate the unused argument conv_adjust to regridding_main
 - Eliminated check_grid, and moved these tests into regridding_main
 - Eliminate check_remapping_grid, and replace it with calls directly to
   check_grid_column inside of remapping_main

All answers are bitwise identical, but there are new public interfaces and
changes to the arguments of other public interfaces.

  Continued the refactoring of code related to ALE_main, including:

 - Added the new public subroutine pre_ALE_adjustments
 - Added the new arguments h_new, dzRegrid, and PCM_cell to ALE_main
 - Added the new arguments h_new and dzRegrid to ALE_offline_main
 - Moved the code copying the new thickness to the primary thickness out of
   ALE_main and into step_MOM_thermo

All answers are bitwise identical, but there are new public interfaces and
changes to the arguments of other public interfaces.

  Move the post_data call for dzRegrid to immediately follow its calculation in
regridding main.  All answers and diagnostics are bitwise identical, but because
this commit changes the order in which some diagnostics are posted, the TC
testing falsely reports a failure due to changed diagnostics.

  Replaced ALE_main with ALE_regrid and calls to ALE_remap_tracers and
ALE_remap_velocities from step_MOM_thermo.  Also eliminated ALE_main_offline and
ALE_offline_tracer_final, as they can now be replaced with a calls to ALE_regrid
and a call to ALE_remap_tracers.  Also added unit descriptions to a number of
comments describing variables in MOM_ALE.F90.  All answers are bitwise
identical, but there are 3 new publicly visible subroutines (ALE_regrid,
ALE_remap_tracers, and ALE_remap_velocities) and three previous interfaces
(ALE_main, ALE_main_offline and ALE_offline_tracer_final) have been eliminated.

  Revised or extended comments to correct or more clearly document the units of
49 internal variables in modules in the core directory.  All answers are bitwise
identical.

  Added separate restart variables for some of the oblique OBC restart variables
at north-south or east-west faces.  Before this fix, some of the full 3-d arrays
for restarts were being overwritten for oblique OBC segments that join at
adjacent corners on the north-east side of tracer points, as is noted in
github.com/NOAA-GFDL/issues/208.  The discussion surrounding this issue
confirms that there are cases using oblique OBCs (like some North-West Atlantic
cases) that do not past the restart reproducibility tests.  Some halo updates
were also corrected, in accordance with the introduction of these new variables
and their proper staggering locations.  In a number of places, the proper
case-sensitive horizontal indexing convention, as described in
github.com/NOAA-GFDL/MOM6/wiki/Code-style-guide#soft-case-convention, is now
being used.  This commit also corrected the comments describing a number of the
variables in the OBC_segment_type to make it clearer where they are discretized.

  This changes the names of oblique OBC-related variables in the restart files,
but since the previous version did not reproduce across restarts, there does not
seem to be any point in retaining those incorrect answers.  All answers in the
MOM6-examples test suite are bitwise identical, but this will change (fix)
solutions with oblique OBCs and OBC_RAD_VEL_WT < 1.

  Added two new arrays (OBCmaskCu and OBCmaskCv) to the ocean_grid_type and
dyn_horgrid_type to mask out values over land or at open boundary condition
points.  Without open boundary conditions, these arrays are identical to
mask2dCu and mask2dCv.  These arrays are used in some of the lateral
parameterization modules to zero out certain gradient-dependent fluxes at open
boundary points.  With these changes, the Bering test case solutions no longer
exhibit any dependence on whether DEBUG_OBC is true or false or the value of
OBC_SILLY_THICK, so this commit should help to address some of the issues
discussed in github.com/NOAA-GFDL/issues/208.  All answers are bitwise
identical in the MOM6-examples test cases, but they change for some other tests
that use open boundary conditions more extensively, and there are new arrays in
some transparent types.

  Eliminated OBC arguments that are no longer used by three internal routines in
MOM_lateral_mixing_coeffs.  All answers are bitwise identical.

This patch moves the local `use_tides` flag of the split RK2 solver,
which tracks the TIDES parameter, into the RK2 control structure (CS),
and uses it to conditionally call `tidal_forcing_end().

The `tidal_forcing_init()` function conditionally allocates the
spherical harmonic fields, but `tidal_forcing_end()` is always called,
which will in turn always attempt to deallocate the fields.  In some
compilers, under certain levels of debugging, this can raise a runtime
error.

…ain-2022-10-27

GFDL to main 2022-10-27

…20221112