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

Main to NCAR (30 July 2024) #291

Merged
merged 158 commits into from
Aug 9, 2024

Conversation

gustavo-marques
Copy link
Collaborator

Incorporating the changes from GFDL to main, as described in mom-ocean#1631.

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 mom-ocean#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 mom-ocean#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.
marshallward and others added 25 commits May 29, 2024 15:59
  Corrected an argument to the call to ALE_remap_scalar() when h_in_Z_units is
true in MOM_initialize_tracer_from_Z(), to avoid the problems documented in
github.com/NOAA-GFDL/issues/589.  A comment was also added explaining the
logic of what is going on in this fork of the code.  This commit will change
answers with some generic tracers that are initialized from a Z-space input
file, restoring them to previous values that worked previously (before about
Feb. 1, 2024 on dev/gfdl) in Boussinesq configurations without dimensional
consistency testing, but in a new form that does pass the dimensional
consistency testing for depths and thicknesses.  All answers are bitwise
identical in any cases that do not use generic tracers.
  Refactored the code in soliton_initialization.F90 to more accurately reflect
the nondimensionalization that was applied in developing this test case.  This
change includes reading in the maximum depth and beta, and using them to
calculate the equatorial deformation radius and the external gravity wave
speed.  The references to the papers describing the test case in the module were
added to the doxygen comments describing the routines in this module.   This
commit also includes adding comments documenting the nature and units of all the
internal variables in this module.  There are two new arguments each (param_file
and just_read) to soliton_initialize_thickness and soliton_initialize_velocity
to accommodate these changes, bringing them into line with the interfaces for
other similar user initialization routines, and MOM_initialize_state was changed
accordingly.  This change could change answers in general, but in the specific
examples that use this code, both beta and the external wave speed are
deliberately set to 1 in MKS units, so this commit does not change answers for
that specific case.
  Calculate the rescaled pseudo-salt flux when KPP is in use inside of
pseudo_salt_tracer_column_physics(), rather than using fluxes%KPP_salt_flux.
This commit also corrects for the fact that salinity and pseudo-salt have
scaling that differs by a factor of US%S_to_ppt, which was not previously being
taken into account.  All answers are bitwise identical when no dimensional
rescaling is being used, but they will change (and be corrected) when salinity
is being rescaled.
  Eliminated the KPP_salt_flux element of the forcing type, which is no longer
being used after the revisions to pseudo_salt_tracer_column_physics.  Also moved
5 allocatable arrays (KPP_NLTheat, KPP_NLTscalar, KPP_buoy_flux, KPP_temp_flux
and KPP_salt_flux) out of diabatic_CS, replacing them with ordinary arrays in
diabatic_ALE(), diabatic_ALE_legacy() and layered_diabatic().  This change could
reduce the high-water memory footprint of the model when KPP is in use.  All
answers are bitwise identical, but one element of a publicly visible type has
been eliminated.
  Corrected the descriptions of several conservative temperature and absolute
salinity arguments in the comments for arguments to 4 TEOS10 equation of state
routines, and eliminated the commented out code to mask out negative salinities
in 8 TEOS10 routines.  Only comments are changed, and all answers are bitwise
identical.
  Added the routine convert_MLD_to_ML_thickness to consistently convert the
mixed layer depths (in units of [Z ~> m]) back to the mixed layer thicknesses
(in [H ~> m or kg m-2]), with proper error checking and handling of the
non-Boussinesq case.  The body of this routine was taken from duplicated
blocks of code in MOM_mixedlayer_restrat.F90.  This is not tested directly
in this PR, but it was tested via revisions to MOM_mixedlayer_restrat.F90 that
are included in a subsequent commit.  All answers are bitwise identical, but
there is a new publicly visible interface.
  Pass mixed layer thickness, rather than mixed layer depth, to
ideal_age_tracer_column_physics, and determine the number of layers within
the mixed layer in count_BL_layers in thickness units rather than depth units.
To accommodate these changes there is now a call to convert_MLD_to_ML_thickness
inside of call_tracer_column_fns.  The thermo_var_ptrs type argument (tv) to
ideal_age_tracer_column_physics and count_BL_layers are no longer used and
have been removed.  All answers are bitwise identical in Boussinesq mode, but
in non-Boussinesq mode there are changes in the passive ideal age tracers at
the level of roundoff.  There are also changes in the units of arguments and
the number of arguments to a public interface.
  Provide the mixed layer thickness, as well as the mixed layer depth as
arguments to mixedlayer_restrat.  The code that had previously been used to
convert between the two has now been removed to the new external function
convert_MLD_to_ML_thickness.  To accommodate these changes, a new element,
h_ML, was added to the vertvisc type, and a call to convert_MLD_to_ML_thickness
was temporarily added in step_MOM_dynamics just before the call to
mixedlayer_restrat.  All answers are bitwise identical, but there
are changes to a public interface.
  Moved the call to convert mixed layer depths into mixed layer thicknesses
from right before the call to mixedlayer_restrat into the various diabatic
routines just after they are calculated.  This code rearrangement changes
answers in non-Boussinesq mode because the layer specific volumes will have
evolved between these two calls, but it is consistent with the mixed
layer thicknesses as used in the boundary layer parameterizations.

  A new argument was added to bulkmixedlayer to return the mixed layer
thickness that was already being calculated.  The previous argument Hml
was renamed to BLD and changed from a pointer into a simple array.  Both
are intent(inout) rather than intent(out) to preserve values in halos.

  This commit also revises the logic around the allocation of visc%h_ML and
its registration as a restart variable, including handling cases where an older
restart file is being read that includes MLD but not h_ML.  Because the mixed
layer depth is used in almost cases, CS%Hml in the control structure for
MOM.F90 was changed from a pointer to an allocatable, with values of 0 in those
cases (e.g., when ADIABATIC is true) where it is not used.

  In addition, the argument Hml to the various diabatic routines was changed to
MLD to more clearly reflect that it is a depth and not a thickness.  Outside of
diabatic, BML is only used to provide the tracer point boundary layer depths
to the calling routines, so it does not need a halo update.  Instead, all of
the halo updates for the various elements of visc that do need halo updates
after the diabatic calls are collected into one place for efficiency and more
accurate timings.

  The scale arguments of 34 checksum calls were revised from GV%H_to_m to
GV%H_to_MKS for more accurate checksums in non-Boussinesq configurations by
avoiding multiplication by an reference specific volume, and instead only
rescaling by an integer power of 2.

  All answers are bitwise identical in Boussinesq mode, but in non-Boussinesq
mode there are changes in answers in cases that use the boundary layer
thicknesses obtained from the boundary layer parameterizations in subsequent
calculations, such as mixed layer restratification with some options.  There are
also changes to the arguments of several publicly visible routines.
  Pass visc%h_ML to applyBoundaryFluxesInOut, which avoids the need for the
thickness conversions inside of that routine.  Answers are bitwise identical for
Boussinesq cases, but they will change in non-Boussinesq cases because the layer
specific volumes have changed between the time the boundary layer thicknesses
have been calculated and when they are being used in the brine plume
parameterization.

  Also use visc%h_ML in place of calls to KPP_get_BLD or energetic_PBL_get_MLD
in step_MOM_dyn_split_RK2, hor_bnd_diffusion and neutral_diffusion_calc_coeffs.
This change requires a new vertvisc_type argument to tracer_hordiff,
hor_bnd_diffusion and neutral_diffusion_calc_coeffs.  Boussinesq answers are
bitwise identical, but non-Boussinesq answers change because they now use the
actual model specific volumes rather than some prescribed constant value to to
the relevant unit conversions.

  The logic in set_visc_register_restarts was also updated so that visc%MLD and
visc%h_ML are allocated and registered for restarts in the cases where they will
be used.

  All answers are bitwise identical in Boussinesq mode, but in non-Boussinesq
mode there are changes in answers in cases that use neutral tracer diffusion
that excludes the boundary layer, horizontal tracer diffusion in the boundary
layer or FPMIX.  There are also changes to the arguments of several publicly
visible routines.
  Added the optional argument h_BL to call_tracer_column_fns, and then pass
visc%h_ML to this routine as a part of the calls from the diabatic routines.
This argument is not provided in adiabatic mode or in offline tracer mode, which
is why it has been made optional.  All answers are bitwise identical in
Boussinesq mode but ideal age tracers can change in non-Boussinesq mode due to
the updates to the layer specific volumes between the calculation of the
boundary layer properties and the calls to call_tracer_column_fns.
diag_send_complete() is removed from disable_averaging(), due to
potential performance costs over sweeping through files and variables
for synchronization across nodes or other thread-like backround
operations.

It also prevents potential errors in drivers where
diag_manager_set_time_end was unset.

We will come back to this issue and work out a better strategy for
diagnostic synchronization.
Three FMAs were introduced when the BBL convex L(:) computation were
moved to a separate function.  Parentheses have been added to disable
these FMAs.

Calculation of the cell width fractions `L` was consolidated and moved
into a new function, `find_L_open_concave_analytic` (later renamed to
`find_open_L_open_concave_trigonometric`).  When compiled with Intel
and using the following flags,

    -O2 -march=core-avx2 -fp-model source -qno-opt-dynamic-align

additional FMAs are added to the function, specifically the following
expressions:

    crv_3 = (Dp + Dm - 2.0*D_vel)

    L(K) = L0*(1.0 + ax2_3apb*L0)

    L(K) = apb_4a * (1.0 - 2.0 * cos(
        C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)
    )

(Third expression FMA was added inside of `acos()` by the compiler)

The crucial flag seems to be `-march=core-avx2` which controls the
introduction of FMAs and explains why the problem went undetected.

The legacy answers were restored by adding new parentheses to the above
equations.

Since this function is not generally considered reproducible and is
provided for legacy support, the loss of performance is not considered
to be an issue.
Five fields associated with KPP were previously initialized during
allocation (using `source=0.`).  It was found that the allocations were
unnecessary, and the fields were redefined as local variables within
functions, without the initialization to zero.

Although the uninitialized points were masked and unused in the rest of
the model, they were still able to change the checksums and were being
reported as regressions.  This patch restores the original checksums by
re-implementing the zero initialization.

Thanks to Bob Hallberg (@Hallberg-NOAA) for diagnosing this issue.
Due to answer changes in multiple production experiments, this patch
introduces a new parameter, USE_HUYNH_STENCIL_BUG, which sets the tracer
advection stencil width to its previous incorrect value.  This parameter
replicates the previous method, which keeps `stencil` at 2 when PPM:H3
is used.

----

The logical parameters are

* P = CS%usePPM
* H = CS%useHuynh
* E = USE_HUYNH_STENCIL_BUG
* R = CS%useHuynhStencilBug

The three expressions of interest:

1. P & ~H
2. P
3. P & ~R

(1) is the original incorrect expression, (2) is the fixed expression,
and (3) is the proposed update.

What we want:

If E is false, then (3) should reduce to (2).  If E is true, then (3)
should reduce to (1).

R is computed as follows:

* R = False  (from derived type initialization)
* if (H) R = E  (from get_param call)

This is equivalent to R = H & E, and (3) is P & ~(H & E).

* If E is False, then P & ~(H & False) = P.
* If E is True, then P & ~(H & True) = P & ~H

So this flag should replicate both the previous and current behavior.
  This commit restores the effectiveness of the runtime parameter
USE_WRIGHT_2ND_DERIV_BUG in determining whether a bug is corrected in the
calculation of two of the second derivative terms returned by
calculate_density_second_derivs_elem() with the "WRIGHT" equation of state,
recreating the behavior (and answers) that are currently on the main branch of
MOM6.  To do this, it adds and calls the new routine set_params_buggy_Wright()
when appropriate, and adds the new element "three" to the buggy_Wright_EOS type.
When the bug is fixed, buggy_Wright_EOS%three = 3, but ...%three = 2 to
recreate the bug.  This commit does change answers for cases using the "WRIGHT"
equation of state and one of the "USE_STANLEY_..." parameterizations from those
on the dev/gfdl branch of MOM6, but in so doing it restores the answers on the
main branch of MOM6.  There is also a new publicly visible subroutine.
  Moved a check for whether to initialize an uninitialized visc%h_ML field from
visc%MLD outside of the test for the value of CS%mixedlayer_restart.  This
change will prevent answer changes between code versions for certain cases that
initialize from a restart file, use the melt_potential field and have
MLE_USE_PBL_MLD = True but also have MIXEDLAYER_RESTRAT = False.  This commit
corrects a very specific oversight that was introduced when the roles of
visc%h_ML and visc%MLD were separated, and it can change answers (reverting them
to answers from older versions of the main branch of MOM6) in very specific
cases where MOM6 is initialized from a restart file from an older versions of
the code, including the 5x5-degree test case in the dev/emc regression suite.
  Added checksum calls for the melt_potential, ocean_mass, ocean_heat and
ocean_salt elements of the surface state in MOM_surface_chksum if these fields
are allocated for more comprehensive debugging.  Also added the symmetric
optional argument to the call to MOM_surface_chksum form extract_surface_state
so that all of the surface velocity values that could contribute to the ocean
surface velocities that are seen by the coupler are checksummed.  All solutions
are bitwise identical, but there are enhancements to the MOM6 debugging
capabilities.
The argument missing_scale was missing in the calls to read_Z_edges.
This patch adds missing_scale=1.0 to these calls. In the future,
we might want to consider adding an option to pass a factor to
scale the output tracers from the units in the input file.
@gustavo-marques
Copy link
Collaborator Author

@mnlevy1981; please check this commit 7afbb6d

I add to add the argument missing_scale was missing in the calls to read_Z_edges. This patch adds missing_scale=1.0 to these calls. In the future, we might consider adding an option to pass a factor to scale the output tracers from the units in the input file.

@gustavo-marques gustavo-marques marked this pull request as ready for review August 9, 2024 19:57
Copy link
Member

@alperaltuntas alperaltuntas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tested again, this time using the pr_mom suite, and encountered no answer changes.

@alperaltuntas alperaltuntas merged commit e96a46e into NCAR:dev/ncar Aug 9, 2024
8 of 10 checks passed
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.