Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GFDL to main 2024-05-16 #1627

Closed
wants to merge 107 commits into from

Conversation

marshallward
Copy link
Collaborator

Once again, GFDL is submitting a truckload of PRs to the main branch. There are several new features and bugfixes, but they are diverse enough that I should probably let them speak for themselves.

With regards to testing, the intermittent failures in documentation, code coverage and (so-called) performance testing ought to be fixed now. 🤞

Thanks to @Hallberg-NOAA for organizing this batch (not to mention the dozens of contributions!)

Features

Bugfixes

Non-Boussinesq fixes (+Boussinesq refactoring)

Diagnostics

Input parameter handling

Refactoring

Variable documentation

Build system

Code maintenance

Contributors

adcroft and others added 30 commits November 22, 2023 17:41
- 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)
  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.
Hallberg-NOAA and others added 24 commits April 10, 2024 13:16
  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.
  Eliminated unneeded rescaling factors when setting the surface Stokes drift
spectrum variables (UStk_Hb and VStk_Hb) from forces%[UV]Stkb in
Update_Surface_Waves.  This bug appears to have been introduced when merging in
code that updates the Stokes drifts that was added on one branch (dev/ncar)
separately from changes in how these Stokes drifts are scaled in the
mech_forcing type on anther branch (dev/gfdl), although it might just have been
an oversight during refactoring.  All answers are bitwise identical when no
dimensional rescaling is being used, but answers could change for some uncommon
cases using the Stokes drift spectrum when dimensional rescaling is applied.
Setting `BBL_USE_TIDAL_BG` in OM5_b003 revealed non-conservation and
non-reproducibility across layouts. A scalar variable was correctly
being set based on the tidal amplitude but was being used after a break
between loops.

- Replaced the scalar variable U_bg_sq with a vector u2_bg(:) that
  is set to either a constant or the square of the tidal amplitude.
- The module parameter CS%drag_bg_vel was not set if using the tidal
  amplitude but WAS being used for conditionals. We now set it to
  1e30 because it should NOT be used or impact the solution. Setting
  it to zero, or anything meaningful would allow code to use it
  inadvertently. This "constant" does not need to satisfy any
  dimensional testing.
Print the line and its tokenization in parse_perf.py if the tokenization
fails for any reason.  This is most likely to occur if the fourth token
is not an integer.

This is possibly a platform-dependent issue, and not something easily
replicated on any machine we are using, so it is easier to simply
integrate it into the codebase.
  Added descriptions or added or corrected the units in the descriptions of 72
variables in the Neverworld_initialization and basin_builder modules.  This
includes minor refactoring to avoid reusing scalar variables with completely
different units.  Also refactored tidal_bay_set_OBC_data by moving the
multiplication by rescaling factors at one point so that the units of a variable
agree with their documented units.  All answers are bitwise identical.
- Adds maximum latitude above which to reduce the background diffusion
  to the minimum value.
  Added standard-format unit descriptions for 34 real variables in comments in
the subroutine vertPFmix.  Another 6 comments were added noting dimensionally
inconsistent expressions in vertFPmix.  Only comments are changed and all
answers are bitwise identical.
  Refactored 4 routines (int_density_generic_pcm, int_density_generic_ppm,
int_spec_vol_generic_pcm and int_spec_vol_generic_plm) in density integrals for
greater computational efficiency by doing fewer calls to the equation of state
routines but calculating entire rows of densities at subgrid locations with each
each call, replicating what was already being done int_density_dz_generic_plm.
To accomplish this, a number of variables now use larger arrays than
previously.  The total computational cost of the non-Boussinesq pressure
gradient force calculation was more than 50% greater with the previous code in
some tests.  All answers are bitwise identical.
…n and in model run

Adds INIT_BOUNDARY_EXTRAP parameter and function, so that REMAP_BOUNDARY_EXTRAP can just be used in the initialisation and not in the model dynamics, as required for an ice shelf in ALE mode.

Background:
When an ice shelf is initialised in ALE mode, MOM_state_initialization first initialises the ocean thickness and T/S with topography, but without an ice shelf. It then calls the function trim_for_ice, which takes in input of the ice shelf pressure, and in turn calls cut_off_column_top which truncates the columns (and adds vanishing layers) to account for the pressure of the ice depressing the free surface. The output of this is a grid with the ice shelf boundary treated in a ALE z-coordinate fashion with vanishing layers at the surface, and @adcroft's upcoming changes will make T and S consistent to the grid. However, after adding the ice shelf, we still need to regrid and remap onto our desired ALE coordinate, for example sigma coordinates. This can be done with regrid_accelerate in MOM_state_initialization or by the block controlled by REMAP_AFTER_INITIALIZATION in MOM.F90

We want the initialisation to use boundary extrapolation in the remapping so that the boundary cells preserve the desired initialisation profile when remapped to their target coordinate. This is controlled by REMAP_BOUNDARY_EXTRAP, except this runtime parameter modifies the control structure for the dynamic part of the model run too, which is bad because it means extrema are not preserved in the ALE remapping and can lead to non-sensible values of tracers.

Changes:

* Added runtime parameter INIT_BOUNDARY_EXTRAP and set the ALE control structure to use that instead of REMAP_BOUNDARY_EXTRAP at first. The default of INIT_BOUNDARY_EXTRAP is REMAP_BOUNDARY_EXTRAP to preserve answers (though it is not recommended to use REMAP_BOUNDARY_EXTRAP during model dynamics...)
 * After ALE regrid/remap in MOM.F90 (final initialisation step using ALE), reset the boundary extrapolation flag in the ALE control structure to go back to REMAP_BOUNDARY_EXTRAP.

This PR combined with @adcroft's upcoming changes should make the initialisation of the ocean thickness and TS in an ALE ice shelf cavity perfect, when the new flag INIT_BOUNDARY_EXTRAP = True and old flag REMAP_BOUNDARY_EXTRAP = False.
  Added the option of using rotationally symmetric expressions in
KPP_smooth_BLD().  This is enabled by setting the new runtime parameter
KPP%ANSWER_DATE to a value of 20240501 or higher.  For now the default is to
use the older symmetry-breaking expressions, but the default value for
KPP%ANSWER_DATE should be revised to use DEFAULT_ANSWER_DATE.  By default, all
answers are bitwise identical but there is a new entry in some
MOM_parameter_doc files.
  Refactored 10 diagnostics related to terms in the kinetic energy budgets and
2 diagnostics of nonlinear relative-vorticity related accelerations to respect
rotational symmetry.  These diagnostics are mathematically equivalent but will
change at roundoff due to changes in the order of arithmetic.  Only diagnostics
are changed, and the solutions themselves are bitwise identical in all cases.
  Fix several bugs when MEKE_GM_SRC_ALT is True, including corrections to two
dimensional rescaling factors, and optionally fix a bug that sets a limit only
on positive slopes but leaving negative slopes unlimited.  This bug is corrected
when the new runtime parameter MEKE_GM_SRC_ALT_SLOPE_BUG is false.  Additionally
there is a new runtime parameter MEKE_GM_SRC_ANSWER_DATE that specifies the use
of rotationally symmetric expressions for PE_release_h when it is set to
20240601 or higher, but it should be noted that rotational symmetry also
requires that MEKE_GM_SRC_ALT_SLOPE_BUG is false.  Four new checksum calls were
also added to verify the correctness of the calculation of MEKE%GM_src when
MEKE_GM_SRC_ALT and DEBUG are true.  By default, all answers are bitwise
identical but there are two new runtime parameters in some MOM_parameter_doc
files.
Vector u2_bg(:) was introduced in commit 6c44c5f (Fix for using
tidal amplitude in determining the BBL thickness, 2024-04-04). It
replaces a scalar in set_viscous_ML, but was not defined, causing
crashes. This commit sets the values of u2_bg(:) to match those given
given in set_viscous_BBL.
The makedep tool was updated to handle parentheses in preprocessor
expressions.

The expression `#if (defined a)` could not be parsed due to poor ad-hoc
handling of parentheses, making it impossible to build the AM2-based
coupled models with makedep,

The tool has has been significantly overhauled to include better overall
unary operator support.

* `defined` is now more of an operation than an exception, since it is
   pushed to the stack like any other operator.

* Parentheses are now handled as "operators", with `)` triggering an
  operator stack push and `(` conducting the actual operation (in this
  case simply returning the contents).

* In order to handle the possibility of macros within parentheses,
  macros are now evaluated only immediately before used in expressions,
  rather than the instant they are first encountered.

This redesign has allowed for many edge cases to be consolidated into
the general purpose parser, greatly simplifying the code.
  This commit includes a series of changes that get MOM6 to reproduce across
processor counts, layouts and restarts when USE_QG_LEITH_VISC = True and also to
implement related changes to the QG Leith code to be consistent with not making
the Boussinesq approximation when in non-Boussinesq mode.  These changes
include:

  - Adding a thermo_var_ptrs and a timestep argument to to the interfaces to
    horizontal_viscosity, initialize_dyn_split_RK2 and initialize_dyn_split_RK2b

  - Correcting a bug in the powers of H_subroundoff in the denominators of some
    expressions for the effective thicknesses associated with interfaces,
    thereby correcting some dimensionally inconsistencies in these expressions

  - Adding the new publicly visible routine calc_QG_slopes

  - Adding the arguments slope_x and slope_y to calc_QG_Leith_viscosity and
    avoiding the use of the VarMix%slope_[xy] elements in that routine

  - Using thickness_to_dz to consistently convert layer thicknesses to the
    vertical distances across layers and provide this as an argument to
    calc_QG_Leith_viscosity

  - Determining vertical distances in calculating the vertical derivatives of
    the slopes consistently with avoiding the Boussinesq approximation when in
    non-Boussinesq mode

  - Increasing some loop extents in calc_QG_Leith_viscosity

  Note that several of the expansions of the halo sizes in updates or in some of
the calculations are not necessary if an extra halo update is included for
slope_x and slope_y and dz after the calls to calc_QG_slopes and thickness_to_dz
in horizontal_viscosity.  This extra halo update has also been included in a
comment for later reference.

  This commit includes several changes to publicly visible interfaces, and it
will change answers with USE_QG_LEITH_VISC = True (which had previously
triggered a fatal error), but answers and output are bitwise identical in other
cases.
- This is required for our BOEM-funded Arctic project. It needs
  something similar for SIS2 as well.
For reasons not entirely clear, deallocation of the RK2 control
structure caused a segmentation fault when compiled by Nvidia.  A deeper
investigation suggested that the compiler was attempting to deallocate a
nullified pointer during automatic cleanup of RK2->SAL->sht.  Apparently
deallocation of a NULL pointer is an error (even though free() is not
supposed to care).

By redefining `sht` as an allocatable component rather than a local
component of the SAL_CS instance, it seemed to satisfy the compiler that
nothing was needed and the error went away.

There are still some lingering questions behind the cause of the
segfault, but for now I am going to put this under "compiler bug".

This patch also initializes some of the logical flags in the SAL_CS
type.  This is possibly unnecessary, but it is consistent with the
general rules of safety followed in other MOM6 derived types.
This patch fixes errors in the parser of perf output.  Previously,
each record was assumed to be separated by spaces, but this failed for
more generic records (usually from C++) which included signatures (such
as `f(a, b)`) or templates (`f<a, b>`).  Nested constructs were also
possible.

This is fixed by introducing a simple tokenizer which extracts <, (, and
whitespace from the output , then rebuilds the records by combining any
whitespace which appears inside of delimiters.

This patch should hopefully resolve the CI errors in GitHub Actions.
@marshallward
Copy link
Collaborator Author

AFAIK these regressions are expected. These are due to introduction of the bit-reproducible cube root function, and rewrite of certain diagnostics to make them rotationally symmetric, in preparation for FMA-reproducibility.

However, I really do not have a clue where these MacOS errors came from:

   internal compiler error: in gfc_walk_array_ref, at fortran/trans-array.cc:11870
  Please submit a full bug report, with preprocessed source (by using -freport-bug).
  See <https://github.com/Homebrew/homebrew-core/issues> for instructions.
  make[2]: *** [netcdf_io.o] Error 1
  make[1]: *** [fms/build/libFMS.a] Error 2
  make: *** [build/deps/lib/libFMS.a] Error 2

I haven't seen this before. I hope that it's temporary.

@Hallberg-NOAA
Copy link
Collaborator

Hallberg-NOAA commented May 16, 2024

Depending on how DEFAULT_ANSWER_DATE is set, in order to recover the previous solutions you might need to explicitly set ML_RESTRAT_ANSWER_DATE = 20231231 (or lower) (for NOAA-GFDL#545) and EPBL_ANSWER_DATE = 20231231 or lower (for NOAA-GFDL#541) to recover the previous solutions. The change to using reproducible answers in the diagnostics in NOAA-GFDL#605 can not be robustly recovered for any set of parentheses, and given that these diagnostics differ from before only at round-off, we opted not to provide a run-time flag to recover the previous answers. So, yes, depending on DEFAULT_ANSWER_DATE is set there are expected regressions. There are also expected regressions with fully non-Boussnesq configurations, but the nobody-is-using-it-yet exception that allows this will be ending very soon.

By contrast, the MacOS compiler time failures were not expected, but they appear to be somehow related to using a newer version of FMS. The rest of this MOM6 commit should still work with older versions of FMS that might not have this problem.

@marshallward
Copy link
Collaborator Author

This does not appear to be about MacOS but rather about GFortran 14.1. I reproduced this on my home machine:

mpifort -DPACKAGE_NAME=\"FMS\" -DPACKAGE_TARNAME=\"fms\" -DPACKAGE_VERSION=\"\" -DPACKAGE_STRING=\"FMS\" -DPACKAGE_BUGREPORT=\"https://github.com/NOAA-GFDL/FMS/issues\" -DPACKAGE_URL=\"\" -DHAVE_MPI=1 -DHAVE_GETTID=1 -DHAVE_SCHED_GETAFFINITY=1 -DHAVE_MPI=1 -Duse_libMPI=1 -DHAVE_STDIO_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_STDINT_H=1 -DHAVE_STRINGS_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_UNISTD_H=1 -DSTDC_HEADERS=1 -DHAVE_NETCDF_H=1 -Duse_netCDF=1 -g -O0 -fcray-pointer -fdefault-real-8 -fdefault-double-8 -fallow-invalid-boz -fallow-argument-mismatch  -c ../src/fms2_io/netcdf_io.F90 -I../src/fms2_io/include
register_variable_attribute.inc:56:60:

internal compiler error: in gfc_walk_array_ref, at fortran/trans-array.cc:11870
0x1ffd19a internal_error(char const*, ...)
	???:0
0x6f365f fancy_abort(char const*, int, char const*)
	???:0
0x844964 gfc_walk_expr(gfc_expr*)
	???:0
0x884767 gfc_conv_procedure_call(gfc_se*, gfc_symbol*, gfc_actual_arglist*, gfc_expr*, vec<tree_node*, va_gc, vl_embed>*)
	???:0
0x8833b3 gfc_conv_expr_reference(gfc_se*, gfc_expr*)
	???:0
0x886ed3 gfc_conv_procedure_call(gfc_se*, gfc_symbol*, gfc_actual_arglist*, gfc_expr*, vec<tree_node*, va_gc, vl_embed>*)
	???:0
0x8d252c gfc_trans_block_construct(gfc_code*)
	???:0
0x8c7980 gfc_trans_select_type(gfc_code*)
	???:0
0x8d252c gfc_trans_block_construct(gfc_code*)
	???:0
0x8c0caa gfc_trans_if(gfc_code*)
	???:0
0x87dcd7 gfc_generate_function_code(gfc_namespace*)
	???:0
0x839c89 gfc_generate_module_code(gfc_namespace*)
	???:0
0x7f1c4f gfc_parse_file()
	???:0
Please submit a full bug report, with preprocessed source (by using -freport-bug).
Please include the complete backtrace with any bug report.
See <https://gitlab.archlinux.org/archlinux/packaging/packages/gcc/-/issues> for instructions.
make[2]: *** [Makefile.dep:235: netcdf_io.o] Error 1
make[2]: Leaving directory '/home/marshall/gfdl/mom6/.testing/deps/fms/build'
make[1]: *** [Makefile:45: fms/build/libFMS.a] Error 2
make[1]: Leaving directory '/home/marshall/gfdl/mom6/.testing/deps'
make: *** [Makefile:342: deps/lib/libFMS.a] Error 2

I don't think this is a MOM6 issue. It just happens that our MacOS image is using a newer version than our Ubuntu image.

@marshallward
Copy link
Collaborator Author

Although I don't think this is really an FMS problem, I have opened the issue with them:

NOAA-GFDL/FMS#1527

I suppose we can either pray that GCC pushes an emergency fix (unlikely), or get FMS to release a modified version which works around the problem. This might require a new patch to MOM6 pointing to this new version.

... assuming that the new version works without any issues.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

10 participants