*+Refactor, correct and augment the MOM6 EoS code #331
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Hallberg-NOAA
added
bug
|
This PR is substantial and changes answers. It should not be merged into dev/gfdl until after the pending voluminous PR from dev/gfdl has been created. |
Hallberg-NOAA
changed the title
Codecov Report
@@ Coverage Diff @@
## dev/gfdl #331 +/- ##
============================================
+ Coverage 37.11% 38.12% +1.01%
============================================
Files 264 269 +5
Lines 74204 75622 +1418
Branches 13750 13911 +161
============================================
+ Hits 27538 28831 +1293
- Misses 41582 41588 +6
- Partials 5084 5203 +119
... and 17 files with indirect coverage changes 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
|
I thought I'd try WRIGHT_FULL in the Arctic. The FMS2 build is no problem, the FMS1 build is giving me trouble with both gfortran 8.3.0 and 11.3.0. The first hint of trouble is: compared to: The Bering domain doesn't have this problem. Perhaps it's time to try Ifort again. |
|
Gah, bitten by FMS1's hatred of NaNs! My merged runoff file had no NaN, except as _FillValue. |
Hallberg-NOAA
force-pushed
the
EOS_parentheses
branch
from
cc29ca7 to
e5926bd
Compare
Hallberg-NOAA
force-pushed
the
EOS_parentheses
branch
from
e5926bd to
dcc4b93
Compare
|
@kshedstrom just checking, is this OK to merge? Not sure from your comment. |
|
Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/18829 ✔️ 🟡 |
|
Yes, it's fine, my fault. |
|
@Hallberg-NOAA I have a small cosmetic suggestion. It took me a bit time to realize that the "red" Wright EOS was in fact a "reduced" Wright EOS. If it were me, I'd just use REDUCED inplace of RED. But feel free to veto this suggestion. |
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.
Refactored MOM_EOS_Roquet_rho and MOM_EOS_Roquet_SpV to work directly with conservative temperatures in [degC] and pressures in [Pa] rather than normalizing them as in the original Roquet publication. However, the coefficients are still set using the values directly copied from that paper, but rescaled where they are declared as parameters, enabling (or requiring) compilers to precalculate them during compilation. These changes are mathematically equivalent but will change answers at roundoff with these two equations of state, but they are not believed to be in use yet.
Added the overloaded interface calculate_TFreeze_TEOS_poly to MOM_TFreeze to use the 23-term polynomial expression from TEOS-10 for the freezing point in conservative temperature as a function of pressure and absolute salinity. This gives results that agrees to within about 5e-4 degC with the algorithm used by calculate_TFreeze_TEOS10, which calls the gsw TEOS10 code that does an iterative inversion of a balance of chemical potentials to find the freezing point (see the TEOS10 documentation for more details). Also added testing for the freezing point calculations to the EOS_unit tests via the new internal subroutine test_TFr_consistency. This new freezing point calculation is invoked by setting TFREEZE_FORM = TEOS_POLY. By default, all answers are bitwise identical, but there are some minor changes in the comments in some MOM_parameter_doc files, and there are several new interfaces.
Added the new module MOM_temperature_convert, which contains the elemental functions poTemp_to_consTemp and consTemp_to_poTemp to convert potential temperature to conservative temperature and the reverse. These routines are mathematically equivalent to the TEOS-10 functions gsw_ct_from_pt and gsw_pt_from_ct, but with some refactoring and added parentheses to help ensure identical answers across compilers or levels of optimization. Also added the new subroutines pot_temp_to_cons_temp and prac_saln_to_abs_saln, and added the new optional argument use_TEOS to convert_temp_salt_for_TEOS10, and cons_temp_to_pot_temp and abs_saln_to_prac_saln. The equivalency between the new code and their gsw_ counterparts is demonstrated in new tests in the new function test_TS_conversion_consistency, which in turn is called from EOS_unit_tests. All answers are mathematically equivalent, but because of the choice to use the new code by default there could be changes at the level of roundoff in some cases that use conservative temperature as their state variable but initialize it from potential temperature. There are not any such cases yet in the MOM6-examples test suite, nor are there believed to be any such MOM6 configurations that are widely used. This commit introduces a new module and several new functions or subroutines with public interfaces.
Updated _Equation_of_State.dox to reflect the new options for the equation of state and freezing point calculations.
Eliminate use_TEOS optional arguments that were recently added to cons_temp_to_pot_temp and 4 other thermodynamic variable conversion functions, along with calls to gsw_pt_to_ct and similar conversion functions. All answers in the MOM6-examples test suite are bitwise identical.
Removed calculate_density_array from the overloaded public calculate_density interface, and similarly for the other EOS calculate_..._array routines, to help standardize how they are called. Calculate_density_derivs_array is the one exception is because it is being called from SIS2 and has to stay publicly visible for now. Additionally, the scalar and 1-d versions of the calculate_stanley_density routines were refactored to just use calculate_density and calculate_density_second_derivs call and avoid any EoS-specific logic, while the unused routine calculate_stanley_density_array is eliminated altogether. All answers are bitwise identical, including in extra tests that use the stanley_density routines.
Revised the setting EQN_OF_STATE to select the Wright equation of state with the reduced-range fit to "WRIGHT_REDUCED" (instead of "WRIGHT_RED") for greater clarity, in response to a comment in the review of the pull request with this sequence of code revisions. All answers are bitwise identical, but this changes the text for a recently added input parameter and it leads to changes in some comments in the MOM_parameter_doc files.
Hallberg-NOAA
force-pushed
the
EOS_parentheses
branch
from
dcc4b93 to
07e0602
Compare
|
Thanks for the suggestion to replace "WRIGHT_RED" with "WRIGHT_REDUCED". I have modified the parameter setting as suggested. |
marshallward
approved these changes
Apr 23, 2023
There was a problem hiding this comment.
Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/18879 ✔️ 🟡
Really amazing work by @Hallberg-NOAA to detect and sort out a wide spectrum of subtle issues in the equations of state.
This was referenced Mar 21, 2024
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR includes a series of commits that add 4 new equation of state (EoS) modules, add and deploy unit testing to the EoS modules, correct actual bugs in every one of the pre-existing modules, adds missing routines to the previous modules, and adds parentheses to most of the EoS modules (but not LINEAR or WRIGHT) to specify the order of arithmetic so that they will reproduce answers across compilers and levels of optimization. By default, answers are bitwise identical when EQN_OF_STATE = LINEAR or EQN_OF_STATE = WRIGHT, but a survey of the community suggested that the other equations of state were not being used yet, so answers will change due to both bug-fixes and changes in the order of arithmetic for all the other equations of state. Despite all of these changes, answers are bitwise identical for all of the cases in the MOM6-examples test suite and the TC tests, and the same is likely to be true of other currently used MOM6 configurations.
The new equation of state modules or settings for EQN_OF_STATE include:
The WRIGHT and TEOS10 equations of state give identical answers to what was there before.
The LINEAR equation of state for density is identical to what was there before, but a bug was fixed in the calculation of specific volume with a reference value was fixed; that but probably would have led to immediate failures of any case that used it.
The UNESCO and NEMO equations of state have been extensively refactored or had internal inconsistencies corrected and missing routines added, but the expressions for density are mathematically equivalent to what was there before.
The commits in this PR include: