GFDL to main (2026-07-27) - #1708
Conversation
Fix the 3-equation iteration for the buoyancy flux between the ocean and an overlying ice-shelf when ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX is true and SHELF_3EQ_GAMMA it false. This code now uses proper bounding of the self-consistent solution, avoiding further amplifying the fluxes in the cases when the differences between the diffusivities of heat and salt to make the buoyancy flux destabilizing for finite turbulent mixing. Both the false-position iterations and the (appropriately chosen) Newton's method iterations have been extensively examined and determined to be working correctly via print statements that have subsequently been removed for efficiency. Previously, the code to determine the 3-equation solution for the buoyancy flux between the ocean and an ice shelf had been skipping iteration altogether or doing un-bounded Newton's method iterations with a sign error in part of the derivative, including taking the square root of negative numbers, leading to the issue described at #945. That issue has now been corrected and can be closed once this commit has been merged into the dev/gfdl branch of MOM6. This commit also changes the names of the runtime parameters to correct the ice shelf flux iteration bugs from ICE_SHELF_BUOYANCY_FLUX_ITT_BUG and ICE_SHELF_SALT_FLUX_ITT_BUG to ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX and ICE_SHELF_SALT_FLUX_ITT_BUGFIX to avoid confusion with other ..._BUG parameters where `true` is to turn the bugs on, whereas here `true` fixes them. The old names are retained via `old_name` arguments to the `get_param()` calls, so no existing configurations will be disrupted by these changes. Additionally, an expression to determine a scaling factor to limit ice-shelf bottom slopes in `calc_shelf_driving_stress()` was refactored to avoid the possibility of division by zero. This commit will change (and correct) answers for cases with ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX set to true, but as these would often fail with a NaN from taking the square root of a negative value, it is very unlikely that any such configurations are actively being used, and there seems little point in retaining the previous answers. No answers are changed in cases that do not use an active ice shelf. Co-authored-by: Alistair Adcroft <adcroft@users.noreply.github.com>
- rotated OBC%segment%num_fields needs to be set.
Fix a bug that the recently changed default answer date for TIDES_ANSWER_DATE is not properly applied to MOM_self_attr_load. TIDES_ANSWER_DATE is used in MOM_self_attr_load to check if SAL_USE_BPA is used after a timestamp, so its default should be consistent with MOM_PressureForce_FV.
Thanks to both @alperaltuntas and @marshallward who noted that a PPM limiter has the expression `( u2 - u1 ) * ( u1 - u0 ) <- 0.0` which is interpreted as `( u2 - u1 ) * ( u1 - u0 ) < -0.0a. Needless to say, the intended code was `( u2 - u1 ) * ( u1 - u0 ) <= 0.0`. The same typo was copied to three files. The high-order estimate of edge value was previously bounded by (u2,u1) or (u1,u0). The missed conditions of either `( u2 - u1) == 0.` or `( u1 - u0 ) == 0.` would then have been caught by the subsequence test for an interior extrema. Thus, I think the cell was still limited to PCM appropriately. However, the typo obscured the intention of the limiter and I was lucky it still worked.
This commit allows the frequency-dependent drag to be implemented in tensor form, by incorporating the off-diagonal components of the wave drag tensor into the MOM_wave_drag module.
Recon1d_PLM_WLS provides a piecewise linear reconstruction where the slope is the "best" fit as determined by volume-weighted least squares. The reconstruction is NOT limited by neighboring cells. Therefore, this reconstruction is NOT useful for vertical remapping or grid generation. It is instead intended for the pressure gradient calculation; the idea is to disconnect the PLM slope from the values in vanish(ing) layers which appear to be the source of pressure-gradient errors over topographic slopes in z*-coordinate tests. Because the normal limiters do not apply, the only test I could think of was to check that the least squares fit was actually correct. The documentation explains how this was checked (which took a while due to round-off challenges with the loss function).
Corrected the descriptions of variable units in 64 comments spread across 16 files, including a dozen instances where "arbitrary" was misspelled. All answers are bitwise identical and only comments were changed.
Updated the values of about 21 parameters (many of which are repeated across TC test cases) used in the TC testing to test the most recent versions of code that is selected with ANSWER_DATE flags and to avoid testing the buggy versions of code that is regulated by _BUG flags. This includes some changes to broaden the range of equations of state that are being tested and to test some newer versions. This does change the details of the TC tests, but they should (and do) still pass TC regression tests across code versions.
* Added frazil to ice shelf The frazil mass flux to the ice-shelf base is calculated by multiplying frazil energy [J m-2] by the inverse of the timestep times the latent heat of fusion [kg J-1 s-1]. This frazil mass flux is incorporated as a negative water flux from the ice shelf. This negative water flux then acts to add the frazil mass to the ice shelf base (MOM_ice_shelf.F90/change_thickness_using_melt) and remove it from the ocean surface as evaporation (MOM_ice_shelf.F90/add_shelf_flux). Note frazil is reset to zero at the start of each therm timestep in MOM.F90/step_MOM. Some additional changes were also made to how the ice-shelf flux factor is implemented, so that is only scales ice-shelf melt without affecting the frazil mass flux. * Fixed a commented line where fluxes%water_flux should be ISS%water_flux
The PLM reconstruction used within the pressure gradient force now supports the weighted least squares approach for slope estimation. In a catastrophic version of seamount/z where vanished layers slightly inflate, the regular finite volume PLM method is sensitive to values in the vanished layers and leads to a feedback that causes en error growth (spontaneous motion). The PLM-WLS method is insensitive to the vanished layers and in the same test leads only to round-off level noise in the flow.
* Spatially varying bottom drag coefficient The spatially varying bottom drag coefficient can be specified by providing a map of the spatially varying scaling factor. * Spatially varying bottom drag coefficient Fixed the inconsistency at open boundaries when CDRAG_MAP is true.
In a number of cases, total resting column thickness is calucated as
G%bathyT + G%Z_ref, which is largely correct but for wetting, i.e.
G%bathyT < 0. This commit makes a correction for seven cases with this
potential bug.
There is no answer changes if no wetting points are used and G%Z_ref is
zero.
List of modules/processes affected:
* MOM_barotropic
* affects only surface stress when BT_NONLIN_STRESS is False.
* MOM_wave_speed
* h2 calculations in
* subroutine internal_tides_init
* subroutine int_tide_input_int
* subroutine tidal_mixing_init
* MOM_lateral_mixing_coeffs
* MOM_MEKE
In commit b8c807b, we made the test for SSH penetrating the sea floor when using BT_LIMIT_INTEGRAL_TRANSPORT because we thought it could never happen. Unfortunately, floating-point round off allows violations and we were hitting the now fatal error. This commit calculates the precision we can expect for the current SSH and then if the ocean thickness has become negative within this precision, we reset to zero thickness. This should not change answers in that BT_LIMIT_INTEGRAL_TRANSPORT is a new option, and if anyone was using it they would have encountered a FATAL, and this fix does not alter any positive thicknesses.
When debugging with all run-time tests turned on, the integer `num_lines` was flagged as used but uninitialized when being passed to `broadcast()`. I don't think the code was wrong, just that the checks expected the "inout" argument to be set on all processors when the purpose of `broadcast()` is to take the value from root_PE and send to everyone else. I don't know why this hadn't been detected before - maybe compiler version. The fix is trivial and has no impact on production codes or answers.
`oldfn` was not initialized when used in a logical test. This did not matter for numerical results; the logical expression always evaluated to the False correctly due to other parts of the expression. Nevertheless, this variable was technically used uninitialized and a debugging executable doesn't get past this. Hence the fix.
When debugging the ice sheet configuration, a non-zero barotropic transport could not be reconciled with the layer transports because the derivative of net layer transports was zero (d/dv hu). This arose due to all layer flows pointed from vanished to thick so that their marginal thicknesses were individually zero. Adding a floor to the marginal thickness allows the solver to find the adjustment that does reconciles the two estimates. I've made this optional via parameter CONT_USE_H_MARG_MIN, and with default of False. If this situation had occurred before, we surely would have had a crash so it's likely that always applying this floor would not change answers. However, there's the weak possibility that a teeny-tiny transport, smaller than H_subroundoff, has existed in a run and then this answer would change. With the default of False we can be sure there are no answer changes, but it is recommended to use this option for safety.
The CHANNEL_DRAG option was using a harmonic mean to interpolate adjacent bottom depths at velocity points to vorticity points. However, this is not well behaved when the bottom depth is negative (i.e., above sea level), as was noted as a part of PR #975. This commit adds the new runtime parameter CHANNEL_DRAG_SHELFBREAK_DEPTH to set a depth below which a harmonic mean bottom depth is still used to mimic a continental shelfbreak profile, but above which a simple arithmetic mean is used to interpolate bathymetry to vorticity points for use with CHANNEL_DRAG. The expressions vary continuously with depth and avoid the previous problems with division by zero or a badly formed harmonic mean. By default, all answers are bitwise identical in any cases that worked previously, but cases with oceans (or Great Lakes) in basins with bottoms that are above sea-level should now work sensibly when CHANNEL_DRAG is enabled. There is a new runtime parameter in some cases.
Corrected the incorrect or inconsistent unit descriptions of 28 variables, added descriptions of the units of 4 others, and corrected the non-standard syntax (e.g. backwards or in the wrong order) in the description of 35 variables, scattered across 27 files. Only comments are changed and all answers are bitwise identical.
* Fix for ice-shelf friction velocity bugs Fixed an incorrect area used to calculate cell-centered ocean surface velocity under the ice_shelf, which can impact the calculation of ice-shelf friction velocity. Added missing flags to some allocate_surface_state calls so that sfc_state%taux_shelf and sfc_state%tauy_shelf are allocated. This is required for the surface-stress-based (rather than surface-velocity-based) calculation of ice-shelf friction velocity. Also added taux_shelf and tauy_shelf as diagnostics for the surface stress under the ice shelf. * Removed unneeded taux_shelf and tauy_shelf diagnostics * Added ustar_from_vel_bugfix flag, which if true, fixes the ustar from ocean velocity bug
This patch undoes a coupling of the FMS infra layer to the MOM6 framework code. In the current FMS infra layers, the `get_extern_field_info()` and `init_extern_field()` functions require content defined in `src/framework`. This prevents the development of new independent infra layers, which much also depend on infra-agnostic content. In particular, the FMS2 implementation of `get_extern_field_axes()` relies exclusively on the framework function, `get_var_axes_info()`. Both infras also return the `axes_info` type, a MOM-specific framework-level descriptor, rather than the infra `axistype`. This patch resolves these inconsistencies. * `axis_info` no longer appears at infra-level. All relevant functions now reference `axistype`. * `src/framework/MOM_io.F90` now provide functions for translating `axistype` to `axis_info`. Some specific changes are summarized below. * `get_external_field_info` is now a framework-level function of `MOM_interpolate.F90` , using infra-level implementations of `get_extern_field_(size|axes|missing)`. Each is now explicitly defined at the infra-level. * The FMS2 `get_external_field_axes` is now an entirely new function, and is largely a duplicate of `get_var_axes_info()`. The major difference is that it returns a list of `axistype`. It also replaces the fixed x-y-z fetch with a slightly more generic list of axes. (It still requires at least three dimensions, however.) * `set_axis_data` is only used internally by the FMS2 infra. It is included in FMS1 but raises an nonimplementation error. There is one minor API change. * The `name` argument was added to `get_axis_data`. It is now the second argument, to match the style of existing functions, and size was moved to the third argument. Other minor framework references have been removed. * `MOM_error` and `FATAL` now refernce their `MOM_error_infra` equivalents. * `lowercase`, which was previously only defined in FMS1, has been added to the FMS2 infra. Note that this is a duplication of the function in `src/framework/MOM_string_functions.F90`.
When a z-coordinate diagnostic grid is specified via the "PARAM" method of coordinate definition, then the number of levels was always the same as the main model. This commit fixes this by first allowing for upto a 1000 levels in the new grid, checking for the actual requested size, and then allocating to that size. It appears we have no examples using this mode, which is probably how this bug has persisted so long. This "PARAM" method of specifying grids is being used in a range of new CMIP7 diagnostics in both MOM6 and COBALT.
…#1003) * Init all sponge tendency diag IDs to -1 immediately * No need to reset to -1 since initialized when declared * Move init_ALE_sponge_diags to after all tracers have been set up
These two references to members of a pointer don't seem to be hit except under special circumstances but nevertheless I ran in to them when debugging an unrelated problem. There are two references to members of `diag%axes` that assume `diag%axes` are associated, but in the specific case I was debugging this was not the case.
Five vertically integrated diagnostics are requested in CMIP7. These ultimately are to be for four vertical intervals (0-300m, 300-700m, etc.) but we will handle that through addition of a 4-level diagnostic grid, configured at run-time. This commit handles the conversion from temperature or salt to heat content or salt content (by mass) and registers a "vertically extensive" quantity so that the diagnostics know to re-integrate rather than remap. Changes: - Added diagnostics absscint, pfscint, scint, chcint and phcint - Moved registration of temp_int and salt_int to within an existing `if (use_temperature)` block - Made public 2 GSW conversion functions in MOM_EOS
…roducing_sum (and therefore, mpp_sum) is called. Previously, several 2-D arrays were each being passed within their own reproducing_sum calls, which is now avoided by consolidating the 2-D arrays into one 3-D array that is passed to a single reproducing_sum call.
…ux calculation. Needed for runs without frazil.
|
@jiandewang @chengz2 There were changes to values inside of masked land cells. This is the reason for the reported CI regressions. Do you know if there are any changes to solutions? (@chengz2 I assigned the NCAR review to @alperaltuntas but I can swap it to you if preferred.) |
|
@marshallward thanks for the explanation which solved my puzzle. Let me run UFS full regression suite and get back to you. |
|
For more details on the masking-based changes in the the Coriolis and advective accelerations at the boundary points, please see the conversation of the relevant commit at NOAA-GFDL#1029, and in particular the change to the masking that is applied at the end of |
|
Mea culpa, my comparator was too AGCM focused and didn't look at It looks like Also, |
|
this PR code didn't really change UFS answers although I will need to create a new baseline for UFS. The only change is two variables cau/cav in restart file and it didn't have any impact on all other model componets (atmos, wave, and ice). Even in ocean all other variables in restart files remain unchanged. In UFS we do cmp for restart files for all components, this is why I need a new baseline. |
|
We are finding that these changes do change solutions for our OBC pan-Antarctic configuration. I'll start digging, but any pointers would be much appreciated. The configuration is here: https://github.com/ACCESS-NRI/access-om3-configs/tree/dev-MC_panan4km_jra_ryf. ADDED: the |
|
@dougiesquire , thanks I am looking into this and some answer changes Inseong reported with the KMUG regional configurations. The CEFI/GFDL regional configurations are not changing answers (99% certainty). I used most of your settings in a test case, but there are a few differences I was not able to resolve (mainly because I do not have all the files needed): Internal tide dissipation is off, salinity restoring is off, a different ALE coord config file (but still z*, if you could share the ocean_vgrid.nc file that would be great), and INIT_LAYERS_FROM_Z_FILE = False. Do you see the answer changes even in runs that start from a restart instead of initial conditions? |
|
Thanks @theresa-cordero. You can download our vertical grid file here.
No - there are no answer changes when starting from a restart.
Unfortunately the pan-Antarctic configuration is our only regional configuration that's mature enough to be part of our routine testing. But we have got some other configurations under development that I could test. |
Great, that confirms my suspicion that the answer changes is probably because of differences in the OBC initialization with I will modify my test to use |
|
It turns out the cold/warm start behaviour and the fact the issue only occurred with our OBC config were both red herrings somewhat. I tracked the answer change down to 0f26d69 ("Add floor to "h_marg" in continuity_PPM"). This commit claims to not change answers by default, but that's only true if I'm happy to submit a PR with a fix if that's useful. |
mathomp4
left a comment
mathomp4
left a comment
There was a problem hiding this comment.
@sinakhani did a longer more comprehensive test and is happy. NASA approves.
|
@dougiesquire If you can submit a PR with a bug fix to this branch that would be great! Do we need a bug fix flag or do your answers stay the same with CONT_USE_H_MARG_MIN=True? |
I went to submit a PR but I see @Hallberg-NOAA has beaten me to it: NOAA-GFDL#1165 |
Set h_marg_min to about -1.057 lightyears when CONT_USE_H_MARG_MIN is False, to recover previous answers from the main branch of MOM6. Before CONT_USE_H_MARG_MIN was added, there was no floor on the marginal thicknesses, but the time that this parameter was added in github.com//pull/996 in December, 2025, a physically sensible floor of 0 was added for the marginal thicknesses for the continuity transports when it was not true. The marginal thicknesses is the thickness at a particular velocity in the sub-gridscale reconstruction that is integrated to get the time-integrated transport. However, there appear to be some cases where setting this floor to 0 changes answers from those with the main branch of MOM6. At this time we do not know how these negative marginal thicknesses might arise; perhaps they arise from CFL values that exceed 1 while iterating to reconcile the baroclinic and barotropic transports. However, resetting the floor to a huge negative value essentially ensures that no floor is applied, and empirically the answers from the main branch are recovered.
|
since @Hallberg-NOAA just added one more commit, let me re-run UFS regression test with focus on OBC cases |
|
The commit that we just added to this pull request should correct all known sources of answer changes relative to the current main version of MOM6. (Thanks go to @dougiesquire and Inseong Jang for working out this solution and verifying that it does restore the previous answers.) It turns out that the source of the answer changes arose from setting a floor of 0 on the marginal thicknesses used during reconciliation between the barotropic solver and the sum of the layer thicknesses in the continuity equation. We do not believe that there is a direct connection with the OBC code. Once again, we DO expect there to be changes in the restart files in the CAu and CAv fields at the coastlines, but these do not impact the exact reproducibility across restarts, nor do they change any solutions. |
dougiesquire
left a comment
dougiesquire
left a comment
There was a problem hiding this comment.
ACCESS-NRI approves
|
my USF regression with @Hallberg-NOAA's yesterday latest commit is done. Resutls are as expected and I will have a clear reason for UFS to have a new baseline. Thanks @herrwang0 for help me figured out the OBC setting and thanks @sanAkel and @DeniseWorthen for going over all the code changes line by line. |
Thank you @jiandewang |
20 participants
This merge includes 135 PRs from the NOAA-GFDL development branch. Major updates include open boundary conditions, ice shelf dynamics, new bottom drag parameterizations, ANN-based mesoscale parameterization, PLM-WLS reconstruction, diagnostics downsampling, and numerous bugfixes.
There are two potential answer-changing issues:
cuberootreplaces a non-reproducible power operator in CHANNEL_DRAG ifSET_VISC_ANSWER_DATEis20260704or newer. The default has changed to20260703to prevent unexpected answer changes, but users with explicit dates may be affected.Time differences (via FMS) were previously always reported as positive, which could introduce unwanted effects such as forcings running backwards in time. This is now fixed, but may impact solutions which relied on negative times being treated as positive.
Open Boundary Conditions (OBC)
Ice Shelf
flux_tmpNOAA-GFDL/MOM6#1096 (@OlgaSergienko)New Features / Parameterizations
Reconstructions / Numerics
Tracers
Bugfixes
Diagnostics
Runtime Parameters / Configuration
Build / Infrastructure
finalto avoid slow GCC builds" NOAA-GFDL/MOM6#1125 (@marshallward)finalto avoid slow GCC builds NOAA-GFDL/MOM6#1118 (@marshallward)Code Cleanup / Style
Documentation
Merge / Sync PRs