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

Merge latest main (Feb 28, 2024) #271

Merged
merged 152 commits into from
Mar 6, 2024

Conversation

gustavo-marques
Copy link
Collaborator

No description provided.

marshallward and others added 30 commits July 7, 2023 10:48
The autoconf Python interpreter search was slightly modified to search
for Python even if $PYTHON is set to an empty string.  This is done by
unsetting PYTHON if it is set but empty, then following the usual macro.

This was required since `export PYTHON` in a Makefile will create the
`PYTHON` variable but will assign it no value (i.e. empty string).
This causes issues in some build environments.

The backup `configure~` script was also added to the developer
`ac-clean` cleanup rule.
  Refactored 6 files in the ALE directory to calculate the nominal depth in
thickness units in a single place (it is done in regridding_main now) and pass
it to the various places where it is used, in a preparatory step to modify how
this calculation is done in non-Boussinesq mode.  There are new arguments to
several publicly visible routines, including:

- Add non_depth_H arguments to hybgen_regrid, build_zstar_grid,
  build_sigma_grid, build_rho_grid, build_grid_HyCOM1, build_grid_adaptive,
  build_adapt_column and build_grid_arbitrary

- Add optional zScale arguments to build_zstar_grid and build_grid_HyCOM1

- Add unit_scale_type arguments to regridding_main, ALE_regrid_accelerated and
  ALE_offline_inputs

  Also eliminated an incorrect rescaling GV%Z_to_H facto when calculating the
total column thickness from the layer thicknesses when an ice shelf is used
with a Hycom grid.  This would have caused dimensional consistency testing to
fail.

  Added the new runtime parameters HYBGEN_H_THIN, HYBGEN_FAR_FROM_SURFACE
HYBGEN_FAR_FROM_BOTTOM, and HYBGEN_DENSITY_EPSILON to set previously hard-coded
dimensional parameters used in the Hybgen regridding code and store these
values in new variables in hybgen_regrid_CS.  Two of these are no longer passed
to hybgen_column_regrid as separate parameters.  By default these new runtime
parameters recover the previous hard-coded values.

  Also eliminated an unused block of code in build_rho_column.  Several
comments documenting variables or their units were also added.

  All answers are bitwise identical, but there are 4 new runtime parameters
that would appear in some MOM_parameter_doc files and there are changes to the
arguments to 11 routines.
As described in issue mom-ocean#372, I would like to be able to create restart files that
contain information about the particle location. These files will be written at
the same time as other restart files. I cannot add these calls directly to the
driver, because the driver does not have information about the particle location.
We have added save_MOM6_internal_state as a subroutine in MOM.F90, and
we added calls to this subroutine from each of the drivers. We hope this will
allow for more new packages to write restart files in the future.

Co-authored by Spencer Jones <[email protected]>
  Added the integer valid_SpV_halo to the thermo_var_ptrs type to indicate
whether the SpV_array has been updated and its valid halo size, to facilitate
error detection and debugging in non-Boussinesq mode.  Tv%valid_SpV_halo is set
to the halo size in calc_derived_thermo or after a halo update is done to
tv%SpV_avg, and it is set to a negative value right after calls that change
temperatures and salinities (such as by ALE remapping) unless there is a call to
calc_derived_thermo.  Tests for the validity of tv%SpV_avg are added to the
routines behind thickness_to_dz, with fatal errors issued if invalid arrays
would be used, but more tests could perhaps be used in any parameterization
routines where tv%SpV_avg is used directly.  Handling the updates to tv%SpV_avg
this way helps to avoid unnecessary calls to calc_derived_thermo, which in turn
has equation of state calls that can be expensive, while also providing
essential verification of new code related to the non-Boussinesq code.  These
tests can probably be commented out or removed for efficiency once there is a
full suite of regression tests for the fully non-Boussinesq mode of MOM6.  In
addition, a new optional debug argument was added to calc_derived_thermo which
can be used to triggers checksums for the variables used to calculate
tv%SpV_avg.  One call to calc_derived_thermo was also added just before the
initialization call to ALE_regrid that will be needed with the next commit, but
does not change answers yet.  All answers are bitwise identical, but there is a
new element in a transparent and widely used type and a new optional argument to
a public interface.
  Use RHO_KV_CONVERT instead of RHO_0 to set the non-Boussinesq version of
GV%m_to_H, so that there is a mechanism for testing the independence of the
fully non-Boussinesq mode from the Boussinesq reference density.  With this
change, GV%Z_to_H is not guaranteed to be equal to (GV%Z_to_m*GV%m_to_H), with
the latter expression preferred when setting parameters. By default the two
parameters are the same, and they will probably only ever differ in testing the
code.  All Boussinesq solutions are bitwise identical, but there are differences
in the description of RHO_KV_CONVERT that will appear in MOM_parameter_doc
files.
  Add new arguments to 7 routines that will be needed for the non-Boussinesq
capability, but do not use them yet, so that there will be fewer cross file
dependencies as the various changes are being reviewed simultaneously.  The
impacted interfaces are MEKE_int, vertvisc_coef, sumSWoverBands, KPP_calculate,
differential_diffuse_T_S, set_BBL_TKE, and apply_sponge In the three
step_MOM_dyn_... routines and in calculateBuoyancyFlux1d, this change includes
calls to thickness_to_dz to calculate the new vertical distance arrays that will
be passed into vertvisc_coef or sumSWoverBands.  The only place where the new
arguments are actually used is in sumSWoverBands and set_opacity where the
changes are particularly simple.  All answers are bitwise identical, but there
are new non-optional arguments to seven publicly visible routines.
* Restore functionality for reading slices from 3d volumes in MOM_io

 - The recent MOM_io modifications in support of FMS2_io accidentally
   removed support for reading on-grid data (same horizontal grid as
   model) k-slices. This is needed in some configurations in the model
   state initialization.

* Add FMS1 interfaces

* Additional patches to enable reading ongrid state initialization data

 - read local 3d volume rather than attempting to slice ongrid
   data vertically.

 - Related bugfixes in MOM_io
- We were reading KV_ML_INVZ2 without logging, then checking for KVML
  and finally logging based on a combination of the two. This had the side
  affect that we get warnings about not using KVML even if KVML was not
  present.
- The fix checks for KVML first, and then changes the default so that
  when KVML=1e-4 is replaced by KV_ML_INVZ2=1e-4 we end up with no
  warnings and KVML can be obsoleted safely.
  Note: this commit alone does not remove all warnings from the MOM6-examples
  suite because we still need to fix the MOM_input that still use KVML
- KVML needs to be unscaled since it is the default for KV_ML_INVZ2
- tc3 used KVML and has been corrected.
  Use the new runtime parameter RHO_PGF_REF instead of RHO_0 to set the
reference density that is subtracted off from the other densities when
calculating the finite volume pressure gradient forces.  Although the answers
are mathematically equivalent for any value of this parameter, a judicious
choice can reduce the impacts of roundoff errors by about 2 orders of
magnitude.  By default, RHO_PGF_REF is set to RHO_0, and all answers are bitwise
identical.  However, there is a new runtime parameter that appears in many of
the MOM_parameter_doc.all files.
The message that a file is being created was issued as a WARNING
when we all agree it should really be a NOTE. Depth_list.nc is
read if it is present to avoid recomputing a sorted list, but the
absence of the file is not an error and does not warrant a warning.

Changes:
- Changed WARNING to NOTE.
- Removed MOM_mesg from imports since it wasn't being used.
The interpolation scheme for state-dependent diagnostic coordinates
was incorrectly registering as the same parameter as the main model.
This meant it was never possible to change the interpolation scheme
from the default (which was not the same as the main model).

Fix registers the generated parameter name which was always computed
but not used. A typical example of the generated parameter is
"DIAG_COORD_INTERP_SCHEME_RHO2".
  Fixed a bug in which wave_speed_init was effectively discarding any values of
mono_N2_depth passed to it via the optional argument mono_N2_depth, but also
changed the default value of RESOLN_N2_FILTER_DEPTH, which was previously being
discarded, to disable the monotonization and replicate the previous results.
There were also clarifying additions made to the description how to disable
RESOLN_N2_FILTER_DEPTH.  This will change some entries in MOM_parameter_doc
files, and it will change solutions in cases that set RESOLN_N2_FILTER_DEPTH to
a non-default value and have parameter settings that use the resolution function
to scale their horizontal mixing.  There are, however, no known active
simulations where the answers are expected to change.
  Revised the calculation of gprime and the coordinate densities (GV%Rlay) in
fully non-Boussinesq mode to use the arithmetic mean of adjacent coordinate
densities in the denominator of the expression for g_prime in place of RHO_0.
Also use LIGHTEST_DENSITY in place of RHO_0 to specify the top-level coordinate
density in certain coordinate modes.  Also made corresponding changes to the
fully non-Boussinesq APE calculation when CALCULATE_APE is true, and eliminated
an incorrect calculation of the layer volumes in non-Boussinesq mode using the
Boussinesq reference density that was never actually being used when
CALCULATE_APE is false.

  This commit will change answers in some fully non-Boussinesq calculations, and
an existing runtime parameter is used and logged in some new cases, changing
the MOM_parameter_doc file in those cases.
  Refactored thickness_diffuse when in non-Boussinesq mode to avoid any
dependencies on the Boussinesq reference density, and to translate the volume
streamfunction into the mass streamfunction using an appropriately defined
in-situ density averaged to the interfaces at velocity points.  This form
follows the suggestions of Appendix A.3.2 of Griffies and Greatbatch (Ocean
Modelling, 2012) when in non-Boussinesq mode.  Thickness_diffuse_full was also
revised to work properly in non-Boussinesq mode (and not depend on the
Boussinesq reference density) when no equation of state is used.

  As a part of these changes, the code now uses thickness-based streamfunctions
and other thickness-based internal calculations in MOM_thickness_diffuse.  For
example, the overturning streamfunctions with this change are now in m3/s in
Boussinesq mode, but kg/s in non-Boussinesq mode.  These changes use a call to
thickness_to_dz to set up a separate variable with the vertical distance across
layers, and in non-Boussinesq mode they use tv%SpV_avg to estimate in situ
densities.  Additional debugging checksums were added to thickness_diffuse.

  The code changes are extensive with 15 new or renamed internal variables, and
changes to the units of 9 other internal variables and 3 arguments to the
private routine streamfn_solver.  After this change, GV%Rho, GV%Z_to_H and
GV%H_to_Z are no longer used in any non-Boussinesq calculations (12 such
instances having been elimated).  Because some calculations have to be redone
with the separate thickness and dz variables, this will be more expensive than
the original version.

  No public interfaces are changed, and all answers are bitwise identical in
Boussinesq or semiBoussinesq mode, but they will change in non-Boussinesq mode
when the isopycnal height diffusion parameterization is used.
* Salt data structures

* First steps at brine plume: pass info from SIS2

* The brine plume parameterization,

- including now passing the dimensional scaling tests.

* Fix problem when running Tidal_bay case with gnu.

* Avoiding visc_rem issues inside land mask.

Tweaking the brine plume code.

* Using the proper MLD in the brine plumes

- it now works better on restart

* Always including MLD in call to applyBoundary...

- I could move it up and make it not optional.

* Adding some OpenMP directives to brine plumes
This commit brings the drifters interface up-to-date with the current version of the drifters package, which requires h (layer thickness) to calculate the vertical movement of particles. The interfaces in the code and in config_src/external are updated to pass this information to the drifters package.
  Pass dt_kappa_smooth to calc_isoneutral_slopes and vert_fill_TS in units of
[H Z ~> m2 or kg m-1] instead of [Z2 ~> m2] for consistency with the units of
other diffusivities in the code and to reduce the depenency on the  Boussinesq
reference density in non-Boussinesq configurations.  In addition to the changes
to the units of these two arguments, there is a new unit_scale_type argument to
vert_fill_TS and MOM_calc_varT and a new verticalGrid_type argument to
MOM_stoch_eos_init.  The units of 4 vertical diffusivities in the control
structures in 4 different modules are also changed accordingly.

  All answers are bitwise identical in Boussinesq mode, but they can change for
some non-Boussinesq configurations.  There are new mandatory arguments to three
publicly visible routines.
  Added a comment justifying the use of a fixed rescaling factor for the
diffusivity used in vert_fill_TS.  All answers and output are identical.
  Added the new public interface find_ustar to extract the friction velocity
from either a forcing type argument, or a mech_forcing_type argument, either
directly or from tau_mag, and in non-Boussinesq mode by using the time-evolving
surface specific volume.  Find_ustar is an overloaded interface to
find_ustar_fluxes or find_ustar_mech_forcing, which are the same but for the
type of one of their arguments.  For now, the subroutines bulkmixedlayer,
mixedlayer_restrajt_OM4, mixedlayer_restrat_Bodner and mixedlayer_restrat_BML
are calling find_ustar to avoid code duplication during the transition to work
in fully non-Boussinesq mode, but it will eventually be used in about another
half dozen other places.

  All Boussinesq answers are bitwise identical, but non-Boussinesq answers will
change and become less dependent on the Boussinesq reference density, and there
is a new publicly visible interface wrapping two subroutines.
  Changed the units of the optional mono_N2_depth argument to wave_speed,
wave_speed_init and wave_speed_set_param in thickness units instead of height
units.  Accordingly, the units of one element each in the diagnostics_CS and
wave_speed_CS and a local variable in VarMix_init are also changed to thickness
units.  The unit descriptions of some comments describing diagnostics were also
amended to also describe the non-Boussinesq versions.  Because this is
essentially just changing when the unit conversion occurs, all answers are
bitwise identical, but there are changes to the units of an optional argument in
3 publicly visible routines.
  Added the new runtime parameter BT_RHO_LINEARIZED to specify the density that
is used to convert total water column thicknesses into mass in non-Boussinesq
mode with linearized options in the barotropic solver or when estimating the
stable barotropic timestep without access to the full baroclinic model state.
The default is set to RHO_0 and answers do not change by default.  This new
parameter is used in non-Boussinesq mode with some options in btcalc and
find_face_areas, when LINEARIZED_BT_CORIOLIS = True or BT_NONLIN_STRESS = False,
and in the unit conversion of the ice strength with dynamic pressure.

  Also cancelled out factors of GV%Z_to_H in MOM_barotropic.F90 to simplify the
code and reduce the dependence on the value of GV%Rho_0 in non-Boussinesq mode.
This involved changing the units of 4 variables in the barotropic_CS type, 3
internal variables in btstep and an internal variable in barotropic_init to use
thickness units.  The rescaled internal variable mass_to_Z was also replaced
with the equivalent GV%RZ_to_H.  There are also 4 new debugging messages.  Also
modified the units of the gtot_est argument to match those of pbce.  There is a
new element in barotropic_CS.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode.  Additionally there is a new runtime parameter that will
appear in some MOM_parameter_doc files.
  Use thickness_to_dz to convert layer thicknesses to depths in 4 tracer modules
(DOME_tracer, dye_example, ideal_age_example and nw2_tracers) so that this
conversion is done correctly in non-Boussineq mode, and there is no longer any
dependency on the Boussinesq reference density in that mode.

  This change includes the addition of a thermo_var_ptrs argument to 5 routines
(initialize_DOME_tracer, initialize_dye_tracer, dye_tracer_column_physics
ideal_age_tracer_column_physics and count_BL_layers) and changes to the units of
some internal variables, and the addition of 6 new 2-d or 3-d arrays with the
vertical distance across layers.  An unused param_file_type argument to
initialize_DOME_tracer was also eliminated.

  Comments were also added to describe the units of 5 of the variables in the
ideal age tracer control structure and 7 internal variables in that same module,
and there was some minor cleanup of the formatting cf calls in
tracer_flow_control_init.

  There was some minor refactoring in the ns2_tracers module to use SZK_(GV)
instead of SZK_(G) to declare the vertical extent of some arrays, and the
vertical indexing convention for interfaces in nw2_tracer_dist was revised from
starting at 0 to start at 1 for consistency with all the other code in MOM6.

  Also moved the code to do halo updates for the physical model state variables
and call calc_derived_thermo before calling tracer_flow_control_init, because
some routines there are now using the layer average specific volume to convert
between thicknesses and heights when in non-Boussinesq mode.

  All answers in Boussinesq mode are bitwise identical, but these passive tracer
modules have slightly different answers in non-Boussinesq mode.  There are
changes to the non-optional arguments to 4 public interfaces.
  Changed a recently added OMP directive for plume_flux from private to
firstprivate to reflect how this variable is actually used.  This bug was
introduced with PR mom-ocean#401, but was causing sporadic failures in some of our
pipeline tests with the intel compiler (essentially due to initialized
memory when openMP is used) for subsequent commits.
This patch merges the internal `save_restart` function with the new
`save_MOM6_internal_state` function into a new general MOM restart
function.

It also makes an effort to eliminate `MOM_restart` as a driver
dependency, narrowing the required MOM API for existing and future
drivers.

Also removes the `restart_CSp` argument from `MOM_wave_interface_init`,
since it appeared to be used for nothing.
MOM simulations typically abort of the restart directory (usually
RESTART) are absent.  This patch adds POSIX support for mkdir() and
creates the directory if it is missing.
Using inquire() to check for directory existence is not possible, since
at least one compiler (Intel) does not consider directories to be files.

The inquire call is replaced with a C interface to the POSIX stat()
function.  We do not fully emulate the behavior of stat, but we use its
return value to determine existence of directories.  This provides a
more reliable method for identifying the existence of the directory.

This should resolve many of the observed problems with RESTART creation
in coupled runs.
  Cancelled out factors of GV%Z_to_H in MOM_hor_visc.F90 to simplify the code
and reduce the dependence on the value of GV%Rho_0 in non-Boussinesq mode.  This
involved changing the units of 3 internal variables in horizontal_viscosity and
one element in the hor_visc_CS type to use thickness units or their inverse.
Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode.
  Revised the units of 12 vertvisc_type elements to be based on thicknesses, so
that vertical viscosities (in [H Z T-1 ~> m2 s-1 or Pa s]) are stored as dynamic
viscosites when in non-Boussinesq mode, with analogous changes to the diapycanl
diffusivity (now in [H Z T-1 ~> m2 s-1 or kg m-1 s-1]).  Similarly changed the
units of the 2 Rayleigh drag velocity elements (Ray_u and Ray_v) of the
vertvisc_type from vertical velocity units to thickness flux units and to more
accurately reflect the nature of these fields.  The bottom boundary layer TKE
source element (TKE_BBL) was also revised to [H Z2 T-3 ~> m3 s-3 or W m-2].
This commit also adds required changes to the units of the viscosities or
shear-driven diffusivities returned from KPP_calculate, calculate_CVMix_shear,
calculate_CVMix_conv, Calculate_kappa_shear, Calc_kappa_shear_vertex,
calculate_tidal_mixing and calculate_CVMix_tidal.

 Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
  Rescaled diapycnal diffusivities passed as arguments in non-Boussinesq mode,
to be equivalent to the thermal conductivity divided by the heat capacity,
analogously to the difference between a kinematic viscosity and a dynamic
viscosity, so that the new units are [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

  This includes changing the units of 4 arguments to set_diffusivity;
3 arguments each to calculate_bkgnd_mixing, add_drag_diffusivity,
add_LOTW_BBL_diffusivity, user_change_diff, calculate_tidal_mixing and
add_int_tide_diffusivity; 2 arguments to KPP_calculate, calculate_CVMix_conv,
compute_ddiff_coeffs, differential_diffuse_T_S, entrainment_diffusive,
double_diffusion, add_MLrad_diffusivity, and calculate_CVMix_tidal; and one
argument to energetic_PBL.

  The units of 36 internal variables were also changed, as were a total of
29 elements in various opaque types, including 8 elements in bkgnd_mixing_cs,
2 in diabatic_CC, 3 in tidal_mixing_diags type, 1 in user_change_diff_CS,
9 in set_diffusivity_CS type, and 6 elements in diffusivity_diags.

  Two new internal variables were added, and several redundant GV%H_to_Z
conversion factors were also cancelled out, some using that
GV%H_to_Z*GV%Rho0 = GV%H_to_RZ.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
  Changed the units for TKE arguments to [H Z2 T-3 ~> m3 s-3 or W m-2] for
find_TKE_to_Kd, add_drag_diffusivity, calculate_tidal_mixing and
add_int_tide_diffusivity, with similar changes to the units of 21 diagnostics
or internal variables in the same routines and in add_LOTW_BBL_diffusivity and
add_MLrad_diffusivity.  Dozens of unit conversion factors were also cancelled
out with these changes, including using that GV%Z_to_H/GV%Rho_0 = GV%RZ_to_H and
that GV%Rho_0*GV%H_to_Z = GV%H_to_RZ for both Boussinesq or non-Boussinesq
configurations.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
cspencerjones and others added 26 commits October 28, 2023 13:08
The directory, time and timestamp variables are needed by the particle code in order to write better restart files. I have changed the particles_save_restart interface to add these variables. I have also removed the option to pass temperature and salinity to particles_save_restart, because these variables are not useful for restart.
the module in now able to read in tidal velocities for different
tidal harmonics and distribute the energy and distribute TKE input over
the different vertical modes. This involves upsizing dimensions of
several arrays and mofiying some API. internal_tide_input_CS is
promoted to public to facilitate the passing of energy input to
MOM_internal_tides
  Removed the unused (and unusable) routine build_grid_arbitrary.  This routine
could not have been used because it had a hard-coded STOP call, and comments in
it indicated that it should have been deleted in July, 2013.  The run-time
parameter setting that would have triggered a call to this routine has been
retained for now, but with a fatal error message explaining that this routine
has not been implemented.   All answers are bitwise identical in any cases that
ran before.
  Removed the unused routine rescale_grid_bathymetry.  This routine was added in
August 2018 as a part of the development of the depth unit conversion and
dimensional consistency testing, but it is no longer being called now that this
conversion is essentially complete (and it has not been called by the code in
several years).  For the original commit that first added this code, see
github.com/mom-ocean/MOM6/commit/ddc9ed1c33a1b7357b213929118ecaa19ae63f9f.  All
answers are bitwise identical.
  Corrected dimensional rescaling bugs in the spherical harmonics SAL code.  An
issue with horizontal length scaling was corrected by using G%Rad_Earth_L in
place of G%Rad_Earth in spherical_harmonics_init.  There are new optional
tmp_scale arguments to calc_SAL and spherical_harmonics_forward to allow the
rescaling to be undone before calling the reproducing sums.

  This commit also modifies the call to the reproducing sums in
spherical_harmonics_forward so that all real or imaginary components are
calculated with a single call, which reduces the cost of the SAL calculation
reproducing sums from about 6.7 times the cost with non-reproducing sums to just
5.5 times as much in testing with the tides_025 test case.

  There is also code added to avoid NaNs arising from a square root operating
on a negative argument from a 32-bit integer roll-over when a very large
number of harmonics components (more than 1024 x 1024) are unadvisedly being
used.

  While this commit corrects the dimensional scaling when HARMONICS_SAL is true,
all answers are bitwise identical when no rescaling is used or when the
spherical harmonics SAL is not used.  There are new optional arguments to two
publicly visible interfaces.
  Added standard-format unit descriptions for 31 real variables in comments
scattered across 14 modules in the core, tracer, and both parameterizations
directories.  Only comments are changed and all answers are bitwise identical.
  Save tv%p_surf in the restart file when USE_PSURF_IN_EOS is true so that the
diagnosed potential energy written to the ocean.stats files after a restart
matches the energy written at the end of the previous run-segment in certain
non-Boussinesq configurations, including the Baltic test case.  Because
p_surf_EOS is a non-mandatory restart field, there is no problem restarting the
run from a restart file created by an older version of the model. The solutions
themselves are bitwise identical. This change requires that tv%p_surf is treated
as an allocatable pointer to its own array rather than being used as a pointer
to the p_surf element of the fluxes or forces structures so that it can be
registered as a restart field.  At some point tv%p_surf could be converted into
an allocatable array instead of a pointer, but this would require more extensive
code refactoring.  All answers are bitwise identical.
…ename specified by new parameter ICE_SHELF_ENERGYFILE). Currently, this file outputs the kinetic energy and mass of the ice sheet according to the smae parameters used to write the ocean.stats file (TIMEUNIT, ENERGYSAVEDAYS, and ENERGYSAVEDAYS_GEOMETRIC).
A few bug fixes so that the GL_couple=.true. option works correctly. Setting GL_couple=.true. will determine the grounding based on ocean column thickness rather than the typical the hydrostatic equilibrium condition. This has the advantage of accounting for changes in sea level, tides, etc. However, it has the disadvantage of not working with the same thoroughly-tested sub-element grounding line parameterization used for the hydrostatic condition. Instead, it accounts for sub-element grounding line movement by, during the SSA solution, using a grounding mask averaged over all ocean (sub)steps that completed since the last SSA solve. Unlike the hydrostatic sub-element parameterization, the dependence of the GL_couple=.true. scheme on grid resolution has not yet been determined. Qualitatively similar grounding line retreat/advance behavior is achieved with both approaches for MISOMIP IceOcean1 on a 2km grid, but GL_couple=.true. results in a rougher grounding line position with less retreat. Note that this commit also fixed a bug in applying the hydrostatic grounding line approach without its sub-element parameterization (though the sub-element parameterization should also be used anyway).
* This commit fixes a bug where restarts were not bitwise identical for coupled ice sheet/ocean runs. This required adding visc%taux and visc%taux (stress on the ocean under ice shelves) to the ocean restart file. Furthermore, ice shelf geometry-related variables needed to be updated and their halo cells filled (within subroutine update_ice_shelf) before (rather than after) calculating shelf fluxes and pressure to the ocean (subroutine add_shelf_flux). Also fixed a bug to calculate ISS%mass_shelf properly in ice-shelf cells that are only partially-filled (ice-shelf area < cell area) due to advection of the ice-shelf front.

* Modified the scheme that attempts to enforce a constant sea level in coupled ice-sheet/ocean runs, where balancing fluxes are applied to part/all of the open ocean to offset the sea level effects caused by changes in ice sheet mass on the ocean. The old scheme assumes the entire ice sheet is floating, and is retained here for the case where CS$override_shelf_movement==.true. and CS%mass_from_file (this approach is used for the MISOMIP tests without a dynamic ice sheet). The new scheme (which is needed for the MISOMIP tests with a dynamic ice sheet, e.g. IceOcean1r and IceOcean1a), accounts for the more general case where some of the ice sheet is grounded. In either case, the ocean balancing fluxes are applied over the entirety of the ice-sheet-free ocean by default. However, if the new parameter CONST_SEA_LEVEL_MISOMIP==.true., the balancing flux is only applied where x>=790 km, which is the sponge region for the MISOMIP tests.

The new scheme also requires calculation of ice sheet dHdT, which is now calculated and optionally saved as a useful diagnositc field for any ice sheet simulation.

* Eliminated array syntax and added inverse dt variables where needed
Makedep can now exclude prescribed directories in the directory tree
used to generate the file lists.

This is required for projects which may not follow normal development
processes, such as the FMS test programs.
The target (regression) configure step did not use a --with-framework
flag and would always build with FMS1, even if FRAMEWORK was set to
fms2.  This patch adds the flag to its configure step.

This patch also does some refactoring of the MOM_ENV and MOM_FCFLAGS
setup rules.  Values common to all rules are set externally, and
additional values for individual rules are appended.

Variable syntax also follows Makefile format (spaces around =) rather
than POSIX shell (no spaces).
The default FMS build in ac/deps is updated to 2023.03.

FMS source now includes a suite of test programs which require explicit
preprocessing macros, which can complicate out makedep-based build when
those macros are not present.  To avoid this, the new "skip" flag has
been added to the makedep build.

The skip flag should not cause errors or other issues in older versions
of FMS which do not have the excluded directory (though perhaps that
could or should change in the future).
  Revised the recently added brine-plume code to avoid using a vertical sum and
to use tv%SpV_avg to determine the brine plume mixing thickness instead of
GV%Z_to_H when in non-Boussinesq mode. Several internal brine plume expressions
now work in units of H instead of Z, and a total of 5 unit conversion factors
were eliminated. This will change certain non-Boussinesq calculations to avoid
any dependency on the Boussinesq reference density, but some Boussinesq answers
may also be changed and made more robust by avoiding the answer-indeterminate
sum() function.
* +Fix for issue mom-ocean#506, tracer OBC bug

 - it only happens in the advection for certain flow directions,
   inflow from OBC plus along-boundary flow.

* Tracer OBCs need more of an h halo update.

- This one should finally fix the processor count issues with OBCs.

* Correct the "if" statement.

* +Adding more halo points to an exchange

 - This will change answers if you start with a non-zero velocity.
   You need three halo points (or maybe cont_stencil) for the
   continuity solver.
 - Also trying to put in some initial DEBUG_OBC code.

* Fix some DEBUG_OBC logic

* Writing to temporary arrays for tres_xy

 - Way to trick some compiler.

* Another fix for DEBUG_OBC

* Fixing the whalo troubles for grids that are 2 wide/long.

* Exchange all the h_new points.

 - without this, because we're remapping all the tres points, it
 dies in debug mode on bad h_new values.

* Trying a different exchange

 - as it was, it was passing my tests, failing the auto-tests.

* Fixing the DEBUG_OBC logging

* Putting the logging statement back.

 - Making an error more verbose too.
- Added simple single-thread program to invoke EOS_unit_tests.F90
- Added not-as-simple program to time the cost of calculate_density()
  and calculate_spec_vol() for both scalar and array APIs
  - Placed in new directory config_src/drivers/timing_tests/
- Renamed MOM_unit_test_driver.F90 to test_MOM_file_parser.F90
- Updated .testing/Makefile
  - Added list of programs in config_src/drivers/unit_tests
    - These are added to BUILDS if DO_UNIT_TESTS is not blank.
      (DO_UNIT_TESTS was an existing macro but it might be uneeded)
    - These programs are compiled with code coverage
  - Added list of programs in config_src/drivers/timing_tests
    - These programs are compiled with optimization and no coverage
  - Fixed rule for building UNIT_EXECS (which did not re-build properly
    because the central Makefile was trying to model the dependencies
    even though those dependencies are in the build/unit/Makefile.dep)
  - Added convenient targets build.unit, run.unit, build.timing, run.timing
- Timing tests currently time a loop over 1000 calls (so that the resolution
  of the CPU timer is not too large) and 400 samples to collect statistics
  on timings. On gaea c5 this takes about 10 seconds.
  - The results are written to stdout in json.
- Added placeholder build and run of timing_tests to GH workflow.
  - Enabled for [push,pull_request]
  - We probably will not be able to use timings from GH but I still want
    to exercise the code we know the timing programs aren't broken by a
    commit.
- Also added driver for string_functions_unit_tests
  Refactored ALE_remap_velocities to separate the code setting the thicknesses
at velocity points from the code that actually does the remapping.  This
includes the creation of the new public routines ALE_remap_set_h_vel and
ALE_remap_set_h_vel_via_dz and the replacement the pair of tracer point
thickness arguments to ALE_remap_velocities and remap_dyn_split_RK2_aux_vars
with a pair of the old and new thicknesses at the velocity points and the
elimination of several arguments to these routines that are no longer being
used.  There are also new internal routines ALE_remap_set_h_vel_partial and
ALE_remap_set_h_vel_OBC to apply modifications to the velocity point thicknesses
with OBCs and one runtime option.  The runtime variable REMAP_UV_USING_OLD_ALG
has effectively been moved from MOM_ALE.F90 to MOM.F90, although it is still
being read in MOM_ALE_init for use with the accelerated regridding during
initialization.  All answers are bitwise identical, but there are two new
public interfaces and changes to the arguments to two other public interfaces
and a run-time parameter was moved between modules resulting in changes to
some MOM_parameter_doc files.
  This commit fixes a bug to ensure that OBC%debug is always being set when OBCs
are in use.  After a recent commit, the value of OBC%debug was not being set and
was being left in an indeterminate whenever DEBUG_OBC=False or DEBUG=False and
DEBUG_OBC was unspecified.  This would lead to the model writing out a number of
checksums in some cases with open boundary conditions enabled, depending on what
value was inherited by the uninitialized OBC%debug.   Although this does not
change any answers, it will avoid a problem that will write out a large volume
of undesired output and greatly slow down configurations with open boundary
conditions.
The most recent NCAR -> GFDL merge created an error (courtesy of myself)
which left the CFC concentration units in the post_data call, even
though these are now handled at registration.

This patch restores this expression and removes the unit conversion from
post_data.
This patch modifies select calculations of PR#1616 in order to preserve bit
reproducibility when FMA optimization is enabled.  We add parentheses and
reorder terms in selected expressions which either direct or suppress FMAs,
ensuring equivalence with the previous release.

We address two specific equations in the PR.  The first is associated
with vertical friction coupling coupling coefficient.  The diff is shown
below.

-  a_cpl(i,K) = Kv_tot(i,K) / (h_shear*GV%H_to_Z + I_amax*Kv_tot(i,K))
+  a_cpl(i,K) = Kv_tot(i,K) / (h_shear + I_amax*Kv_tot(i,K))

The denominator is of the form `a*b + c*d`.  A compiler may favor an FMA
of the form `a*b + (c*d)`.  However, the modified equation is of form
which favors the `a + c*d` FMA.  Each form gives different results in
the final bits.

We resolve this by expliciting wrapping the RHS in parentheses:

  a_cpl(i,K) = Kv_tot(i,K) / (h_shear + (I_amax*Kv_tot(i,K)))

Although this disables the FMA, it produces the same bit-equivalent
answer as the original expression.

----

The second equation for TKE due to kappa shear is shown below.

-  tke_src = dz_Int(K) *(((kappa(K) + kappa0)*S2(K) - kappa(k)*N2(K)) - &
-                              (TKE(k) - q0)*TKE_decay(k)) - &
+  tke_src = h_Int(K) * (dz_h_Int(K)*((kappa(K) + kappa0)*S2(K) - kappa(k)*N2(K)) - &
+                              (TKE(k) - q0)*TKE_decay(k)) - &
    ...

The outer equation was of the form `b + c` but is promoted to `a*b + c`,
transforming it to an FMA.

We resolve this by suppressing this FMA optimization:

  tke_src = h_Int(K) * ((dz_h_Int(K) * ((kappa(K) + kappa0)*S2(K) - kappa(k)*N2(K))) - &
                         (TKE(k) - q0)*TKE_decay(k)) - &
        ...

----

The following two changes are intended to be the smallest modification
which preserves answers for known testing on target compilers.  It does
not encompass all equation changes in this PR.  If needed, we could
extend these changes to similar modifications of PR#1616.

We do not expect to support bit reproducibility when FMAs are enabled.
But this is an ongoing conversation, and the rules around FMAs should be
expected to change as we learn more and agree on rules of
reproducibility.
@alperaltuntas alperaltuntas merged commit 13fbeb7 into NCAR:dev/ncar Mar 6, 2024
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.