Support for time-varying input of ice sheet surface mass fluxes from file#1705
Closed
MJHarrison-GFDL wants to merge 164 commits into
Closed
Support for time-varying input of ice sheet surface mass fluxes from file#1705MJHarrison-GFDL wants to merge 164 commits into
MJHarrison-GFDL wants to merge 164 commits into
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 NOAA-GFDL#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 mom-ocean#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.
…mom-ocean#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.
`cuberoot()` uses the `modulo()` intrinsic to compute the exponent of the rescaled value of its input `a` in the range [0.125,1). However, there are platforms (e.g. NVIDIA GPUs) which do not implement `modulo()`. Review of this function also revealed some unnecessary complexity which was simplified by using the `sign()` intrinsic. This patch updates the `rescale_cbrt()` implementation, and also documents the underlying mathematics more carefully.
Converted the internal subroutine add_xyz_method() in MOM_diag_mediator into the function xyz_method() for greater generality and code-reuse. It no longer takes a diag_type argument, instead working directly with the contents of its axes_grp type argument. All answers are bitwise identical, and all publicly visible interfaces remain the same.
Set downsampled masks using the same stencils as is used for the masked downsampling of the underlying data. This required extensive revisions to the downsample_mask_2d() and downsample_mask_3d() routines, including a new method argument to indicate how the data is to be downsampled analogous to the argument that was already being provided to the downsample_field routines. All solutions and primary diagnostics are bitwise identical, but this commit corrects a glaring inconsistency in the masking of downsampled diagnostics.
Eliminated the MSK downsampling method and the MSK option code in the downsample_field() routines because this option did not take in to account the points that are actually used in downsampling fields that are not at tracer points. Instead the downsampled mask is determined by calling downsample_field() on an array of ones to determine where there are unmasked downsampled points. This change only applies in cases where an optional mask argument is being supplied to a post_data() call (as is the case for some static fields). As post_data() is currently used in MOM6, a mask argument only is passed in for diagnostics at tracer points so it is possible that there are not actually any instances where this will change any answers.
At one point, fill_miss_2d had an optional `prev` array that could be used to fill in values that were not horizontally connected to other points at the same level, and it issued a warning when that array was not provided. However, `prev` was always being provided, so in 2022 as a part of some more extensive refactoring and clean-up of MOM_horizontal_regridding.F90 this optional argument was made mandatory, and the extra logical test around the warning block became identical to the first. As noted in github.com/NOAA-GFDL/issues/929, this lead to a block of code that would clearly never be reached. This extra block of code has now been removed. All answers are bitwise identical.
Added an else branch in thickness_diffuse_init to set CS%Use_KH_in_MEKE to false when USE_MEKE is false, thereby avoiding the possible problems with segmentation faults due to an uninitialized variable, as documented in github.com/mom-ocean/issues/1696. It will remain a fatal error if `USE_KE_IN_MEKE` is in an input file but `USE_MEKE = False` and `FATAL_UNUSED_PARAMS = True`. All answers that worked previously will be bitwise identical, but this could fix problems with a segmentation fault arising from an uninitialized logical variable.
* Change loop order from m-[ij]-k to m-k-[ij]. * Rename u_L/v_L to L_in/L_out for clarity. * Add flux_to_res and face_area temporary variables for efficiency * Hoist dir and OBC%tres_[xy] writes outside the tracer loop. * Remove unused local variable ntr.
Copy resrv_lfac_in/out from segment%field (looked up by fd_id at runtime in update_segment_tracer_reservoirs) into the segment tracer registry. The copy can be populated via new optional arguments added to register_segment_tracer, which is only used by BGC tracers. resrv_lfac_in/out is used in subroutine register_segment_tracer, which updates tracer reservoir. Logically, segment%field is the "raw" data, internal to OBC module, while segment%tr_Reg is what is seen by the rest of the model. Storing resrv_lfac_in/out in the registry frees the need of segment%field data in update_segment_tracer_reservoirs. The change eliminates the need to use fd_id to identify whether resrv_lfac_in/out should be the default (1.0). Correspondingly, fd_index optional arguments is removed from register_segment_tracer, as its no longer needed. The change also paves the way to use resrv_lfac_in/out to fold the special uniform concentration case in the upcoming commit.
Modified the offline tracer code to always allocate the offline mixed layer depth when offline tracers are being used, and set it to HMIX_FIXED when READ_MLD is false. In addition, MLD_VAR_NAME is now passed to update_offline_from_files(), whereas it was previously just being set and ignored. As a follow up, it might be useful to add the option of setting the mixed layer in offline tracer mode dynamically using diagnoseMLDbyEnergy() or diagnoseMLDbyDensityDifference(), but that is beyond the scope of this commit. This commit addresses the bugs identified at github.com/NOAA-GFDL/issues/1100. All answers are bitwise identical in all cases that worked previously, but some offline tracer cases that had been failing previously with a segmentation fault should work properly now.
…nservation, and adding time integration (mom-ocean#1088) * Add ability to use different ML depths with brine plume Codes in 3 options, MLD_003, MLD_EN1, and H_EPBL, all of which can be passed to brine plume param to set the plume depth * Update brine plume scheme to conserve salt and time integrate - Brine plume scheme failed dimensional consistency because it was not multiplied by the time step. The units of plume_flux in the description are updated, which exposes why it failed the time dimension consistency test. - The brine plume scheme was moved after the application of other fluxes that can add salt/heat/mass to the ocean in applyboundaryfluxesinout. This change made it easier to track that it only moves salt in the vertical and doesn't create/destroy salt. - The brine plume salt flux is now multipled by dt to convert from a rate to an amount. This has a large impact on the behavior, but is approximately compensated by an unrealated bug in SIS2. - The brine plume scheme did not conserve salt because it could neglect salt if the mixed layer depth wasn't deeper than the center of the bottom grid cell, possibly due to round off. A new condition is added to insert salt into the bottom layer if it is thick enough to accept salt. This significantly improves conservation, though round-off level errors are still possible. - Complete a stub for a brine plume diagnostic, which can be used to verify how the brine plume scheme moves salt in the vertical. - Adds an option to increase the vertical scale for brine plumes proportionally to the chosen mixed layer depth. - Updated some parameter names and formatting related to the MLD used in the brine scheme. * Fix OMP directives for brine plume helpers in MOM_diabatic_aux * Updates for conservation in brine plume parameterization - The brine plume salt flux contribution is no longer added to the top level and subtracted later in the brine plume part. Rather, the salt that is added in applyboundaryfluxes now has the brine plume salt flux removed. The salt flux is then simply applied by adding it back in during the brine plume part. - The code to handle flux deposition within the mixed layer was rewritten to ensure the salt flux is entirely used by the bottom-most finite thickness layer. - The flux is now prescribed by analytical integral over the layer instead of midpoint quadrature. - Debugging options are enhanced to better track and report on salt conservation. * Make sure salt added/removed in brine plume is always computed * Fix formatting and OMP directive * Clean up brine plume commit, standard out, and debug options * Fix OMP associated with brine plume * Fix OMP associated with brine plume (again) --------- Authored-by: brandon.reichl <brandon.reichl@noaa.gov>
This commit provides the MOM6 diag_mediator with information about the locations and orientations of any open boundary condition segments at velocity points, and then uses this to select the correct thicknesses to use for vertical remapping of diagnostics onto alternative diagnostic grids (e.g., z or rho2). It also makes a change to properly interpolate the thicknesses that are used as weights for downsampling velocities. Previously a one-sided estimate of the thicknesses was being used, which did not satisfy rotational consistency. The specific changes here include: . The addition of the new publicly visible routine diag_mediator_set_OBC_info() and a call to this routine from initialize_MOM() when OBCs are being used. This routine stores information about the location and orientation of OBC segments. It does not need to be called when OBCs are not being used. . The addition of two new integer arrays to the diag_ctrl type with values of 0, -1 or 1, depending on the presence and orientation of OBCs at velocity faces. . The addition of OBC_u and OBC_v arguments to 3 publicly visible routines in MOM_diag_remap (diag_remap_do_remap(), vertically_reintegrate_diag_field() and vertically_interpolate_diag_field()) and to their private counterparts (do_remap vertically_reintegrate_field() and vertically_interpolate_field()). These arrays are used to select the appropriately interpolated thicknesses at velocity points to use as weights for remapping. In tests of these changes in the loop current test case, the changes to the remapped velocities and transports at OBC points are subtle but real, while in the interior (away from OBCs) the diagnostics are unaltered. All solutions are bitwise identical and no primary diagnostics are changed, but there are changes to remapped diagnostics in cases with open boundary conditions. Within this commit this a new publicly visible subroutine, several changes to the interfaces of 3 publicly visible routines and the addition of two new 2-d arrays to the diagnostics control structure.
* *Add framework support for ANN mesoscale streamfunction Extend MOM_isopycnal_slopes to output density gradients (drdx_u, drdy_v, drdz_u, drdz_v) as optional arguments from calc_isoneutral_slopes. Extend MOM_thickness_diffuse to accept externally computed streamfunctions (Sfn_unlim_u_3D, Sfn_unlim_v_3D) in thickness_diffuse_full. Wire up the new ANN streamfunction module in MOM.F90 with initialization, computation, and restart calls. * *Add ANN-based mesoscale streamfunction module Introduces MOM_meso_sfn_ANN, which predicts density-flux-based mesoscale streamfunctions using a neural network. Features include configurable stencil windows for local density and velocity gradient inputs, runtime-tunable clamping thresholds, boundary proximity masking, and a fallback GM_simple mode for testing. The ANN outputs buoyancy fluxes that are converted to streamfunctions via the local density gradient magnitude. This is answer-changing only when USE_THICKNESS_DIFFUSE_ANN=True. The parameter defaults to False, so no existing experiments are affected. * *Clean up ANN mesoscale streamfunction review feedback Address Pavel Perezhogin's review and failing CI checks. Review: - Swap ANN_apply_vector_orig -> ANN_apply_vector_oi (faster path) - Gate MOM_thickness_diffuse limiters on use_meso_sfn_ANN so the dev/gfdl baseline is unchanged when the ANN is off - Refactor the streamfunction clamp to act on the velocity-scale Upsilon (Ferrari et al. 2010) before multiplication by dy_Cu/dx_Cv, giving a grid-independent cap. Rename the runtime parameter MESO_SFN_SFN_CLAMP to MESO_UPSILON_CLAMP and shift the default from 1.0e6 volume-transport units to 1.0e4 m2 s-1 velocity-scale. Unit cleanup: - Density flux variables (Fx_u, Fy_v, Fx_c, Fy_c) relabeled from [R L-1 ~> kg m-4] to [R L T-1 ~> kg m-2 s-1] - Volume streamfunction dummy args sfn_u/sfn_v relabeled [L2 T-1] -> [Z L2 T-1 ~> m3 s-1] to match Sfn_unlim_u_3D - mag_grad / mag_grad_floor switched from [R L-1] to [R Z-1] (the computation already produces [R Z-1] via the Z_to_L rescale of drdx_u; the floor comparison only works correctly under internal rescaling when both sides share the same internal unit) - Fix conversion= factors on all 11 meso_sfn_* diagnostic fields; previous values were copy-pasted from sfn_unlim_x and did not match the registered quantities - Add scale= argument to mag_grad_floor and MESO_UPSILON_CLAMP get_param calls for rescale-safe behavior - Replace MESO_SFN_GRAD_MAG_EPSILON parameter with derived h_neglect-style locals rho_grad_neglect [R L-1] and vel_grad_neglect [T-1] (1e-30 scaled per dimension) CI fixes: - build-openmp: add present_drdx_u/drdx_u/drdz_u and present_drdy_v/drdy_v/drdz_v to the two !$OMP parallel do shared(...) clauses in calc_isoneutral_slopes - check-style-and-docstrings: split multi-variable id_* declarations in MESO_SFN_ANN_CS so each has its own doxygen !< comment Answer-changing only when USE_THICKNESS_DIFFUSE_ANN=True. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * *Audit cleanup and CI fixes on meso_sfn_ANN Act on the broader audit of MOM_meso_sfn_ANN and fix the two remaining CI failures. Audit-driven changes: - Remove the GM_simple and GM_ann model-type branches entirely; they were testing scaffolding. The only production path is nondim_rhoF_ann. This also removes the MESO_SFN_ANN_TYPE runtime parameter, the meso_sfn_ann_model_type CS field, and the now-unused local Kappa. - Add a startup bounds check for MESO_SFN_ANN_WINDOW: the stencil shift must be <= 3 to stay inside the halo=3 provision in meso_sfn_ANN_compute, and the window must be an odd integer. - Set and check the CS%initialized flag so using the module without calling _init produces a clear FATAL instead of later undefined behaviour. - Rename lowercase runtime params to the MOM6 uppercase convention: meso_sfn_ann_coeff/type/window/file -> MESO_SFN_ANN_COEFF, MESO_SFN_ANN_WINDOW, MESO_SFN_ANN_FILE (type removed). - Expand the module docstring with a Ford-style summary of the pipeline and an explicit warning that the dimensionalization in _compute is an implicit contract with the training procedure. - Document the extended stencil bounds in the corner-point velocity gradient loop. CI fixes: - build-openmp: add Sfn_unlim_u_3D / Sfn_unlim_v_3D to the shared(...) clauses of the two !$OMP parallel do blocks in thickness_diffuse_full. - check-style-and-docstrings: wrap the Fy_c declaration comment onto a second line so the line stays under the 120-char limit. Breaking parameter changes (branch not yet released): MESO_SFN_ANN_TYPE removed; the four lowercase params listed above renamed to uppercase. Answer-changing only when USE_THICKNESS_DIFFUSE_ANN=True. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Resolve memory leak (deallocate) Use faster ANN inference interface (ANN_apply_array_sio). 18% acceleration, regression is the same * Changes to pass MOM6 testing system. The tests passed: test.grid test.layout test.restart test.openmp test.nan test.dim 1) Trailing space is removed 2) Physical dimensionality is added to all float variables 3) Loop bounds are adjusted to get rid of out of bound memory access and pass test.grid 4) Scaling with physical units (i.e., structure US) is adjusted in multiple places to pass test.dim 5) Maximum allowed window size is limited to 3, which was tested now. Larger windows will likely needed calling pass_var (not yet implemented) * Add proper interaction of buoyancy ANN with use_stored_slopes and SKEBs parameterizations * Change default value of MESO_UPSILON_CLAMP parameter Default values of limiting parameters are adjusted based on the diagnosis of CESM ocean model output. MESO_UPSILON_CLAMP reduced from 1e4 m2 s-1 to 15 m2 s-1. This choice guarantees that the maximum overturning per one grid point and water column is approximately 1 Sv on CESM grid (before limiting and smoothing performed in MOM_thickness_diffuse.F90). The previous default value could allow up to ~600 Sv per grid point, which is an enormously large number. MESO_SFN_FLUX_CLAMP, which is ~100 kg m-2 s-1, is likely never activated. Still, this number is reasonable as maximum predicted fluxes in CESM are ~20 kg m-2 s-1. MESO_SFN_MAG_GRAD_FLOOR, which is 1.0e-10 kg m-4, is a very realistic number for CESM. The density gradients are approximately in a range from 1e-10 up to 1e-1, with a couple of grid points having 1e-12. * Integrate pass_vector for sfn_u and sfn_v Adding an exchange of halo points for the stream function predicted by the neural network. This halo exchange has no effect in simple geometries, but gaurantees Temperature/Salinity conservation on a tripolar grid, where rotational invariance is required to achieve conservation through the pole. Exact rotational invariance is not supported by the neural network. Rotational invariance is learned from data approximately, and is preserved less strictly than the numerical precision. Exchange of a streamfunction guarantees that there is no jump in the thickness flux at the grid point adjacent to the pole. * Encapsulate use_meso_sfn_ANN in thickness_diffuse_CS Address Hallberg's review comments about opaque-type encapsulation (NOAA-GFDL#1109). - MOM.F90: Drop the if/else around the two thickness_diffuse calls and pass u, v unconditionally. The previous gating on CS%thickness_diffuse_CSp%use_meso_sfn_ANN peeked inside the opaque thickness_diffuse_CS type. The runtime check inside thickness_diffuse at line 503 (CS%use_meso_sfn_ANN .and. present_vel) still gates the ANN call correctly; present_vel is now always true but the conjunction is harmless. - MOM_thickness_diffuse.F90: Make use_meso_sfn_ANN private (was incorrectly declared public). No external readers remain after the MOM.F90 refactor. Expected to be bitwise identical based on inspection of the runtime gate at MOM_thickness_diffuse.F90:503 (the conjunction CS%use_meso_sfn_ANN .and. present_vel still routes the same way). Not locally verified by test.repro / test.regression; CI will check. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Clarify drdx_u/drdy_v docstrings as gradients along z surfaces Per Hallberg's review feedback on NOAA-GFDL#1109: emphasize in the Doxygen comments that the new drdx_u and drdy_v output arguments of calc_isoneutral_slopes are evaluated along surfaces of constant z (not along isopycnals or model interfaces). All answers are bitwise identical. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * *Apply MOM6 soft-case convention in MOM_meso_sfn_ANN Address Hallberg's style review on NOAA-GFDL#1109: - Loop variables: rename lowercase k -> K on interface-range loops (k=2,nz and k=1,nz+1) in meso_sfn_ANN_compute, center_grad_rho, and center2uv. Rename lowercase j -> J on v-point loops (j=js-1,je with v-point array access in the body). Update one inner u-loop (do i=is-1,ie -> do I=is-1,ie) to match the body's I usage. - Array indexing: flip interface-staggered array accesses (drdx_c, drdy_c, drdx_u, drdy_v, drdz_u, drdz_v, Fx_c, Fy_c, Fx_u, Fy_v, sfn_u, sfn_v, var1_c, var2_c, var1_u, var2_v) to use uppercase K on their third index. Layer-staggered velocity arrays (dudx, dudy, dvdx, dvdy) keep lowercase k to signal "the layer below interface K". - mag_grad expression: rewrite from the FMA-paranoid (US%Z_to_L*drd_*)*(US%Z_to_L*drd_*) form to Hallberg's preferred US%Z_to_L**2*drd_***2 form (he explicitly noted the extra parentheses are not needed for FMA/symmetry). Answer-changing only when USE_THICKNESS_DIFFUSE_ANN=True. The case-only changes are bit-identical by language definition (Fortran is case-insensitive). The simplified mag_grad expression is bit-different under FMA-aware compilers. Since the feature defaults to .false. and has not been released, no _BUG flag is required. Not locally verified; CI will check. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Drop optional from u,v in thickness_diffuse; fix v docstring Address Hallberg's follow-up comments on NOAA-GFDL#1109: - Drop `, optional` from the u and v dummy arguments of thickness_diffuse. Both callers in MOM.F90 unconditionally pass u,v after the previous opaque-type refactor (98038be), so the optional attribute is no longer needed. - Remove the now-redundant `present_vel` variable (declaration and assignment), and simplify the ANN call gate from `if (CS%use_meso_sfn_ANN .and. present_vel)` to `if (CS%use_meso_sfn_ANN)`. - Restore the truncated `v` Doxygen description and tighten whitespace on both u and v declaration lines. All answers are bitwise identical: `present_vel` was always true after the prior refactor, so the conjunction always took the same branch. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Make find_eta halo size conditional on use_meso_sfn_ANN Per Hallberg's review on NOAA-GFDL#1109: avoid the cost of the wider halo in find_eta for default users. The ANN streamfunction needs the interface elevation `e` with a wider halo to support its stencil. Previously the halo bump was unconditional. Now it is gated on `CS%use_meso_sfn_ANN` so that default users (USE_THICKNESS_DIFFUSE_ANN=False) retain the original `halo_size=1` call. Default-off users now match upstream dev/gfdl exactly at this call site. ANN-on users are unchanged relative to the previous branch state. Hallberg also recommends verifying halo updates by comparing runs at different MPI process counts (e.g., mpirun -n 63 vs -n 64); that test has not yet been run. CI may catch some issues; the cross-process-count check is a separate validation step. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com> Co-authored-by: Pavel <pperezhogin@gmail.com> Co-authored-by: Pavel Perezhogin <35234405+Pperezhogin@users.noreply.github.com>
Building MOM_diag_mediator.F90 with GCC and -O2 became much slower after
diagnostic buffers were added.
```
commit 6f6975a
Author: Andrew Shao <andrew.shao@hpe.com>
Date: Wed Mar 18 12:17:43 2026 -0700
+Extend diag_mediator to allow the piecemeal posting of diagnostics (mom-ocean#809)
```
Before this commit, MOM_diag_mediator.F90 built in about 15 seconds.
After this commit, build time increased to more than 170 seconds.
The slowdown was isolated to derived types with multiple allocatable
array components, not to abstract types, inheritance, or type-bound
methods. A minimal form of the triggering pattern is shown below.
```
type :: diag_buffer_2d
real, allocatable :: buffer(:,:,:)
integer, allocatable :: ids(:)
end type
```
The component did not need to be referenced in MOM_diag_mediator.F90.
It was enough for diag_buffer_2d to appear as a component of axes_grp.
Somehow, this was resolved by defining an explicit dummy finalizer in
`diag_buffer_2d`.
```
type :: diag_buffer_2d
! ...
contains
final :: finalize_diag_buffer_2d
end type diag_buffer_2d
subroutine finalize_diag_buffer_2d
type(diag_buffer_2d), intent(inout) :: this
end subroutine finalize_diag_buffer_2d
```
With this change, build time returns to about 15 seconds.
* The finalizers are intentionally empty. Allocatable components are
automatically deallocated by the language, so no explicit deallocate()
calls are needed here.
Finalizers should only be used for custom cleanup, such as external
resources or manually managed pointer targets.
* `diag_buffer_[23]d` contains an array of `buffer_[23]d` objects. The
language standard notes that components are finalized before the
parent.
* No similar compile-time benefit was seen from adding dummy finalizers
to buffer_[23]d, so those types are left unchanged.
Previously, the nonlinear part of basal friction was evaluated cell-wise, while the linear part was evaluated at quadrature points; the resulting inconsistency in the Jacobian prevented quadratic convergence when using Newton iterations, so moved all friction to be evaluated evaluated at quadrature points instead. Viscosity should also be evaluated at the 4 quad points rather than cell-centered, so changed the default for NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS from 1 to 4 (which is more accurate anyway). These changes allow us to achieve the expected quadratic convergence when using Newton. With these changes, several deprecated cell-averaged ice-shelf CS fields were removed: basal_traction, newton_umid, newton_vmid, newton_drag_coef. Subroutine calc_shelf_taub was converted to return area-averaged basal shear stress, now for diagnostic purposes only. Also removed incorrect 0.5 factor from newton_visc_factor
-Added NEWTON_CONJUGATE_GRADIENT_TOLERANCE to set a separate tolerance for Newton iterations compared to Picard. Helpful when using a higher convergence tolerance for initial inexact Newton tolerances, or a lower tolerance is used for "exact" Newton. -Can now set NEWTON_AFTER_TOLERANCE<=0 to only use Newton iterations to solve the SSA, bypassing Picard iterations entirely. This works well if the initial guess is close to the solution. Further developed the Inexact Newton Eisenstat-Walker (EW) scheme: -EW now evolves (correctly) based on the change in L2 norm of the stress residual -New parameters NEWTON_EW_GAMMA and NEWTON_EW_ALPHA to specify the EW gamma and alpha instead of hardcoding them -Added new EW safeguard based on PETSc, choice 3, along with parameters: -NEWTON_EW_SAFETY: choose between the 2 safeguards -NEWTON_EW_1_THRES: threshold used for when NEWTON_EW_SAFETY is 1 -NEWTON_EW_ETA_MAX: Max allowed EW eta (between 0 and 1) All of these new parameters will typically work well at their defaults, but are available for fine-tuning Added new SSA convergence criteria based on L2 norms (needed for EW): -NONLIN_SOLVE_ERR_MODE=4: err_mode for convergence based on nonlin residual | F | / | F_0 | -NONLIN_SOLVE_ERR_MODE=5: err_mode for convergence based on relative residual | F | / | tau | -SSA_ADD_REL_RESID = True: If NONLIN_SOLVE_ERR_MODE<=4, this will also ensure convergence of the relative residual (via err_mode=5). While generally useful for ensuring convergence of the solution (err_mode 2 or 3) and the residual simultaneously, in some cases, it is essential. For example, if using a solution-based convergence criterion with inexact Newton, initial Newton iterations may not change the solution much and could trigger the stopping criterion before the problem is actually converged; using SSA_ADD_REL_RESID avoids this issue. -ICE_RR_NONLINEAR_TOLERANCE: Tolerance for err_mode 5 if SSA_ADD_REL_RESID=True. Added several new MOM_mesg calls to track the adaptive inner solver tolerance and the nonlinear relative residual. Added VERBOSITY parameter to ice-shelf driver Ice-shelf efficiency improvements - Precomputed some commonly-used variables - Replace some pass_var()/pass_vector() pairs with do_group_pass() - Indentation Add missing ice_shelf CS%rhow_rhoi substitutions
Set a pointer to the `df_2d_y` optional argument in `register_tracer()`, treating it like `df_2d_x` and the other diagnostic arrays that can be provided as optional arrays to `register_tracer()`. Any 2-d y-diffusion diagnostics impacted by this bug would have reported all 0 values, but now will return sensible values. Almost all (perhaps all) tracers have their standard diagnostics set up via `register_tracer_diagnostics()`, and this bug fix does not impact them. As a result, it is entirely possible that this bug does not impact any code that is used, which is why it was never detected. Additionally, this commit comments out two unused and unimplemented diagnostics (ear and ebr) from MOM_offline_main. For now the registration calls for these diagnostics are being left in comments in case someone actually intends for them to exist, but if they are not actually implemented they should be removed altogether at a later date. These broken diagnostics were detected while scanning for `register_diag_field()` calls that appear to be missing appropriate unit conversion factors. Without an actual implementation, there would have been an error message about missing `post_data()` calls for these diagnostics had they ever been turned on in the diag_table of an offline tracer run, so they can not have been in use in any working simulations. All solutions are bitwise identical, but there will be changes in the available_diags files for offline tracer runs.
When ODA_INCUPD=True, the NOTE "ended updating fields with increments. " was intended to be produced once per run but was instead being printed every time step after the increments were done. The NOTE isn't necessary and so has been deleted. Answers are unchanged.
mom-ocean#1115) * Add capability to use multiple diagnostics downsampling levels in MOM6 - This commit presents an enhancement of diag downsampling workflow for MOM6 - It allows an arbitrary number of downsampling factors in the same run (if they are commensurate with grid processor layout) - The number of downsampling levels should be specified in MOM override parameters E.g., #override NUM_DIAG_DOWNSAMP_LEV = 3 - The downampling factors should be specified E.g., #override DIAG_DOWNSAMP_LEVS = 2,4,5 - The required downsampled domains and axes are then generated at the initialization by MOM_diag_mediator No need for other modules to know anything about the downsampling scheme. - No need for a #define MAX_DSAMP_LEVEL, everything is allocated at runtime - Since everything is allocated at runtime users need to be congnisent of the number of levels they request. - Any downsampling diagnostics can be put in diag_table , only the levels that requested in parameters are handled - The overall dependence of downsampling levels on the grid shape and pe layout remains Only levels should be requested where NIGLOBAL/x_layout and NJGLOBAL/y_layout are divisible by downsampling factor
The time reference for CS%dt_obc_seg_time was initialized from Time (the start of the current run segment) rather than Time_init (the overall simulation start time). After a restart, Time > Time_init, so the first scheduled update could happen at a different phase than it would in a continuous run, breaking restart reproducibility. The fix anchors dt_obc_seg_time to Time_init so that the update schedule is invariant to where in the simulation a restart occurs. The old behavior is preserved by OBC_BGC_TIME_REF_BUG.
When DT_OBC_SEG_UPDATE_OBGC exceeds DT_TRACER_ADVECT, BGC OBC external data (%t) read at initialization cannot be recovered after a restart because it is not saved in restart files. This means the first tracer advection step after a restart may use incorrect boundary data. Issue a FATAL error when this condition is detected on a restart run, and a WARNING on a new run to alert the user that any future restart will fail. Both messages direct users to set DT_OBC_SEG_UPDATE_OBGC to DT_TRACER_ADVECT or 0 until %t is added to restart files.
This commit adds as edge case for tracer length scale at OBC segments.
Runtime parameters OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN/OUT provide
length scales (at domain level, identical for all segments) for tracer
advection and the scale can be positive/finite and zero (i.e.
instant values from interior/externally prescribed, which is the
default). A new case that allows infinite length scale is now added, in
which case the reservoir is frozen. This change helps folding the
OBC_inflow_conc cases (two users cases) into the reservoir update.
* OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN/OUT can now be negative, which is
a sentinel as infinite length scale.
* Internally, Tr_InvLscale is reclassified in a more mathematically
consistent manner:
* For finite length scale, Lscale > 0, Tr_InvLscale = 1 / Lscale.
* For zero length scale, Lscale == 0, Tr_InvLscale = -1.
* For infinite length scale, Lscale < 0, Tr_InvLscale = 0.
* Per-tracer resrv_lfac_in/out works as a tracer specific factor, but
can also override the per-segment Tr_InvLscale in the edge cases.
* resrv_lfac = 0 -> infinite length scale
* resrv_lfac = -1 -> zero length scale
* Add pre-computed per-tracer signed inverse length scales
I_Lscale_in/out to OBC_segment_tracer_type. The sign encodes the update
mode:
* positive = finite length scale
* zero = infinite (reservoir frozen)
* negative = instant-update sentinel.
Other changes:
* Replace the segment-level b_in/b_out switches in
update_segment_tracer_reservoirs with per-tracer mask_a_in/out and add
mask_L_in/out, both derived from I_Lscale.
* Update OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN/OUT parameter descriptions
to document the three-way (finite/infinite/instant).
There are currently two treatment of OBC tracers. 1) Use the full 2D
array of %t and %tres which goes through a per-advection step reservoir
update through update_segment_thickness_reservoirs. 2) Use a simple
constant scalar OBC_scalar in place of the reservoir.
The scalar method is only used by a couple of user cases, but it demands
an if-branch in MOM_tracer_advect and complicates the code structure.
This commit incorporates the scalar method into the full array route.
OBC_scalar is populated into the reservoir array (%tres).
resrv_lfac_in/out are set to zero (infinite distance to the reservoir
from both interior and exterior), so that %tres does not update from
update_segment_thickness_reservoirs.
The old OBC_inflow_conc method is realized by providing an optional
input to register_segment_tracer.
* The t and tres arrays are now unconditionally allocated at
registration in register_segment_tracer.
* In register_segment_tracer
* Remove the OBC_array and fd_index optional arguments from
register_segment_tracer.
* Add a FATAL error when OBC_scalar and resrv_lfac_in/out
are both provided.
* Remove the misused temp/salt_segment_data_exists flags from
OBC_segment_type.
* Remove fd_index from OBC_segment_tracer_type and
OBC_segment_thickness_type. It was used to find resrv_lfac_in/out.
* Add allocated(OBC%tres_[xy]) guards in radiation_open_bdry_conds
* Since tres is now always allocated at registration, remove the
OBC_inflow_conc fallback branches from MOM_tracer_advect and
use tres directly without the allocated() guard.
* Remove OBC_inflow_conc from tracer/thickness types
* Replace loop range for tracer m from OBC%ntr (total number tracers)
to segment%tr_Reg%ntseg (per-segment number of OBC tracers). Currently,
they don't necessarily equal.
* Remove unused Tr attribute from OBC_segment_tracer_type
* Condense the register_segment_tracer header and tr_ptr argument description. * Add a detailed header to update_segment_tracer_reservoirs describing the backward-Euler implicit update, the three I_Lscale regimes (frozen, finite length scale, instant update), and the OBC%tres_x/y restart output. * Add inline comments above the I_Lscale assignment block explaining three regimes and the sentinel clamping rules. Assisted by Claude Sonnet 4.6
This reverts commit 56001bc. The dummy `final` functions in `MOM_diag_buffers` cause segmentation faults in NVIDIA. Rather than work around a potential optimization in GCC, the safer option is to remove this work and revisit it later if the problem is not resolved on the compiler side.
Add an optional Time argument to set_dtbt() so that each call emits a verbosity-3 message on the root PE recording the model date, dtbt, and dtbt_max at the moment of the reset. Pass Time to set_dtbt() from barotropic_init(), step_MOM_dyn_split_RK2(), and step_MOM_dyn_split_RK2b(). The RK2 callers pass Time_local minus one baroclinic dt to report the time at the beginning of the step.
This commit address a potential reproducibility issue across restarts for some niche cases. Previously, DTBT is always recalculated at the first dynamics time step, ignoring DTBT from restart files. The reason is that the initial CS%dtbt_reset_time is set to the most recent multiple of the interval before or equal to current time. So the check (Time_local >= CS%dtbt_reset_time) [Time_local is the end of the step] in step_MOM_dynamics is always true at the first step. In addition, the code for backing up initial CS%dtbt_reset_time by one interval is dead code. The variable calc_dtbt return from barotropic_init() was never used (and it was incorrectly left as .true. for restart runs with a DTBT saved). Add a runtime paraemeter DTBT_RESTART_BUG (default=True) to fix the bug and preserve answers. Fix calc_dtbt in barotropic_init() to be .false. whenever DTBT was read from a restart file. When .false., dtbt_reset_time is advanced by DTBT_RESET_INTERVAL after the current time, so the first step uses the restart value unless a reset is naturally due. Also normalize DTBT_RESET_PERIOD to 0 when it is <= DT, since any period that short triggers a reset every step anyway and the timer logic adds no value.
1) Set opposing ocean pressure to zero if the ice is grounded above sea level. 2) Corrected the ice-front Neumann BC to use cell edge lengths instead of dxT and dyT, in alignment with the finite element SSA solution.
This file contains hidden or 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
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
19 participants
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.
These changes are introduced for the ice sheet solo driver and enable the specification of time-varying surface mass input. Setting MOM_parameter: ICE_SMB_TIME_VARYING = True enables this option.