Latest dev/gfdl#12
Merged
Merged
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.
t17d and t20d are the depth of the 17 and 20 degC isotherms respectively, measured from the surface. Interpolation is done via Recon1d reconstruction and then a new function in the class that solves that inverts the reconstruction. This solver was not in the original class because I had originally only thought it would be used in the new grid generator which interpolates density and not for the reconstructed field. Since some other diagnostics for CMIP will make use of remapping or T, I figured we might as well use the reconstructions to do the interpolation. - adds to new registrations and posts in MOM_diagnostics - adds function `%x(self, k, val)` to the Recon1d class.
Revised the MOM_ice_shelf_diag_mediator code that creates the MOM_IceShelf.available_diags file to follow the same pattern of logging as the rest of MOM6. The specific changes include moving the module name(s) to their own line, adding static diagnostics to the variables that are logged in the available_diags file (however, the ice shelf code does not have any static diagnostics yet), logging the axes of the diagnostic, and adding cell metrics to describe how a diagnostic might be reduced and logging these as well. All solutions are bitwise identical, but there are new entries and added information in the MOM_IceShelf.available_diags files.
Modified set_IS_axes_info() to follow the pattern from the rest of MOM6 and take the units for the horizontal axes from the grid type. In most cases, this has the effect of changing the units from "degrees_E" and "degrees_N" to "degrees_east" and "degrees_north" to follow CMIP protocols, and to mirror changes that were made to the standard MOM6 framework code in 2017. The names of the gridspace axes were changed from `xB`, `yB`, `xT` and `yT` to `Iq`, `Jq`, `ih` and `jh` to mirror the names in MOM6. The long descriptions of these array axis variables were also updated for clarity when `USE_INDEX_DIAGNOSTIC_AXES` is true. All solutions are bitwise identical, but there are changes to the metadata and axis variable names in some output files, and some runtime parameters that had been read from both `set_IS_axes_info()` and `set_grid_metrics()` in MOM_grid_initialize.F90 are now only read in the latter routine. As a result of these changes, `set_IS_axes_info()` no longer takes a `param_file_type` argument and the grid argument is now `intent(in)`.
Added a missing dimensional scaling argument in the test for when to write harmonic analysis diagnostics in HA_accum. Also used the new unscale argument in 2 calls to real_to_time in HA_init and used the new scale argument in one call to time_to_real in HA_accum. Answers and diagnostics are bitwise identical when scaling is not applied, but this does correct a problem where some diagnostics may have been written at the wrong time previously when dimensional rescaling of time is being used.
Edited the conversion factors for the 'ale_u2' and 'ale_v2' diagnostics to use an equivalent form that makes it a little more obvious that the conversion factors correspond to the declared units. All answers are bitwise identical.
This commit clarifies a number of comments in the MOM_diag_mediator. It also corrects multiple instances of spelling errors or non-standard use of white space in comments. There was also some rearrangement of module use statements and the elimination of other unneeded module use calls. Only comments or spacing are changed, and all answers are bitwise identical.
A locally allocated array was uninitialized which is harmless for most configurations that use time-filtering in the vertical grid generation. The time-filter uses land masks and sets things to zero on land. Without the time-filter, we encountered non- reproducible behavior on the root PE. I also de-allocated a local variable (Rb) used on another if-branch in the same routine.
This commit refactors update_OBC_segment_data using the previously
introduced new structure segment%field component. Specifically, the code
that calculates normal_vel, normal_trans, normal_trans_bt,
tangential_vel gradient and external values of tracer / thickness
reservoir, is updated. The changes include,
* The previous loop + if structure is replaced by direct indexing of the
physical fields.
* Tidal velocity and SSH code is refactored to cji structure (c for
constituent), rather than an embedded if-branch within [ji] loop of a
[ij]c structure.
* As a result, segment has three new 2D arrays, tidal_v[nt] and
tidal_elev.
* Introduce orientation-agnostic looping indices [ij][se]_seg, currently
only used for tracer cell fields.
* Add tr_index in OBC_segment_data_type (segment%field) to facilitate
reservoir update.
* Add Newton iterations for the ice-shelf velocity solution The previous Shallow Shelf Approximation solution used an iteration-on-viscosity scheme with Picard iterations only. This commit adds the capability to change to Newton iterations once a certain tolerance is met, which should reduce the number of iterations needed to reach convergence. This tolerance to change to Newton is set by the new parameter NEWTON_AFTER_TOLERANCE. This parameter defaults to the same value of ICE_NONLINEAR_TOLERANCE so that Newton is not called by default. The Newton scheme is paired with the Eisenstat & Walker 1994 approach for adapting the CG tolerance to help the Newton scheme achieve quadratic convergence without over-tightening at later Newton steps. To use the adaptive CG tolerance, set NEWTON_ADAPT_CG_TOLERANCE=True (default is True). * Newton iteration PR fix Style changes, parentheses, removed string comparison from loops * Pass_var optimization and FMAs Made some changes to pass_var 'complete' arguments for Newton-related variables. Also removed some Newton-related pass_vars and a CS%doing_newton=.false. that were made before entering the the SSA solution loop (at this point, Newton is already false, and there is no need to update Newton-related halos yet). Added parentheses for rotational symmetry with FMAs.
Replaced "end if" with "endif" in 35 places in 10 files and "end do" with "enddo" in 12 places in 4 files to comply with the MOM6 style guide. In addition, in 23 places in 7 files, code like `if(test)` were replaced with `if (test)`, also to align with the MOM6 style guide. All answers are bitwise identical.
Removed 41 trailing semicolons (MOM6 is in Fortran, not C). In 406 other lines of code the white space around semicolons was standardized to ' ; ' to match the other ~14000 lines of code with semicolons joining Fortran statements on the same line. In a few cases there were also some improvements to indents related to these changes. With these changes, all semicolons that could be replaced with a newline and still have correct Fortran follow the ' ; ' pattern. These changes are scattered across 86 files. Most of these changes involve only white space, and all answers are bitwise identical.
Create a new subroutine read_OBC_segment_data out of a very chunky update_OBC_segment_data. The new subroutine enables better control on whether the update is associated with new external data through IO or just recalculation (with evolved model internal state, i.e thickness, tides). Misc * Use orientation-agnostic loop indices for segment%h and segment%htot * Add [ij]_offset_in to loop over interior cells * Field segment%Cg, which is never used, is now removed.
Create a new subroutine initialize_OBC_segment_reservoirs to initialize thickness and tracer reservoir of OBC segments. This part of code was original part of update_OBC_segment_data, but it is only called during initialization.
This commit fixes a bug in particle advection of particles that are advected by uh and vh. When the tracer timestep doesn't match the dynamics timestep, these particles need to be advected using the tracer timestep. The previous version of the code sometimes leads to unphysical particle movement, and this commit fixes the issue.
Added the new interfaces `symmetric_sum()` and `symmetric_sum_unit_tests()` to the MOM_array_transform module, and added a call to `symmetric_sum_unit_tests()` to `unit_tests()`. The two functions wrapped by the `symmetric_sum()` interface are intended to be a general template for rotationally symmetric sums, but they demonstrably work as intended (as demonstrated by the passing unit tests) for any rectangular arrays. The new driver in test_MOM_array_transform.F90 tests these new arrays. Although `symmetric_sum()` is not used yet in the MOM6 code apart from the tests, it was developed to address problems with failures of rotational symmetry in the diagnostic downscaling routines, and will be used there soon. All answers are bitwise identical, but there are two new publicly visible interfaces.
Resolves merge conflict due to new PLM_WLS test.
This commit refactors normal velocity/transport calculation in update_OBC_segment_data, using the orientation-agnostic loop ranges. Note that there is still an orientation-branch because the grid spacing G%dyCu and G%dxCv are different. * The field segment%h is removed, as we can just use h directly. * Comments are added to better illustrate the subroutine structure
Remove inaccurate warning message in logs when OBC_TEMP_SALT_NEEDED_BUG is false but temperature and salinity are needed and provided.
This commit fixes the bug that OBC segment data does not update if only value mode is used for all segments. Previously, there were two issues, 1. Subroutine update_OBC_data, which is a wrapper over subroutine update_OBC_segment_data is only called if global flag OBC%update_OBC is True. This flag, however, is only turned on if there OBC is specified by files. 2. Inside of update_OBC_data, subroutine update_OBC_segment_data is only called if global flags OBC%any_needs_IO_for_data or OBC%add_tide_constituents is True, which again prohibits OBC update with "value". Example: if normal velocity is specified by value, the transport is calculate once during initialization using initial thickness. And the transport is never updated through the run. And the model will restart with a different transport. A new runtime parameter OBC_VALUE_UPDATE_BUG is added to fix/recover this bug. In addition, OBC%update_OBC is set to True if initialize_segment_data is called. a separate procedure is added to set any_needs_IO_for_data, for clarity.
Added diag_mediator capabilities that had been developed and used in the SIS2 and MOM_ice_shelf versions of the diag mediator, in preparation for some of these separate versions to be combined into a single version. This includes adding an overload to the register_scalar_field() interface so that scalars can be registered with or without providing a null axis, retaining the existing interface while adding the one that is widely used with the ice shelf code. The previous routine register_scalar_field() is now register_scalar_field_CS(), but it has the same public interface. This commit also adds the new public interface MOM_diag_send_complete() that calls the FMS routine that flushes the IO buffers. Although this call can be too expensive for production runs, it can be very useful for debugging. Internally, register_diag_field() now checks for standard axes for 2-d fields, analogously to what was previously done for 3-d fields, thereby avoiding the allocation of a separate axis group for each 2-d diagnostic. This has been the case for ice-shelf code for some time, and it will save memory will giving identical results. The testing for incompatible mask and array sizes was consolidated into a single call each in post_data_2d_low() and post_data_3d_low(), paving the way for adding the option of masking static fields. All solutions are bitwise identical and diagnostic output (including metadata) is identical to what was previously output, but there are two new public interfaces.
Renamed the q-point index space axis names to be uniformly Iq and Jq, rather than having names that change depending on whether or not symmetric memory is being used. This commit also changes the long names of the grid space axes to be consistent between the x and y axes and to mirror those that have been adopted in the ice shelf code. The order of the logical tests for index space axes and symmetric memory around the calls to diag_axis_init were swapped to that corresponding cases are next to each other, perhaps making any inconsistencies more obvious. There were also minor changes to the code setting the axis label variables to follow the MOM6 case convention for indices and the spacing conventions for assignments, as described in the MOM6 style guide. All solutions are bitwise identical, but there are changes to the names of the axis label variables and their long descriptions when USE_INDEX_DIAGNOSTIC_AXES is true.
Corrected the cell_methods for 8 families of tracer content tendency diagnostics. Without this change, the 2-fold downscaled tendencies are 4 times too large, and the 3-fold downscaled tendencies 9 times too large. The misleading comment describing the local `conv_units` string variable in register_tracer_diagnostics() was also fixed. All solutions and standard diagnostics are bitwise identical, but there are changes in the values of downscaled diagnostics and changes in metadata.
Makedep preprocessor parsing assumed define macros of the form `#define macro(X) some-macro` but did not support empty macros of the form `#define macro(X)`. Preprocessor flag macros were assumed to be of the form `-DMACRO=VALUE`, and did not support empty macros `-DMACRO`. This patch addresses both issues.
…_vars - To reproduce across restart files requires an additional halo update on diffu and diffv in MOM_dynamics_split_RK2.F90 - To reproduce across layouts with vertex shear option requires a halo update on thickness before remapping the viscosities on the vertex points.
- Setting BBL_USE_TIDAL_BG = True causes the model to fail dimensional consistency test. If you set the spatial scale as m_to_L instead of m_to_Z in MOM_set_viscosity it fixes the issue. - Fixing documentation fo L/T scaling of tideamp - Change suggested Hallberg
#1079) This PR corrects an implementation of a finite-element solver for ice-sheet/shelf velocity. It computes a Jacobian weight J_q in shape functions for quadrature points and applies it to integrals over cells. This correction is significant when the shape of the cell is not rectangular. The same corrections are applied for subgrid parameterizations. Additionally, it fixes missing deallocations of several arrays. Changes made: - nitialize_ice_shelf_dyn: allocates a new variable `Jac` - bilinear_shape_fn_grid: computes a Jacobian determinant (`Jac`) at each quadrature point - CG_action: uses `jac_wt= Jac * IareaT` for computing the ice deformation and basal traction terms - matrix_diagonal: uses `jac_wt= Jac * IareaT` for computing diagonal terms - CG_action_subgrid_basal: computes `jac_sub_wt` using new arguments `dxCv_S`, `dxCv_N`,`dyCv_E` and `dyCv_W` that determine cell-edge spacing - ice_shelf_dyn_end: deallocates `Jac`, `CS%sx_shelf`, `CS%sy_shelf`, `CS%u_flux_bndry_val`, `CS%v_flux_bndry_val`,`CS%calv_mask`, `CS%Phi`, `CS%Phisub`, `CS%PhiC`. Copilot (model Claude Sonet 4.6) assisted with this PR.
…1089) * Vendor sphinxcontrib-autodoc_doxygen into docs/_ext/ This is piece 1 of the sphinx toolchain upgrade. The docs build previously depended on four forks which haven't been updated since 2020. The chain pinned the entire build to Sphinx 3.2.1, requiring a growing list of transitive version ceilings (jinja2<3.1, sphinxcontrib_*<1.0.x, alabaster<0.7.14, setuptools<82). This commit handles the sphinxcontrib-autodoc_doxygen dependency. The upstream (rmcgibbo) has been dormant since June 2021; the fork was effectively a rewrite (~90% of xmlutils.py changed, both documenters replaced, the whole domain switched from cpp to f). Monkey-patching the upstream package was rejected because there was nothing stable to patch against. Instead the fork's code is vendored in-tree at docs/_ext/autodoc_doxygen/ so it can be debugged, blamed, and edited like any other project source. Changes vs the fork at tag 0.7.13: - Renamed doxynamespace.rst template to doxymodule.rst and dropped the unused C++ doxyclass.rst (DoxygenClassDocumenter was already unregistered in the fork). - Removed ~80 lines of dead code: commented-out pdb.set_trace lines, the unregistered DoxygenClassDocumenter class, the dead visit_ref_angus alternative, _import_by_name_original, and the unused try-import of flint. - Dropped Python 2 compatibility cruft: from __future__ imports, from six import itervalues (replaced with dict.values()). - Ported four Sphinx 3 -> 8 API changes: * DoxygenAutosummary.get_items now passes self.bridge (the DocumenterBridge) to documenter constructors instead of self, so Documenter.__init__ finds directive.genopt. * self.warn / self.directive.warn calls (removed in Sphinx 8) replaced with sphinx.util.logging.getLogger(...).warning. * Wrapped env.doc2path in os.fspath to silence the RemovedInSphinx90Warning about str paths. docs/conf.py: prepend _ext to sys.path and rename the extension in the extensions list from sphinxcontrib.autodoc_doxygen to autodoc_doxygen. docs/conf.py also carries a monkey-patch for a parallel-build bug in upstream VACUMM/sphinx-fortran's FortranDomain.merge_domaindata. The upstream method references an undefined `outNames` (typo for `ourNames`) and also unpacks the wrong tuple shape from modules and objects dicts. With -j > 1 the merge silently fails and the f domain ends up empty, losing every Fortran cross-reference. The patch is heavily commented with a TODO to submit an upstream PR and drop the workaround during piece 2 of the upgrade. docs/requirements.txt: drop the sphinxcontrib-autodoc_doxygen git dependency. The other three forks remain pending pieces 2, 3, 4. Validation: built against a new venv.sph8 (Python 3.11 + Sphinx 8.2.3 + sphinx-rtd-theme 3.1.0 + sphinxcontrib-bibtex 2.6.5 + upstream VACUMM/sphinx-fortran master). `make html` produces 220 warnings vs 224 on the Sphinx 3 baseline. The f domain indexes 48 modules and 538 objects. Equation pages are effectively byte-identical (off only by Sphinx 8's section-id hoist optimization and the pilcrow-as-CSS theme change). Module pages are slightly larger because more cross-references now resolve than in the baseline build. * Swap sphinx-fortran fork for upstream VACUMM pinned commit Piece 2 of the sphinx toolchain upgrade. The jr3cermak/sphinx-fortran@1.2.2 fork existed to add Sphinx 3 API compatibility in 2020. Upstream VACUMM/sphinx-fortran has continued evolving since then (Sphinx logging API migration, parallel read support, modern setuptools, removal of the future dependency) and now works with Sphinx 8 directly. Upstream has not cut a PyPI release past 1.1.1 so we pin to a specific master commit for reproducibility. The pinned commit ships a broken FortranDomain.merge_domaindata that loses every f-domain object in parallel builds; the monkey-patch in conf.py from the previous commit handles that. We will revisit that patch and the pin together once an upstream PR is merged (see the TODO(piece-2) marker in conf.py setup). Validation: rebuilt docs with the pinned commit installed into venv.sph8. `make html` produces 220 warnings (identical to the pre-swap sph8 build), the f domain indexes 48 modules and 538 objects, mom.html has 36 internal cross-reference links, and f-modindex.html is generated. No behavior change from the sph8 build validated in the previous commit. * Drop sphinx fork for stock Sphinx 8 Piece 3 of the docs toolchain upgrade plan. This commit removes the last need for the patched sphinx fork and collapses the long list of transitive version ceilings it forced. Changes: - requirements.txt: drop git+https://github.com/jr3cermak/sphinx.git @v3.2.1mom6.4. Pin sphinx>=8,<9 from PyPI. Drop the eight ceiling lines (jinja2<3.1, sphinxcontrib_*<1.0.x, alabaster<0.7.14, setuptools<82.0.0) and the workaround comment around sphinx-fortran's broken requirements.txt; those existed only because Sphinx 3.2.1 dragged them in. Add lxml as an explicit dep (it is required by the vendored _ext/autodoc_doxygen extension to parse Doxygen XML). Keep six because the pinned sphinx-fortran commit imports it at module load time. - conf.py: add a monkey-patch for sphinx.util.math.wrap_displaymath that replicates the functional changes in the jr3cermak fork. The patch detects parts containing `\begin{equation}`, `\begin{eqnarray}`, or `\begin{align}` and emits them verbatim (no outer `\begin{equation}\begin{split}...\end{split} \end{equation}` wrapping). Without this, MOM6's math-heavy sources produce nested LaTeX environments that pdflatex chokes on. The patch is heavily commented with the rationale and a TODO marker for a possible upstream contribution. Verified by unit-testing the function directly (plain math still wrapped, explicit-environment math passed through verbatim) and by inspecting the generated _build/latex/ MOM6.tex: 0 double-wrapped equation/eqnarray/align blocks, 100 + 18 + 17 explicit environments preserved as top-level blocks. - conf.py: tighten exclude_patterns. Sphinx walks the entire source tree by default, which previously caused it to descend into local virtualenvs (`venv.sph3/`, `venv.sph8/`, `venv-3.11/`, etc.) and pick up LICENSE.rst, README.rst, and stray autosummary template files from inside site-packages. Each of those generated a "document isn't included in any toctree" warning and a corresponding stray .html file in the build output. The new exclusions cover venv*, venv-*, _build.*, and _ext/*/templates. Final build: 90 HTML files (down from 148, the lost 58 were all venv junk and template scaffolding) and 122 warnings (down from 220, with the same project- internal warnings as before plus the venv noise removed). - Makefile: add the equation post-processing hook to the html target, gated on UPDATEHTMLEQS=Y. The jr3cermak/sphinx fork carried this as a patch to sphinx/cmd/build.py so it ran after every sphinx-build invocation. Now invoked from the Makefile so it works against stock upstream Sphinx. The existing nortd target already had the same hook with different arguments (-p APIs -b doxygen vs -p html -b sphinx). The forked sphinx-fortran fix and the upstream PR for the math wrapping change are deferred to a follow-up. Validation: built from a fresh venv (venv.sph8.fresh) created by `python3.11 -m venv` and `pip install -r requirements.txt`. The install resolves cleanly with no manual intervention. `make html` succeeds with 122 warnings, indexes 48 modules and 538 f-domain objects, generates f-modindex.html, and produces 36 internal cross- reference links on mom.html. `make html UPDATEHTMLEQS=Y` runs the post-processing script successfully. `make latex` produces 0 double-wrapped math environments. * Drop flint dependency from docs build Piece 4 of the docs toolchain upgrade plan. With this commit the dependencies in docs/requirements.txt are only on maintained upstream sources. flint (marshallward/flint, vendored as jr3cermak/flint at 0.0.1) was historically pulled into the docs build with the intent of patching up Doxygen's incomplete parsing of Fortran functions with `result()` clauses. The fork's sphinxcontrib-autodoc_doxygen imported it at module load time but never actually called it; the import was a no-op behind a `try: import flint; except: pass`. We removed that import as part of piece 1's dead-code cleanup, so the vendored extension no longer touches flint at all. Verified empirically: uninstalled flint from venv.sph8.fresh and rebuilt. The build is byte-identical to the previous one with flint installed (122 warnings, 48 f-modules, 538 f-objects, 36 internal cross-reference links on mom.html, f-modindex.html present, 90 html files). Nothing in the current pipeline depends on flint. The two surviving "flint" mentions in docs/_ext/autodoc_doxygen/ are both in comments describing the original intent, not live code; they are left in place as historical context. If a future enhancement wants to actually fix Doxygen's result() clause parsing, that should be a separate piece of work and would not necessarily resurrect flint specifically. docs/README.md is updated to drop the flint bullet from the requirements list and the flint row from the credits table. Other stale entries in that section (sphinxcontrib_autodox-doxygen, the sphinx and sphinx-fortran fork rows) remain and will need a broader README sweep in a follow-up. * Sphinx upgrade follow-up fixes Consolidates four small fixes that followed the Sphinx 8 upgrade: - Switch all Makefile targets from `sphinx-build -b` to `-M` for consistent parallel build support. - Update docs/README.md to match the modernized toolchain (remove references to the four sphinx* forks, update requirements and credits sections). - Fix crash in autodoc_doxygen visit_image on empty <image> elements (node.text is None when doxygen produces <image/> without caption text). - Add __pycache__ and *.pyc to docs/.gitignore so the vendored extension's bytecode does not appear in git status. * RTD build infrastructure Consolidates four fixes for the Read the Docs build environment: - Override RTD's default sphinx runner with a build.jobs.build.html entry in .readthedocs.yml that runs sphinx-build -M html -j auto, since RTD's high-level sphinx: key has no parallelism option. - Make doxygen_xml path resolution robust to cwd changes: resolve relative paths against app.confdir rather than the ambient cwd, which differs between local builds and RTD. - Switch Makefile html target from -j 4 to -j auto to match RTD. - Enable html_static_path = ['_static'] in conf.py so custom CSS files (autodoxysource.css) are copied into the build output. Previously commented out, which meant app.add_css_file() added the <link> tag but the file was never deployed. * Fix O(N^2) XPath scans in autodoc_doxygen Two instances of the same class of bug in the vendored autodoc_doxygen extension, both scanning the entire merged doxygen XML tree (1230 files, 109 MB with programlisting) on every call: 1. scanNode used `node.xpath('//latexonly')` etc. In XPath, `//` starts from the document root, not from `node`. Since the extension merges all XML into one tree, every call scanned the whole tree. Fix: `//` -> `.//` (descendant-of-self). Recovery: 6m33s -> 48.8s wall clock at full input scale. 2. visit_ref used `get_doxygen_root().findall('.//*[@id=X]')` to resolve each prose <ref> by linearly scanning every element in the merged tree. With XML_PROGRAMLISTING=YES the tree tripled in size, making this the dominant cost (250s of 911s serial). Fix: lazy {id: element} dict built once on first use, O(1) lookup thereafter. Recovery: 911s -> 676s serial. * Add source code browser to Sphinx docs Generate one HTML page per Fortran source file with syntax highlighting and clickable cross-references. Clicking an identifier in the source jumps to its API documentation page; each API entry gets a [source] link back to the source listing. Implementation: - Enable XML_PROGRAMLISTING in Doxyfile_rtd (standalone commit originally, now folded in). - New autodoxysource directive in _ext/autodoc_doxygen/ that walks doxygen's <programlisting> XML, emitting per-line anchors, CSS-classed highlight spans, and pending_xref nodes for identifiers. - Stub generator produces 329 :orphan: source pages under api/generated/source/. - [source] links added to DoxygenMethodDocumenter, DoxygenTypeDocumenter, and DoxygenModuleDocumenter. - CSS styled to match the sphinx_rtd_theme's code blocks (font stack, size, colors, line gutter). - Node count optimization: coalesce text runs in _walk_highlight (10-30 nodes/line -> 3-5), set support_smartquotes=False on the source container to skip the smartquotes transform. Combined with the visit_ref fix, total build time with the source browser is 115s parallel (vs 360s before optimization). * Improve API documentation Consolidates five improvements to the rendered API pages: - Show Fortran type declarations for function/subroutine parameters and derived-type members. Types are rendered as inline code (e.g. ``real, intent(in), optional``) by looking up the <param> elements on the parent <memberdef>. - Move [source] links from the top of each entry to the bottom, so the description and parameters appear first. - Fix functions (as opposed to subroutines) not being registered in the sphinx-fortran domain. format_name() was prepending the return type (e.g. "real") which sphinx-fortran's f_sig_re regex cannot parse, leaving function entries without anchors. - Add Functions and Source Files index pages to the API Reference, with cross-reference links to module pages and source browser pages respectively. - Exclude upstream CVMix sources from doxygen input to avoid warnings from undocumented third-party code. * Fix multi-line math blocks losing indentation in LaTeX Multi-line math content (e.g. \begin{align}...\end{align}) from doxygen XML was emitted as: .. math:: \begin{align} \mathbf{v} = \mathbf{u} + ... \end{align} RST requires directive body content to be indented. The continuation lines at column 0 were parsed as regular text by docutils, producing garbled LaTeX output with "Runaway argument" and "Missing $" errors — 164 LaTeX errors on the Notation page alone. Fix: in visit_formula, detect multi-line math text and emit it as a properly indented block: .. math:: \begin{align} \mathbf{v} = \mathbf{u} + ... \end{align} Single-line math is unchanged (stays on the .. math:: line). After fix: `make latexpdf` LaTeX errors drop from 164 to 10. The remaining 10 are unrelated (9 Unicode Greek characters in source comments that pdflatex cannot render, 1 duplicate label).
Under optimization -O2 with ifort we had some floating-point divergences in the new Recon1d schemes. Two schemes diverged from their counterparts that use the old "select" pathways. One scheme failed self-consistency tests due to inlining allowing some aggressive optimization to occur. Recon1d_PPM_H4_2019: - Add explicit parentheses in end_value_h4() to enforce evaluation order and ensure bit-reproducibility across compilers. Recon1d_PPM_hybgen: - Refactor reconstruct() into a cleaner two-pass approach using new private helpers (hybgen_ppm_coefs, bound_edge_values, check_discontinuous_edge_values) copied from phased-out modules to avoid external dependencies. Bit-for-bit with the old PPM_HYBGEN scheme. - Add average() override that clamps sub-cell averages to [min(ul,ur), max(ul,ur)], reproducing force_bounds_in_subcell behaviour. Recon1d_PLM_WLS: - Fix LS_error(): remove 0.5/0.25 factors and reformulate as squared deviation from the optimal slope. - Also updated Doxygen math which had some errors. - In check_reconstruction(), copy state via assignment (perturbed = this) since an inline call to %reconstruct() on same data was yielding LSB differences.
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
17 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.
No description provided.