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

Merging dev/gfdl #11

Merged
merged 96 commits into from
Jul 25, 2023
Merged

Conversation

OlgaSergienko
Copy link
Owner

No description provided.

adcroft and others added 30 commits April 14, 2023 10:13
- Renamed function from psi(z) to mu(sigma)
- Added comments and units in function mu(sigma)
- Added [numerical] unit tests for mu(z), including special limits,
  special values, and one test value (checked against a python script).
Changes:

- Allow MLE parameterization to see surface buoyancy flux return from
  PBL scheme (affects MOM.F90, MOM_variables.F90:vertvisc_type,
  MOM_diabatic_driver.F90, MOM_set_viscosity.F90)

- Adds the Bodner et al., 2023, parameterization of restratification by
  mixed-layer eddies to MOM_mixed_layer_restrat.F90

  - This is a new subroutine rather than embedded inside the previous
    "OM4" version. It uses different inputs, different parameters,
    filters the BLD differently,

- Renamed mixedlayer_restrat_general to mxiedlayer_restrat_OM4 to better
  distinguish the two versions.

- Added function rmean2ts to extend the resetting running-mean time
  filter used in OM4 to use different time scales when growing or
  decaying. While mathematically the same in the limit of a zero
  "growing" time-scale, the implementation differs in the use of a
  reciprocal instead of division so was not added to the OM4 version.

- Updated module documentation

Co-authored-by: Abigail Bodner <abigail.bodner@gmail.com>
This patch adds the Bodner MLE testing parameters to the tc2.a test.
  Add the combined unit scaling factors Pa_to_RL2_T2 and Pa_to_RLZ_T2 to the
unit_scale_type to rescale pressures and wind stresses.  All answers are bitwise
identical, but there are two new elements in a public type.
  Use the new combined unit scaling factor US%Pa_to_RL2_T2 to rescale input
pressure fields and US%Pa_to_RLZ_T2 to rescale input wind stresses in various
places in the MOM6 code, including in the solo_driver and FMS_cap drivers.
Analogous changes could also be made to the mct and nuopc surface forcing files,
but have been omitted for now.  All answers are bitwise identical.
  Added the new runtime parameter TAUX_MAGNITUDE to set the strength of the
zonal wind stresses when WIND_CONFIG = "2gyre", "1gyre" or "Neverworld", with a
default that matches the previous hard-coded dimensional parameters that were
used to specify the wind stresses in these cases.  Also use US%Pa_to_RLZ_T2 to
rescale wind stresses throughout solo_driver/MOM_surface_forcing.F90.  By
default, all answers are bitwise identical, but there is a new runtime parameter
in the MOM_parameter_doc files for some test cases.
  Correct inconsistent dimensional rescaling of the input values of MLD_EN_VALS,
setting them all to [R Z3 T-2 ~> J m-2] to reflect that these are energies
associated with vertical turbulent mixing.  This fixes a rescaling bug when
these energies are set to non-default values at runtime, but all answers and
output are bitwise identical when no rescaling is used.
  Add better error handling to read_var_sizes when a missing file or missing
variable is provided as an argument.  Without this change the model fails with a
segmentation fault on line 768 of MOM_io.F90 if a bad file or variable name is
provided.  With this change, a useful error message is returned.  All answers
are bitwise identical in all cases that worked previously.
  Redid the scaling of 52 checksum or check_redundant calls for thickness or
transports to use the MKS counterparts of the thickness units (i.e., m and m3/s
or kg/m2 and kg/s, depending on the Boussinesq approximation), rather than
always rescaling them to m or m3/s.  In Boussinesq mode, everything remains the
same, but in non-Boussinesq mode, this means that the model's actual variable
are being checksummed and not a version that is rescaled by division by the
(meaningless?) Boussinesq reference density.  All solutions are bitwise
identical, but some debugging output will change in non-Boussinesq mode.
  Use a conversion factor to rescale the units of masscello, just like every
other diagnostic.  This does not change the diagnostic itself, but it changes
the order of the rescaling and the vertical remapping of this diagnostic onto
other coordinates (like z) or spatial averaging of this diagnostic, which can
change values in the last bits for this diagnostic for Boussinesq models (but
not for non-Boussinesq models, for which the conversion factor is an integer
power of 2).  As a result some of the diagnostics derived from masscello can
differ and this commit nominally fails the TC testing for reproducibility across
code versions.  All solutions and primary diagnostics, however, are bitwise
identical, and even the derived diagnostic calculations are mathematically
equivalent.
  Remove the code to account for unit rescaling within the restart files.  This
rescaling within the restart files has not been used in the code since March,
2022, and the model will work with older restart files provided that they did
not use dimensional rescaling, and even if they did they can be converted not to
use rescaling with a short run with the older code that created them.  Also
removed the publicly visible routines fix_restart_scaling and eliminated the
m_to_H_restart element of the verticalGrid_type; in any cases of non-standard
code using this element, it should be replaced with 1.0.  The various
US%..._restart elements and fix_restart_unit_scaling are being retained for now
because they are still being used in the SIS2 code.

  These changes significantly simplify the code, and they lead to a handful
of constants that are always 1 not being included in the MOM6 restart files.
All answers are bitwise identical, but a publicly visible interface has been
eliminated, as has been an element (GV%m_to_H_restart) of a transparent type.
  Added the new module MOM_EOS_Wright_full to enable the use of the version of
the Wright equation of state that has been fit over the larger range of
temperatures (-2 degC to 40 degC), salinities (0 psu to 40 psu) and pressures (0
dbar to 10000 dbar), than the does the restricted range fit in MOM_EOS_Wright,
which had been fit over the range of (-2 degC to 30 degC), (28 psu to 38 psu)
and (0 to 5000 dbar).  Comments have been added to both modules to clearly
document the range of properties over which they have been fitted.  The new
equation of state is enabled by setting EQN_OF_STATE = "WRIGHT_FULL".  In
addition, the default values for TFREEZE_FORM and EOS_QUADRATURE were changed
depending on the equation of state to avoid having defaults that lead to fatal
errors.  All answers are bitwise identical in any cases that currently work, but
there are new entries in the MOM_parameter_doc files.

  For now, only the coefficients have been changed between MOM_EOS_Wright and
MOM_EOS_Wright_full, but this means that it does not yet have all of the
parentheses that it should, as github.com/mom-ocean/issues/1331 discusses.
A follow up PR should add appropriate self-consistency and reference value
checks (with a tolerance) for the various EOS routines, and then add enough
parentheses to specify the order of arithmetic and hopefully enhance the
accuracy.  Ideally this can be done with the new equation of state before it
starts to be widely used, so that we can avoid needing a extra code to reproduce
the older answers.
  Cleaned up the comments describing the routines and added a proper doxygen
namespace block at the end of the MOM_EOS_Wright and MOM_EOS_Wright_full
modules, based on changes that A. Adcroft had on a detached branch of MOM6.
Only comments are changed, and all answers are bitwise identical.
  Added parentheses to all expressions with three or more additions or
multiplications in the MOM_EOS_Wright_full code, so that different compilers
and compiler settings will reproduce the same answers in more cases.  In doing
this, an effort was made to add the smallest terms first to reduce the impact
of roundoff.  In some cases, the code was deliberately rearranged to cancel
out the leading order terms more completely.  In addition, two bugs had been
identified in calculate_density_second_derivs_wright_full.  These were corrected
and the entire routine substantially refactored with renamed variables to make
the derivation easier to follow and verify.  Apart from the bug corrections in
the calculation of drho_dt_dt and drho_dt_dp, the changes in the expressions
are mathematically equivalent, but they might make the model less noisy in some
cases by reducing contributions from round-off errors.

  Also added comments highlighting two bugs in the drho_dt_dt and drho_dt_dp
calculations in calculate_density_second_derivs_wright in the original
MOM_EOS_Wright code, but did not correct them to preserve the previous answers.
  Created a new module, MOM_EOS_Wright_red, that uses the reduced range fit
coefficients from the Wright EOS paper, but uses the parentheses, expressions
and bug fixes that are now in MOM_EOS_Wright_full.  To use this new module, set
EQN_OF_STATE="WRIGHT_RED". This new form is mathematically equivalent using
EQN_OF_STATE="WRIGHT" (apart from correcting the bugs in the calculations of
drho_dt_dt and drho_dt_dp), but the order of arithmetic is different, so the
answers will differ.  This change is probably as close as we can come to
addressing the issues discussed at github.com/mom-ocean/issues/1331, so
that issue should be closed once this commit is merged onto the main branch.
Also corrected some misleading error messages in MOM_EOS and modified the code
to properly handle the case for equations of state (like NEMO and UNESCO) that
do not have a scalar form of calculate_density_derivs, but do have an array
form.  By default, all answers are bitwise identical.
  Corrected a sign error in calculate_spec_vol_array_linear and
calculate_spec_vol_scalar_linear when a reference specific volume is provided.
This bug will cause any configurations with EQN_OF_STATE="LINEAR" and
BOUSSINESQ=False (neither of which is the default value) to have the wrong sign
of the pressure gradients and other serious problems, like implausible sea
surface and internal interface heights.  This combination of parameters would
never be used in a realistic ocean model.  There are no impacted cases in any of
the MOM6-examples tests cases, nor those used in the ESMG or dev/NCAR test
suites, and it is very unlikely that any such case would work at all.  This bug
was present in the original version of the calculate_spec_vol_linear routines,
but it was only discovered after the implementation of the comprehensive
equation of state unit testing.  This will change answers in configurations that
could not have worked as viable ocean models, but answers are not impacted in
any known configuration, and all solutions in test cases are bitwise identical.
  Added the new publicly visible function EOS_unit_tests, along with a call to
it from inside of unit_tests.  These tests evaluate check values for density and
assess the consistency of expressions for variables that can be derived from
density with finite-difference estimates of the same variables.  These tests
reveal inconsistencies or omissions with several of the options for the equation
of state.  The EOS self-consistency tests that are failing are commented out for
now, so that this redacted unit test passes.  All answers are bitwise identical,
but there can be new diagnostic messages written out.
  Changed recently added doxygen labels in the two newly added EOS_Wright_red and
EOS_Wright_full modules to avoid reusing names that were already being used
by EOS_Wright.  All answers are bitwise identical, but the doxygen testing that
had been failing for the previous 5 commits is working again.
  Corrected numerous issues with the NEMO equation of state so that it is now
self consistent:

- Modified how coefficients are set in MOM_EOS_NEMO so that they are guaranteed
to be internally self-consistent, as verified by the EOS unit tests confirming
that the first derivatives of density with temperature and salinity are now
consistent with the equation of state.  Previously these had only been
consistent to about 7 decimal places, and hence the EOS unit tests were failing
for the NEMO equation of state.

 - Added new public interfaces to calculate_density_second_derivs_NEMO, which
had previously been missing.

 - Added code for calculate_compress_nemo that is explicitly derived from the
NEMO EOS.  The previous version of calculate_compress_nemo  had worked only
approximately via a call to the gsw package

  With these changes, the NEMO EOS routines are now passing the consistency
testing in the EOS unit tests.  Answers will change for configurations that use
the NEMO EOS to calculate any derivatives, and there are new public interfaces,
but it does not appear that the NEMO equation of state is in use yet, at least
it is not being used at EMC, FSU, GFDL, NASA GSFC, NCAR or in the ESMG
configurations.

This commit addresses the issue raised at github.com/mom-ocean/issues/405.
  Added the new public interface calculate_density_second_derivs_UNESCO, which
is an overload for both scalar and array versions, to calculate the second
derivatives of density with various combinations of temperature, salinity and
pressure.  Also added a doxygen block at the end of MOM_EOS_UNESCO.F90 to
describe this module and the papers it draws upon.  Also replaced fatal
errors in MOM_EOS with calls to these new routines.  All answers are bitwise
identical, but there are newly permitted combinations of options that previously
failed.
  Added the new public interface calc_density_second_derivs_wright_buggy to
reproduce the existing answers and corrected bugs in the calculation of the
second derivatives of density with temperature and with temperature and pressure
in in calculate_density_second_derivs_wright.  Also added the new runtime
parameter USE_WRIGHT_2ND_DERIV_BUG to indicate that the older (buggy) version of
calculate_density_second_derivs_wright is to be used.  Most configurations will
not be impacted, but by default answers will change with configurations that use
the Wright equation of state and one of the Stanley or similar nonlinear EOS
parameterizations, unless USE_WRIGHT_2ND_DERIV_BUG is explicitly set to True.

  This commit also activates the self-consistency unit testing with the Wright
equation of state (now that it passes) and limited unit testing of the TEOS-10
equation of state, omitting the second derivative calculations, one of which is
failing (the second derivative of density with salinity and pressure) due to a
bug in the TEOS10/gsw code.  Also added a unit test for consistency of the
density and specific volume when an offset reference value is used.
  Refactored the expressions in MOM_EOS_UNESCO.F90, adding parentheses to
specify the order of arithmetic, starting with the highest-order terms first for
less sensitivity to round-off.  Also added comments to better describe the
references for these algorithms.  Although the revised expressions are all
mathematically equivalent, this commit will change answers for any cases that
use EQN_OF_STATE = "UNESCO".  However, it is believed based on a survey of the
MOM6 community that there are no active configurations that use this equation of
state.
  Refactored the expressions in MOM_EOS_NEMO.F90, adding parentheses to specify
the order of arithmetic, starting with the highest-order terms first for less
sensitivity to round-off.  A number of internal variables were also renamed for
greater clarity, and a number of comments were revised to better describe the
references for these algorithms..  Although the revised expressions are all
mathematically equivalent, this commit will change answers for any cases that
use EQN_OF_STATE = "NEMO".  However, there is another recent commit to this file
that also changes answers (specifically the density derivatives) with this
equation of state, and it is believed based on a survey of the MOM6 community
that there are no active configurations that use this equation of state.
  Added the new equation of state module MOM_EOS_Roquet_SpV with the polynomial
specific volume fit equation of state from Roquet et al. (2015).  This equation
of state has also been added to MOM_EOS, where it is enabled by setting
EQN_OF_STATE="ROQUET_SPV".  Two other new valid settings have been added to
EQN_OF_STATE, "ROQUET_RHO" and "JACKETT_MCD", which synonymous with "NEMO" and
"UNESCO" respectively, but more accurately reflect the publications that
describe these fits to the equation of state.  The EoS unit tests are being
called for the new equation of state (it passes).  By default, all answers are
bitwise identical, but there are numerous new publicly visible interfaces.
  Added the new equation of state module MOM_EOS_Jackett06 with the rational
function equation of state from Jackett et al. (2006).  This uses potential
temperature and practical salinity as state variables, but with a fit to more
up-to-date observational data than Wright (1997) or UNESCO / Jackett and
McDougall (1995).  This equation of state has also been added to MOM_EOS, where
it is enabled by setting EQN_OF_STATE="JACKETT_06".  The EoS unit tests are
being called for the new equation of state (it passes).  This commit also adds
slightly more output from successful EoS unit tests when run with typical levels
of verbosity.  By default, all answers are bitwise identical, but there are
numerous new publicly visible interfaces.
  Added the routine calculate_specvol_derivs_UNESCO to calculate the derivatives
of specific volume with temperature and salinity to the MOM_EOS_UNESCO module.
Also added some missing parentheses elsewhere in this module so that the answers
will be invariant to complier version and optimization levels.  Also revised the
internal nomenclature of the parameters in this module to follow the conventions
of the other EOS modules.  Although the revised expressions are mathematically
equivalent, this commit will change answers for any cases that use EQN_OF_STATE
= "UNESCO".  However, it is believed based on a survey of the MOM6 community
that there are no active configurations that use this equation of state.  There
is a new publicly visible routine.
  Added the new publicly visible subroutine EOS_fit_range and equivalent
routines for each of the specific equation of state modules to return the range
of temperatures, salinities, and pressures over which the observed data have
been fitted.  This is also tested for in test_EOS_consistency to indicate
whether a test value is outside of the fit range, but the real purpose will be
to flag and then figure out how to deal with the case when the ocean model is
called with properties for which the equation of state is not valid.  Note that
as with all polynomial or other functional fits, extrapolating far outside of
the fit range is likely to lead to bad values, but things may not be so bad for
values that are only slightly outside of this range.  However the question of
how far out of the fit range these EoS expressions become inappropriate for each
of temperature, salinity and pressure is as yet unresolved.  All answers and
output are bitwise identical, but there are 10 new public interfaces.
  Removed unused and unnecessary #include <MOM_memory.h> statements from 5
equation of state modules.  All answers are bitwise identical.
  Refactored the specific volume calculations for the WRIGHT_FULL and WRIGHT_RED
equations of states for simplicity or to reduce the impacts of roundoff when
removing a reference value.  Also added code to multiply by the reciprocal of
the denominator rather than dividing in several places in the int_spec_vol_dp
routines for these same two equations of state, both for efficiency and greater
consistency across optimization levels.  These changes are mathematically
equivalent but will change answers at roundoff with these two equations of state,
but they are so new that they can not have been used yet.
  Renamed the module MOM_EOS_NEMO to MOM_EOS_Roquet_rho to more accurately
reflect its provenance, although setting either EQN_OF_STATE = NEMO or
EQN_OF_STATE = ROQUET_RHO will still work for using this code.  All answers
are bitwise identical, and previous input files will still work, but there are
some minor changes in the MOM_parameter_doc files.
Pperezhogin and others added 27 commits June 26, 2023 13:25
The icebergs project now includes drivers and tests which can interfere
with the coupled nolibs build, so we only pass its src directory to
mkmf.
  This commit includes changes to the get_param_real and log_param_real
interfaces to make the units arguments mandatory.  It also adds an optional
unscale argument to the log_param_real interfaces.  Without other changes in the
previous commits, it will cause the MOM6 code to fail to compile.  However, by
itself this commit does not change any answers or output.
- Update actions/checkout from v2 to v3 (suggested at
  NCAR#231 (comment) thanks
  to @jedwards4b)
The FMS2 function `get_unlimited_dimension_name` raises a netCDF error
if no unlimited dimension is found.  This is problematic for legacy or
externally created input files which may have not identifed their time
axis as unlimited.

This patch adds a new function, `find_unlimited_dimension_name` which
mirrors the FMS2 function but returns an empty string if none are found.

This is an internal function, not intended for use outside of the
module.
  Refactors the internal tide code in MOM_internal_tides and MOM_diabatic_driver
to consolidate it in the MOM_internal_tides module and allow the control
structure for that module to be made opaque.  This includes moving the internal
wave speed diagnostics and the call to wave_speeds or other code setting the
internal wave speeds into propagate_int_tide.  The get_param calls for
INTERNAL_WAVE_CG1_THRESH and UNIFORM_TEST_CG were moved from the diabatic module
to the MOM_internal_tides module. The wave_speed_CS and uniform_test_cg were
removed from diabatic_CS and added to int_tide_CS.  The Nb argument to
propagate_int_tide has been made intent inout, as it is now usually set via the
call to wave_speeds in that routine, but for certain tests it could use the
value passed in from diabatic_driver.

  All answers are bitwise identical, but there are changes to public interfaces
and types, and the order of some entries in the MOM_parameter_doc files and the
available_diags files is changed for some cases.
  Add new allocatable tau_mag arrays to the forcing and mech_forcing types to
hold the magnitude of the wind stresses including gustiness contributions.
There is also a new tau_mag diagnostic.  This same information in tau_mag is
being transformed into ustar, but these changes avoid division by the Boussinesq
reference density (GV%Rho0), and allow for a more accurate calculation of
derived fields when in non-Boussinesq mode, without having to multiply and
divide by GV%Rho0.  There is also a new optional tau_mag argument to
extract_IOB_stresses to support these changes.  These new arrays are not being
used yet in the MOM6 solutions, but they are being allocated and populated in
the routines that set the ustar fields, and they have been tested in changes to
the modules that use ustar that will come in a subsequent commit. This commit
also adds the new RLZ_T2_to_Pa element to the unit_scale_type to undo the
scaling of wind stresses and it makes use of it in some of the new code.  All
answers are bitwise identical, but there are new arrays or elements in three
transparent public types.
  Fixed several problems with the recently added Bodner mixed layer
restratification parameterization code.

- Corrected the dimensional rescaling in the expressions for psi_mag by adding a
  missing factor of US%L_to_Z.

- A logical branch was added based on the correct mask for land or OBC points to
  avoid potentially ill-defined calculations of the magnitude of the Bodner
  parameterization streamfunction, some which were leading to NaNs.

- Set a tiny but nonzero default value for MIN_WSTAR2 to avoid NaNs in some
  calculations of the streamfunction magnitude.

- Revised the expression for dd within the mu function in a mathematically
  equivalent way to avoid any possibility of taking a fractional exponential
  power of a tiny negative number due to truncation errors, which was leading to
  NaNs in some cases while developing and debugging the other changes that are
  not included in this commit.  This does not appear to change any answers in
  the existing test cases, perhaps because the mixed layer restratification
  "tail" is not being activated by setting TAIL_DH to be larger than 0.

- Corrected or added variable units in comments in the mixedlayer_restrat
  control structure.

  These could change answers (and avoid NaNs) in some cases with
USE_BODNER23=True, MLE_TAIL_DH > 0 or MLE%TAIL_DH > 0, and there will be changes
to the MOM_parameter_doc files for some cases, but given how recently this code
was added, it is expected that all answers are bitwise identical in the existing
test cases.
This patch adds wrappers to the set_domain and nullify_domain functions
used in FMS1 for internal FMS IO operations.  These are not used in
FMS2, so the wrapper functions are empty.

This is required to eliminate FMS1 IO dependencies in SIS2.
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 #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 <spencerjones@tamu.edu>
  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.
@codecov-commenter
Copy link

codecov-commenter commented Jul 25, 2023

Codecov Report

Merging #11 (af25f3f) into ice_dynamics (fc823f5) will increase coverage by 1.08%.
The diff coverage is 35.86%.

❗ Current head af25f3f differs from pull request most recent head 2f6e86e. Consider uploading reports for the commit 2f6e86e to get more accurate results

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

@@               Coverage Diff                @@
##           ice_dynamics      #11      +/-   ##
================================================
+ Coverage         37.08%   38.16%   +1.08%     
================================================
  Files               264      269       +5     
  Lines             74347    76588    +2241     
  Branches          13784    14085     +301     
================================================
+ Hits              27568    29228    +1660     
- Misses            41684    42106     +422     
- Partials           5095     5254     +159     
Files Changed Coverage Δ
...g_src/drivers/solo_driver/MESO_surface_forcing.F90 0.00% <0.00%> (ø)
...g_src/drivers/solo_driver/user_surface_forcing.F90 0.00% <0.00%> (ø)
config_src/infra/FMS1/MOM_coms_infra.F90 59.47% <0.00%> (-2.86%) ⬇️
config_src/infra/FMS1/MOM_io_infra.F90 22.00% <0.00%> (-1.72%) ⬇️
src/ALE/MOM_hybgen_regrid.F90 0.00% <0.00%> (ø)
src/ALE/coord_adapt.F90 0.00% <0.00%> (ø)
src/ALE/coord_hycom.F90 0.00% <ø> (ø)
src/ALE/coord_rho.F90 13.33% <ø> (+0.59%) ⬆️
src/ALE/regrid_interp.F90 1.47% <0.00%> (-0.03%) ⬇️
src/core/MOM_density_integrals.F90 16.92% <0.00%> (-0.19%) ⬇️
... and 90 more

... and 20 files with indirect coverage changes

@OlgaSergienko OlgaSergienko merged commit 9613191 into OlgaSergienko:ice_dynamics Jul 25, 2023
10 checks passed
OlgaSergienko pushed a commit that referenced this pull request Jul 8, 2024
Applying the first iteration explicitly appears to speed up the cuberoot
function by a bit over 20%:

 Before:
 Halley Final:  0.14174999999999999

 After:
 Halley Final:  0.11080000000000001

There is an assumption that compilers will precompute the constants like
`0.7 * (0.7)**3`, and that all will do so in the same manner.
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.

8 participants