Added option to turn off baroclinic OBC time filtering by kshedstrom · Pull Request #995 · mom-ocean/MOM6

kshedstrom · 2019-09-06T00:49:05Z

This allows us to compute phase speed for tangential flow from tangential flow instead of using the normal flow phase speed. Not sure it helps.

Also, swapped in and out for the tracer reservoirs, to match my expectations.

codecov-io · 2019-09-06T01:05:11Z

Codecov Report

❗ No coverage uploaded for pull request base (dev/gfdl@a8dce64). Click here to learn what that means.
The diff coverage is 0%.

@@             Coverage Diff             @@
##             dev/gfdl     #995   +/-   ##
===========================================
  Coverage            ?   43.16%           
===========================================
  Files               ?      213           
  Lines               ?    62350           
  Branches            ?        0           
===========================================
  Hits                ?    26912           
  Misses              ?    35438           
  Partials            ?        0

Impacted Files	Coverage Δ
src/core/MOM_open_boundary.F90	`31.07% <0%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a8dce64...53b8be5. Read the comment docs.

MJHarrison-GFDL

Removing the time filter is identical to setting gamma to zero. So why add this option?

kshedstrom · 2019-09-06T16:43:50Z

What I'm really adding is the option to compute tangential velocity phase speed from itself rather than reusing the normal velocity phase speed. I can check for gamma==0 to trigger the code instead if you prefer.

* 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 reverts commit 401ec5b.

* 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

* +*Fix 3-equation ice-ocean flux iteration (#972) 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 https://github.com/NOAA-GFDL/MOM6/issues/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> * Fixes shelfwave failure in debug mode - rotated OBC%segment%num_fields needs to be set. * Make sure reversed segments get rotated. * Bugfix for default TIDES_ANSWER_DATE in SAL 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. * Fixes a typo in Recon1d PPM limiter 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. * Frequency-dependent drag in tensor form 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. * Adds a PLM reconstruction scheme using least squares for the slope 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 unit descriptions in 64 comments 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. * *Update TC testing parameters for late 2025 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 (#985) * 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 * Extend the PGF reconstruction to allow PLM-WLS 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 (#983) * 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. * Correction on total column thickness for wetting 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 * Allow overshoot for for grounding test In commit b8c807be327c0, 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. * Initialize and integer only set on root_PE() 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. * Fix an uninitialized float in set_viscous_ML() `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. * Add floor to "h_marg" in continuity_PPM 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. * *+Fix CHANNEL_DRAG with bathymetry above sea level 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 66 unit descriptions in comments 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 (#995) * 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 * (+) Decouple FMS infra from framework 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`. * Added latent heat flux from ice shelf to ocean fluxes * Fixes wrong number of levels in z-coord diags 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. * Fix bug in registration of ALE sponge diagnostics for generic tracers (#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 * Fix reference of (rarely) unassociated pointer 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. * Adds 5 CMIP7 diagnostics for vertically integrated heat/salt content 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 * Optimized the ice-shelf CG scheme by reducing the number of times reproducing_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. * Check that frazil is allocated before adding it to ice-shelf water flux calculation. Needed for runs without frazil. * Added melt_mask for ice shelves * Added melt_mask to ice-shelf restart * comments and units * subroutine ice_shelf_solve_inner: Completed variable descriptions and units; converted cg_halo and max_cg_halo from real to integer * +Add trim_trailing_commas and ints_to_string Copied the function i2s from MOM_diag_mediator into the function ints_to_string in MOM_string_functions, and moved the code removing trailing commas from two places in MOM_diag_mediator into the new function trim_trailing_commas in MOM_string_functions. Because of the duplication of code between MOM6, SIS2 and the MOM6 ice shelf code, these functions would need to be replicated 3 or 6 times without these changes. Also added unit tests of both new functions to string_functions_unit_tests. All answers are bitwise identical but there are two new public functions in MOM_string_functions. * Call trim_trailing_commas from register_diag_field Call trim_trailing_commas from register_diag_field and register_static_field and ints_to_string from trim_trailing_commas and eliminated the now redundant routine i2s. All code functions exactly as before but there is less duplicative code. * Refactor nsten_halo in routine advect_tracer Move nsten_halo out of iteration loop * Fix OBC indexing bug in MOM_tracer_advect Fix a bug that tracers in domain outside of the OBC is falsely updated when then the OBC is in the interior. The bug was due to an indexing error in routine advect_x. * MOM_interpolate: use get_axis_size() The prior version of `get_external_field_info` incorrectly relied on the `size` output of `get_external_field_info_infra` to determine the size of an external field's axes, since all external fields are assumed to be domain-decomposed. Since axis metadata is generally opaque, we have introduced a new infra function, `get_axis_data`, which returns the size of an axis. * ice-ocean-nolib: Fix SIS2 paths Patch to fix the SIS2 paths in the pipeline CI script. Explicitly excludes the icebergs stub, since we are using the actual icebergs model. * Correct the path to the Icepack interfaces The previous attempt to fix the automated no-library build of the ice-ocean model incorrectly specified the path to the Icepack_interfaces. This has now been corrected from `src/SIS2/config_src/external/Icepack_interfaces` to `src/SIS2/config_src/external/Icepack_interfaces` in pipeline-ci-tool.sh. The real mystery here is why the testing on the previous PR actually worked. * Delete unneeded masks args from 25 post_data calls Removed redundant mask arguments from 25 post_data() calls for 2-d arrays that were using masks that would have been set anyway based on the axes of these diagnostics. Explicit masks are only required for arrays that use unusual masks, pass atypically sized arrays (e.g., just the computational domain), or are static diagnostics that do not evolve in time. All answers and diagnostic output are bitwise identical. * Add 2D meanSL field The spatially varying time mean sea level meanSL is used as a reference height to calculate, e.g., time mean ocean column thickness max(meanSL + bathyT, 0.0). This field allows the model run in a domain with spatically varying mean height, e.g. the Great Lakes system. This first commit insulates the changes from the rest of the model. It only adds the field to ocean_grid_type and dyn_horgrid_type, the transcription between the two types, and a routine to read it from a file. The field is not yet used by the rest of the code. * Use meanSL to calcualte mean column thickness This commit uses G%meanSL in 13 modules. The change is essentially replacing G%bathyT + G%Z_ref with G%meanSL + G%bathyT. Note that this does NOT mean parameter G%Z_ref is replaced by G%meanSL. G%Z_ref is factored in both G%meanSL and G%bathyT and it is kept as a useful consistency testing tool. Another cosmetic change is made by using G%meanSL + G%bathyT, instead of G%bathyT + G%meanSL, which (hopefully) can be easily interpreted as G%meanSL - (-G%bathyT). * Modify max_depth calculation using meanSL max_depth is really used as a maximum static thickness throughout the model, so meanSL needs to be considered. * +Fix how missing values are handled in post_data At no point does MOM6 code actually set arrays passed to the post_data() to have a missing value. Instead a missing value is set in output files entirely by masking. This commit eliminates the logic that would (inaccurately) try to reset fields that seem to match rescaled missing values to the output missing value. The previous code was inaccurate, in that a rescaled field could have taken on the unscaled missing value as a valid data point and still have been incorrectly marked is missing, although the odds of this happening are exceptionally small and it would only be cases with dimensional rescaling where this could have applied. For 2-d diagnostics, this commit eliminates a duplicative array syntax math expression that did exactly what the code now does. All solutions are identical, and because the missing value was not being explicitly it is unlikely that any diagnostics will change. * Remove ice-sheet melting/freezing contribution to fluxes%latent because it is already accounted for in fluxes%sens * +Add G%IdxCu_OBCmask and G%IdyCv_OBCmask Added the new elements `IdxCu_OBCmask` and `IdyCv_OBCmask` to the `ocean_grid_type` and `dyn_horgrid_type` to facilitate the application of no-gradient open boundary conditions at faces with essentially no added overhead. These new arrays are set initially in `set_derived_metrics()` and `set_derived_dyn_horgrid()`, but may be reset in `initialize_masks()` and `open_boundary_impose_land_mask()`. All answers are bitwise identical but there are a pair of new 2-d arrays in two transparent grid types. * Use G%IdxCu_OBCmask in 7 places Modified the code to use `G%IdxCu_OBCmask` and `G%IdyCv_OBCmask` in 7 places each in 6 modules. They are used instead of `G%OBCmaskCu*G%IdxCu` and `G%OBCmaskCv*G%IdyCv`, to which they are equivalent. This change should slightly speed up the model, and as expected all answers are bitwise identical. * Add option to scale tidal amplitude for bottom ustar. (#1016) * Add option to scale tidal amplitude for bottom ustar. - previously we used the tidal amplitude to compute ustar. - The additional factor translates between amplitude and time mean tidal current. - Setting the factor TIDEAMP_FACTOR<0 preserves old answers. * Update tideamp factor implementation for efficiency - factor out the negative "default" value to automatically set to multiply by 1.0 instead of using an if-block. - factor in the c-grid averaging 0.5 to further reduce extra operations, but clearly label the parameter to reflect this. --------- Co-authored-by: brandon.reichl <brandon.reichl@noaa.gov> * Add vertical tracer flux diagnostic for dye tracers (#1022) * Add vertical tracer flux diagnostic for dye tracers - Register vertical flux diagnostic in initialize_dye_tracer - Calculate net vertical flux from entrainment (positive upward) - Post flux diagnostic in dye_tracer_column_physics * changed diagnostic registration to be at interface, made sure boundary fluxes are zero * changed lines 338 and 354 as needed. Fixed accidental space on Line 1. * Regroup MOM_initialize_fixed params in param_doc This commit is meant to fix the issue that all parameters in MOM_initialize_fixed after OBC are logged under module MOM_open_boundary in MOM_parameter_doc. By moving log_version call after OBC, parameters from MOM_initialize_fixed are now logged under three "modules" in MOM_parameter_doc: 1. Parameters before OBC are under module MOM_grid_init, which also (incorrectly) includes topography relatd parameters. 2. module MOM_open_boundary 3. Parameters after OBC are under module MOM_initialize_fixed. The change makes sure OBC parameters are well separated from the other parameters. This is a hack rather than a fix. * Minor open_boundary_config refactor * Make OBC related calls in MOM_initialize_fixed explicitly conditional for readibility. * Early return in open_boundary_config if there is no segment, so that OBC is not allocated and then deallocated. * Refactor gradKE with simpler treatment of OBCs Refactored gradKE() to only take a single layer of velocities and thicknesses and to use G%IdxCu_OBCmask and G%IdyCv_OBCmask to avoid extra loops over the OBC segments for efficiency. With these changes, there is also an expectation that gradKE() will be more likely to be inlined or threaded and hence become more efficient. All solutions are bitwise identical, but there are some diagnostics of accelerations that are fully or partially set to zero at solid (masked) boundaries. * Autoconf: Build dependencies as libraries This patch significantly modifies the Autoconf-based builds to separate external content in `pkg/` from the main MOM6 source code. These changes are meant to address the challenges with analysis of source code which is outside of our control. Such content is now built separately, isolating it from the usual MOM6 code requirements. A secondary benefit is that a submodule checkout is no longer required for ocean-only builds. The patch includes the following specific changes. * Builds relying on `ac/deps/` now build libraries for GibbsSeaWater (`libgsw.a`) and CVMix (`libcvmix.a`) alongside `libFMS.a`. This replaces any content placed externally into `pkg/`. * A new macro, `MOM6_FC_CHECK_LIB`, improves library detection tests by including support for both subroutines and functions, as well as generic argument lists. * Makefiles across the build system (`ac/deps/`, `.testing/`, etc.) have been updated to generate and use the new libraries. These changes are not intended to break existing builds, as summarized below. * `src/parameters/CVMix/` and `TEOS10/` are explicitly excluded from the makedep source trees, but the symbolic links to `pkg/` are retained, and existing builds which assume their presence should still work. * The legacy `AX_FC_CHECK_LIB` macro has been retained, so that it remains available to any external MOM6-examples builds. * The git submodules have not been removed, although the GitHub Actions no longer uses then in testing. The patch also includes the following related minor modifications. * The GitHub Actions FMS build stages now includes GSW and CVMix * `.gitignore` files are modified to target files produced by the build system, rather than globbing of various incidental files. * A non-POSIX shell operation `==` in `configure.ac` has been replaced. * ac: Configure use of pkg/ source `ac/deps/Makefile` and `.testing/Makefile` now provide optional configuration of a target `pkg/` directory, using the PKG macro. To build using the submodule source code, PKG is unset:: make -j In this case, symbolic links to `pkg/` are used to build CVMix and GSW libraries. Modifications to the equivalent files in `src/` will be seen by the library builds. Developers who prefer to edit symbolic links from rogue submodules can still work as before. To build using `libcvmix.a` and `libgsw.a`, unset `PKG`:: PKG= make -j or make -j PKG= To accommodate this, there are now explicit macros `MAKEDEP` to specify the path of `makedep`. This prevents some of the "pathing up" issues that should never have been used, and were exposed by the new symlinks in to `pkg/`. In just about every normal use case, these are pre-configured. Macros for Python and autoconf detection were also modified to check preset `PYTHON` and `MAKEDEP` settings before attempting to search `PATH`. Following autoconf convention, the preset values of `PYTHON` and `MAKEDEP` are not reported. * Deleted about 170 unused internal variables Eliminated about 170 unused internal variables scattered across 26 files. Most of these were deleted outright, but a few that are associated with commented out blocks of code were just commented out. All answers are bitwise identical. * Decrease open_boundary_config indentation This is a follow up to commit 27f4243 (break down for easy review). if (OBC%number_of_segments > 0) is removed from open_boundary_config, so that the indentation level can be reduced for most part of the routine. * Correct the unit of diagnostic "eta_cor" It should be in thickness unit (m for Boussinesq and kg/m2 for non-Boussinesq). The long name description is also slightly changed. * Merge [uv]bt_st and [uv]bt_first in btstep The two pairs are used for different purposes. `[uv]bt_st`` are used for diagnosing barotropic momentum tendency and `[uv]bt_first` are used for OBC. But they are equivalent (initial BT velocity) and neither uses wide halo. This commit merges the two. Extra halo points in assigning `[uv]bt_st` are also removed. * Fixes for flang * Change variable not type * Bugfix for OBC_SPECIFIED_STRAIN option Specified shear strain is not properly applied at OBCs due to a faulty if-block in MOM_hor_visc. A runtime flag OBC_SPECIFIED_STRAIN_BUG is added to address this bug. * Add OBC_VORTICITY_CONFIG and OBC_STRAIN_CONFIG Use string two runtime parameters OBC_VORTICITY_CONFIG and OBC_STRAIN_CONFIG to replace two sets of four mutually exclusive logical runtime parameter for vorticity and strain. The two new config parameters include five options: zero, freeslip, computed, specified and none. OBC_VORTICITY_CONFIG replaces * OBC_ZERO_VORTICITY * OBC_FREESLIP_VORTICITY * OBC_COMPUTED_VORTICITY * OBC_SPECIFIED_VORTICITY OBC_STRAIN_CONFIG replaces * OBC_ZERO_STRAIN * OBC_FREESLIP_STRAIN * OBC_COMPUTED_STRAIN * OBC_SPECIFIED_STRAIN The old parameters are still allowed and a warning is given when they are specified in MOM_input. However, MOM_parameter_doc only logs the two new config parameters. * Add clocks timing MOM_end and MOM_save_restart Added the new clocks 'Ocean MOM_end' and '(Ocean MOM save_restart)' to record the time spent during the MOM part of the end of run cleanup and the time spent saving restarts. These are both included in the time under the 'Ocean' clock, which now also includes the time spend in 'Ocean Initialization'. The addition of these clocks required that the calls to `MOM_infra_end()` be moved after the call to `MOM_end()` in the solo_driver and mct_driver code, as was already the case for the FMS and nuopc drivers. All answers are bitwise identical, but there is additional timing information written at the end of the run for many cases. * Add missing spaces in 53 strings Added missing spaces to strings that are concatenated across line breaks in 53 places, and rearranged strings in another 14 places that could appear to a grep as though they likely are missing spaces. Many of the problematic instances were in error messages that are probably never called, but 35 occurred in the descriptions in MOM_parameter_doc files. A total of 32 files were modified, and all solutions are bitwise identical. * Fix 9 misspellings in MOM_parameter_doc files Corrected 9 spelling errors scattered across 39 files. These spelling errors were detected in a review of the MOM_parameter_doc.all files for the OM4_025 test case. The specific errors that were corrected include "periodcially", "BOUSSINSEQ", "fn.", "the the", "diffusivites", "instabilies", "horizonal" and "amd". The syntax of one doxygen comment was also corrected. All answers are bitwise identical, but there are multiple changes in MOM_parameter_doc files. * ac: Conditional testing of GSW and CVMix libs The autoconf tests for existence of GSW and CVMix are now conditional, and are applied if `--with-gsw` and `--with-cvmix` are provided to `configure`. The `srcdir` `ac/` gymnastics have also been removed, since we no longer try to operate from the `ac/` directory and now generally copy `configure.ac` into a build directory before running `configure`. (Neither is ideal, but copying has generally proven more useful.) * Add runtime OBC_TEMP_SALT_NEEDED_BUG flag If true, this new runtime flag stops the model if no temperature and salinity data is provide for any open boundary segment, when the model is run with T/S. This is done by using segment attributes t_value_needed and s_value_needed, which were always set to False. If the bug is on, the external field array %t are always initialized to interior cell T/S by fill_temp_salt_segments and never reset to external value [there is none] in update_OBC_segment_data. A major consequence is the update to reservoir array %t_res are inconsistent across restarts. Default is True (bug is on) for backward compatibility. * Add FATALs in OBC parse_segment_manifest_str Minor refactor of parse_segment_manifest_str to check if 1. total fields > max (unlikely) 2. there are duplicate fields * Surface state diagnostic index bugfix This fixes an error in the index bounds for v-velocity in the uppermost mixed layer. * Update MOM_wave_interface.F90 Fixed the index error and initialization issues for Robust_PGF issue * Add function time_minus_signed for time differences * Use time_minus_signed for time relative to a reference (is < 0 when time < reference) * Add run-time flag USE_ADJUST_INTERFACE_MOTION USE_ADJUST_INTERFACE_MOTION=True (default) now allows us to bypass the s/r adjust_interface_motion() that is otherwise called after building all ALE coordinates. This s/r runs into round-off problems under grounded ice sheets; the total thickness of the ocean might be NK * Angstrom but the math used to make adjustments fails if the surface is depressed a long way from 0. Rather than fix the math, the newer grid generating functions are written to not need "fixing" after the fact. * ALE: Add REGRID_USE_DEPTH_BASED_TIME_FILTER parameter Add a new runtime parameter REGRID_USE_DEPTH_BASED_TIME_FILTER that controls whether depth-based time filtering is always applied or only when old_grid_weight > 0. Previously, filtered_grid_motion() was called unconditionally in all coordinate builders (z*, sigma, rho, HyCOM1, adaptive), even when no filtering was needed, incurring unnecessary computation. The new logical flag (default=.true. for back-compatibility) is exposed via set_regrid_params() and read in ALE_init(). When false, the call to filtered_grid_motion() is skipped if CS%old_grid_weight == 0, avoiding redundant work. The flag is also initialized in initialize_regridding() for non-main-parameter code paths (e.g., diagnostic and ODA remapping). Note: With CS%old_grid_weight=0 and either REGRID_USE_DEPTH_BASED_TIME_FILTER set to True or False, the results are not bitwise-reproducibility even though mathematically they should be identical. This flag is thus needed to bypass the calculations. * ALE: Fix copy-paste error in REGRID_FILTER_DEEP_DEPTH description The parameter description incorrectly referenced REGRID_FILTER_SHALLOW_DEPTH twice instead of REGRID_FILTER_SHALLOW_DEPTH and REGRID_FILTER_DEEP_DEPTH. Fixing this in current PR (separate commit) since the doc files from this module are being changed anyway. Two stones, one bird... * Separate frazil for ice sheets and sea ice (#1027) * Separate frazil for ice sheets and sea ice This commit fixes an issue where sea ice was forming under ice shelves. The fix involves removing the frazil that forms under ice shelves from the ocean surface state before is it passes to SIS2. The approach here is to extract this sub-shelf frazil from the surface state into a separate variable that is used in the ice shelf module to add frazil to the shelf. This also included adding an ice-shelf frazil restart and diagnostic. This approach is compatible with all coupler/MOM6 options regarding timestepping and ordering of model components. It is also compatible with partially-filled ice-sheet cells (see comments in new subroutine MOM_ice_shelf.F90/adjust_ice_sheet_frazil), as long as the ice front does not move over time (as is currently the case for all ice-sheet/ocean fluxes). * Add adjust_ice_sheet_frazil to init routines * Remove frazil restart from ice shelf * Add missing use statements for adjust_ice_sheet_frazil * Enable exchange of runoff DOC from Land to Ice used in CMIP7 ESM4.5 (#1034) * Enable exchange of runoff DOC from Land to Ice - Land component calcutes the DOC (Dissolved Organic Carbon) content flux in runoff (KgC/m2/s) - COBALT needs access to that runoff flux (currently done via data_override from file) - Ocean component needs to make the flux available to COBALT - Ice component needs to make the flux available to ICE - land_ice_flux_exchange needs to exchange the flux from Land to Ice * Enable exchange of runoff DOC from Land to Ice used in CMIP7 ESM4.5 - Apply changes per review for PR #1034 * Apply review requets and fix MOM6 CI errors * Address recent review to make carbon exchange configurable - Nothing will be allocated/avaliable unless a parameter is set in overrides # Allow Carbon Runoff diagnostics #override ALLOW_CARBON_FLUX_EXCHANGE = True * +Extend diag_mediator to allow the piecemeal posting of diagnostics (#809) * Create dynamic buffer types for diagnostics posted piecemeal Some quantities in MOM6 are calculated in subroutines which expect slices of the model's 2d or 3d arrays. Diagnosing these quantities can be challenging because the usual post_data calls expect whole arrays. To solve this problem, the changes here introduce a dynamic buffer that can grow as needed over the course of an simulation. This buffer keeps track of the diagnostic IDs and the index in the buffer. When a slot in the buffer is no longer needed, for example if the whole array has been computed and posted, the slot is marked as "available" for overwriting. The buffer is only allowed to grow if all the currently allocated slots are in use. Any computational cost associated with growing this buffer will only happen in the first few steps of the model as post_data for requested diagnostics is called for the first time in that run. * Implement piecemeal posting of diagnostics The diag mediator has been extended to add a dynamic buffer to each axes group. Three new methods have also been added to enable the piecemeal posting (by column, by point) of a diagnostic and a 'final' method to allow the buffer to be reused later. * Demonstrate piecemeal diagnostic posting in ePBL ePBL calculates the vertical diffusivity column by column. This provides a convenient sanity check of the new piecemeal posting of diagnostics. The original diagnostic Kd_ePBL is done by posting the full 3d prognostic array, whereas a new diagnostic Kd_ePBL_col_by_col posts the same array from within ePBL but does so column by column. * Add docstrings for diag_buffer procedures Fixes failures in the CI due to some procedures and type members not having docstrings. * Conform to style for if statements * Specify that buffer units are arbitrary * Allow setting of the fill value when expanding the buffer * Refactor buffers based on Marshall's suggestions * Actually set the fill value and add tests * Fix issue with partial construction of buffer * Fixes logic error when using piecemeal buffers Two issues were leading to sporadic compilation errors and segmentation faults 1) The post_diagnostics_by_point only looked up the slot in the buffer but would not expand the buffer -> now uses check_capacity_by_id which will expand the buffer if needed 2) post_data_3d_final now conditionally checks to make sure that the buffer slot is not zero before taking any action * Fix missing check on non-negative diagnostic ID Fixes a bug where the Kd_ePBL_col_by_col diagnostic's post_data calls were not being checked before being called. Also, move setting of extents for buffers on native axes to mask generation * Initialize extents after set_axes_info Update the extents for the piecemeal buffers after the axes info is set. Resolves a problem where the extents were not being retained. * Refactor update_OBC_segment_data loop indices This commit changes the indexing variables in open boundary routine update_OBC_segment_data, so that they are consistent with the rest of the model. The changes improve readability and prepares moving part of the subroutine to initialize_segment_data. The changes includes, * Use [ij][se]d and [IJ][se]dB aliases for indexings in segment%HI * Replace [ij][se]_obc with [IJ][se]dB * For interior cells, retire [ij]shift and use [ij]sd that are already saved in segment%HI Other minor changes: * Add a comment that segment%Cg is never used by the model * Change normal_trans_bt in update_OBC_segment_data from an allocatable array to a scalar * Move alloc buffer_dst to initialize_segment_data This commit moves allocating buffer_dst array of segment%field type to routine initialize_segment_data from update_OBC_segment_data. This array is always allocated and exits in the duration of the run. It makes more sense to have the allocation occur in the init subroutine, than a routine called every time step. Other changes * Simplify remapping_core_h call in update_OBC_segment_data. * make dz_stack a fix-sized array than an allocatable * Add subroutine allocate_segment_field_data This commit moves out all allocations from initialize_segment_data to a new subroutine allocate_segment_field_data, which includes buffer_src and buffer_dst. allocate_segment_field_data always assigns other attributes of segment%field, including "on_face", "scale", "use_IO" Minor updates * Remove OBC%thickness_segment_data_exists, which is never used. * Change string "genre" in segment%field to logical "bgc_tracer" * Remove GV input (unused) from subroutine scale_factor_from_name * Rename local variables in initialize_segment_data for clarity * fields -> field_names * fieldname -> varname * +Add & use optional maskT arg to initialize_masks Added the new optional argument maskT to initialize_masks() to allow for the tracer point masks to be specified separately from consideration of the bathymetry. Initialize_masks() was also refactored to set the other masks based on the values of G%mask2dT rather than referring back to the bathymetry. The second part of this commit uses this new argument to set the ice shelf masks for diagnostics in new calls to initialize_masks() from initialize_ice_shelf(). All answers are bitwise identical, but there is a new optional argument to a publicly visible interface. * +Fix bugs in ice shelf top slope calculation This commit fixes bugs with the calculation and diagnostics of the ice shelf top slope. The pending updates to MOM_ice_shelf_diag_mediator give a fatal error if the array size sent to a post_data() call is inconsistent with the mask assigned based on the declared grid location. The 'sh_shelf' and 'sy_shelf' diagnostics were failing this test because their locations were declare inconsistently to be at corner points instead of tracer points. A similar problem was present for 'surf_slope_mag_shelf', but the array used to calculate this diagnostic was also allocated with the wrong array sizes as though it were at corner points. All of these problems have been corrected, and the new version of the ice shelf diag mediator is no longer giving fatal errors related to these diagnostics. While investigating where the shelf-top slopes were actually staggered, it became clear that there were a number of directionally inconsistent bugs in the estimates of the grid spacings used to calculate these slopes. These bugs would not have been evident for a simple Cartesian grid. In the revised code, these bugs are corrected when the new runtime parameter ICE_SHELF_TOP_SLOPE_BUG is set to false. The previous answers are retained via code that has been refactored to make the bugs more obvious. In addition, the extent of several do loops in calc_shelf_driving_stress were reduced to their appropriate extents. By default (unless ENABLE_BUGS_BY_DEFAULT is false), solutions are bitwise identical, but there is a new runtime parameter in the MOM_parameter_doc files and the reported locations of some diagnostics are corrected. * +Enhance ice_shelf_diag_mediator like framework This commit greatly enhances the capabilities of the ice_shelf_diag_mediator by pulling in multiple capabilities that had been added to the standard MOM6 diag_mediator in the 13 years since the ice_shelf_diag_mediator was spawned from it (via the SIS2 diag mediator). The new capabilities include the ability to handle CMOR name options for diagnostics and related metadata, write out scalar diagnostics, record cell methods as metadata, use integer index axes in place of nominal latitudes and longitudes, print out checksums of all primary diagnostics and detect any obsolete diagnostics in diag tables. This commit is a large step toward being able to use the standard MOM6 diag mediator with the ice-shelf code and eliminate the specialized diag mediator for the ice shelves. The new public interfaces with this commit include define_axes_group(), diag_masks_set(), diag_register_area_ids(), MOM_IS_diag_send_complete() and found_in_diagtable(). None of these are exercised as yet by the ice shelf code, but new options that exercise them will appear in future commits. Also the publicly visible but opaque axesType has been renamed to axes_grp, mirroring MOM6. This commit adds the new runtime parameters USE_INDEX_DIAGNOSTIC_AXES, DIAG_MISVAL, DIAG_AS_CHKSUM and CHKSUM_DIAG_FILE to support the new capabilities. The previously exercised interfaces are unchanged, although register_MOM_IS_scalar_field() can be called either with an `axes` argument as is done in the ice shelf code or directly with a `diag_ctrl` argument as in the standard version of MOM6_diag_mediator. For now, all diagnostic and available_diags files written by the ice shelf are identical to what the previous version wrote out, but there is commented out code that will add new metadata and formatting like that generated by the standard MOM6 diag_mediator. It is the nature of this commit that the changes swept up by it are extensive, but with such extensive refactoring this is probably unavoidable and it seems like there would be minimal purpose served by breaking this up into intermediate waypoints. * +Extract find_dz_for_eta from find_eta Created the new publicly visible subroutine find_dz_for_eta() that does the calculation of the integrals that convert the layer thicknesses into height changes across layers and used this in find_eta_2d() and find_eta_3d(). As a part of these changes, find_eta_2d() and find_eta_3d() were refactored to call find_dz_for_eta() to reduce code duplication and to simplify the openMP directives. All answers are bitwise identical, but there is a new publicly visible interface. * +Add ACCURATE_NONBOUS_THICK_CELLO Add the new runtime option ACCURATE_NONBOUS_THICK_CELLO to improve the calculation of the diagnosed layer thicknesses in non-Boussinesq mode by using the same code in find_dz_for_eta() as is use for the free surface. The code to calculate the diagnostics 'e' and 'e_D' were also refactored to use the same thicknesses from find_dz_for_eta(), thereby avoiding extra calls to the equation of state. Depending on which diagnostics are used, this change may speed up the code in some non-Boussinesq cases. By default, all answers and diagnostics are bitwise identical, but there is a new runtime option in the MOM_parameter_doc files for non-Boussinesq runs. * Bugfix for reintegrating diagnostics in non-boussinesq (#1060) * Store dz_begin for reintegrating extensive diagnostics in non Boussinesq mode - We need dz_diag from the source grid to be consistent with dz_diag that was used when constructing the target grid for reintegrating extensive quantities to the diagnostic coordinates. If dz_diag is not consistent with the target grid, there can be a mismatch in total thickness of source and target grids causing the reintegrate_column routine to fail and give errors in diagnostics. - Previously we stored h_diag and recomputed dz_diag from h_diag when it was needed for non-Boussinesq models. However, if the model's T&S evolved between recomputing dz_diag, the algorithm did not reproduce the dz_diag that was used to set up the target grid. Thus the thicknesses passed to reintegrate_column on the source and target grids were inconsistent. - This commit stores dz_begin in addition to h_begin when in non-Boussinesq mode and using Z based diagnostic coordinates. It then uses the stored dz_begin for reintegrate_column, ensuring the source and target grids are consistent. - The new array is only allocated and filled if it will be used. * Typo fix on dz_diag_needed --------- Co-authored-by: brandon.reichl <brandon.reichl@noaa.gov> * Remove stray leading $ * Add setuptools<82 to keep pkg_resources available * Refactor of limit vel function for nicer layout for GPU compute (#1023) This PR modifies vertvisc_limit_vel by refactoring the single CFL test and velocity truncation jki loop into two kji loops for each stage. Much like before, the first loop tests for CFL violations and logs the smallest velocity magnitude for each ij column. The second loop independently checks for and applies potential potential truncations in the compute domain. The third write_[uv]_accel loop is preserved, with a new global control flag. The pattern is repeated for u and v. This split was required to eliminate some performance bottlenecks observed in the GPU vert_friction implementation. Some of the specific steps include: * port refactored limit vel * Revisions to vertvisc_limit_vel PR to dev/gfdl * `trunc_any_array` is promoted to edge-domain in both directions `(SZIB_(G), SZJB_(G), ...)` * Some i/j indexing errors (e.g. `h(i+1,j)` -> `h(i,j+1)`) (Mostly i- code in the j- block) * Removed `!$omp parallel do` directives (NOTE: maybe we undo this one?) * Style fixes * Index case fixes (e.g. `I` -> `i`) * `end do` -> `enddo`, `end if` -> `endif` * Operator/keyword spacings * Trailing whitespace * Removed unused CFL_based_trunc from vertvisc_CS * Update thickness check conditions Replaced the older thickness H_report conditions with the newer omes based on compute domain and using CS%h_[uv]. Also moved the write_u_accel block to mirror dev/gfdl. No idea if there are performance implications but we can check it out. * loop ordering * any to avoid loop * do any write instead of any to avoid copy of array unless necessary * recompute instead of using an array, should be cheaper * address review comments: - truncation moved inside the main big if since they are sequential, limit the scope of the write any variable - fix indentation * trunc any initialized to false with do_any_write * trunc_any set to true if CFL>CS%CFL_trunc in a separate if * Fix variable case in vertical friction calculation * Refactor velocity reporting conditions in MOM_vert_friction --------- Co-authored-by: Marshall Ward <marshall.ward@noaa.gov> Co-authored-by: Marshall Ward <marshall.ward@gmail.com> * +Add time_to_real and scaling in real_to_time Added the new function time_to_real() that wraps time_type_to_real() with the option to rescale the output time into different units than seconds, and added optional unscale and scale arguments to real_to_time() and time_minus_signed(), respectively. Identical changes were made to both the FMS1 and FMS2 versions of MOM_time_manager.F90. All answers are bitwise identical, but there is a new publicly visible interface and new optional arguments to two other publicly visible interfaces. * Remove 5 unused module uses for time_type_to_real Removed unused module use statements for time_type_to_real in 5 modules. All answers are bitwise identical. * Rescale via time_to_real or real_to_time calls Refactored 34 calls to real_to_time() to use the new unscale argument and replaced 16 calls to time_type_to_real() with calls to the new routine time_to_real(), of which 15 use a scale argument. Also use the new scale argument in 5 calls to time_minus_signed(). In four instances, extensive expressions for the end of a diagnostic interval are being stored in new temporary time-type variables. Although these changes lead to some lines that are longer than before, the use of the scale and unscale arguments makes it much easier to detect incorrect rescaling factors. All answers are bitwise identical. * Restructure segment%field of OBC This commit reworks segment%field array. Previously, the array is dynamically decided by available input fields, i.e. segment%num_fields = number of input fields + number of BGC tracers. Problems: 1) indexing of the field is not fixed. 2) There is no direct one-on-one mapping between field array and a dozen of xxx_needed flags in segment. The new approach always creates a field array with its first 13 members being the 13 physical fields, which predetermines the indexing. It also merges the xxx_needed flags to their corresponding fields. Note that only used/available fields have their buffer_src and buffer_dst arrays allocated, so there is no additional memory cost. With the new structure, initialize_segment_data is refactored to have four steps, 1) Initialize 13 physical fields, including toggling required fields using subroutine segment_determine_required_fields. 2) Create a list of available inputs from parameter file and map their indices with the 13 physical fields. 3) Allocate available physical fields. (allocated() = is_available) 4) Determine and allocate BGC fields. Misc * Two new subroutine/function * subroutine segment_determine_required_fields * function find_phys_field_index * The model stops if a variable not in the 13 physical field is given. * The model stops if there is a required input missing and informs which segment and field. * A warning is now issued if a field is provided but not needed. * MAX_OBC_FIELDS (=100) is replaced by NUM_PHYS_FIELDS (=13) * Remove 19 flags and 6 indices in segment type * Remove all x_values_needed flag and 6 indices to locate tidal fields * Remove subroutine rotate_OBC_segment_values_needed * Apply fixed indices to tidal fields in update_OBC_segment_data * Add missing doc comments to undocumented variables in ice_shelf - MOM_ice_shelf_dynamics.F90: split loop_bounds_type members (ish, ieh, - MOM_ice_shelf.F90: add missing = -1 initializer to id_dhdt_shelf A disclosure: Copilot has been used to make these changes. * Adds t17d and t20d diagnostics 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. * +Revise MOM_IceShelf available diags file format 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. * +Change some MOM_IS axis units to CMIP standards 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)`. * Diagnostic output time scaling fix in HA_accum 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. * Use clearer conversion factor for ale_u2 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. * Diag_mediator cleanup 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. * Zero out allocated array in MOM_temp_salt_initialize_from_Z 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. * Refactor update_OBC_segment_data 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 (#1072) * 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. * Change "end if" to endif and "end do" to enddo 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. * Trailing semicolons and semicolon white space 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. * Split reading files from update_OBC_segment_data 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. * Split reservoir init from update_OBC_segment_data 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. * Move particles_run call into step_MOM_tracer_dyn 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. * +Add symmetric_sum functions 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. * Refactor update_OBC_segment_data This commit refactors normal velocity/transport calculation in update_OBC_segment_data, using the …

* Add SHALLOW_ALE_RESOLUTION SHALLOW_ALE_RESOLUTION implements a HYBGEN-style Z-sigma-Z near surface fixed coordinate for HYCOM1. For example the US Navy's GOFS 3.1 HYCOM setup has 41 layers, with the top 14 layers in a Z-sigma-Z configuration. For MOM6 HYCOM1 this is: SHALLOW_ALE_RESOLUTION = 14*1.0,27*0.0 for 14 1m "shallow" layers. Let N_SIGMA be the number of consecutive non-zero entries, typically < NK. When rest depth is shallower than SUM(SHALLOW_ALE_RESOLUTION(1:N_SIGMA)) use SHALLOW_ALE_RESOLUTION. When rest depth is deeper than SUM(SHALLOW_ALE_RESOLUTION(1:N_SIGMA)) use ALE_RESOLUTION. Otherwise use a linear sum of the two weighted by rest depth. The default of all zeros turns this option off, and when off answers are unchanged. The new parmeter SHALLOW_ALE_RESOLUTION is only present when using HYCOM1. * Non-integer HYBRID_MAP values The 2-d REAL map array in HYBRID_MAP usually contains integer values each referencing one profile. It can instead contain non-integer values of the form I+frac, which indicate a weighted sum of profiles: (1-frac) p(I) + (frac) p(I+1). The same profile can be used multiple times, e.g. if 1st profile is also 4th can get profiles between 1 and 2 and between 1 and 3. HYBRID_3D is more general, but HYBRID_MAP covers most practical uses. * indent continuations, source code <= 100 chars * +Add RESCALE_STRONG_DRAG Added the new runtime option RESCALE_STRONG_DRAG, that can be set to true to reduce the barotropic contribution to the layer accelerations to account for the difference between the forces that can be counteracted by the stronger drag with BT_STRONG_DRAG and the average of the layer viscous remnants after a baroclinic timestep. In testing, this new capability eliminates some of the growing instabilities that can occur with an ice shelf and BT_STRONG_DRAG set to true. This commit also adds new diagnostics of the barotropic step viscous remnants and the eta anomalies contributing to barotropic pressure forces, either averaged over the barotropic step or at each barotropic step. By default all answers are bitwise identical, but there is a new runtime parameter and 4 new diagnostics. * Add option to horizontally homogenize the Stokes drift when used via … (#967) * Add option to horizontally homogenize the Stokes drift when used via the dataoverride surfbands procedure. * Add variable description in new method for horizontally averaging Stokes drift. --------- Co-authored-by: brandon.reichl <brandon.reichl@noaa.gov> * Fix calculation of CAv_Stokes diagnostic Corrected a horizontal indexing bug in the calculation of the CAv_Stokes diagnostic, making it rotationally consistent and consistent with the calculation of CAu_Stokes. This bug has been there since the CAv_Stokes diagnostic was originally added. The loop range over which qS is calculated was also reduced to the range over which it is used. All solutions are bitwise identical, but this commit does change the values of a (perhaps infrequently used) diagnostic. * makedep: Update interpreter directive to python3 The interpreter directive ("shebang") of makedep is updated to `python3`, rather than the version-agnostic `python`. Although we never invoke the shebang of the script, there are OS environments out there which will object to any presence of a versionless python. PEP 394 also strongly recommends the adoption of python3 as the executable name, regardless of Py2 support. * Continuity ppm port to gpu (#29) * present_vhbt_or_set_bt_cont: merge couple of loops * meridional_flux_thickness: cpu optimize a bit * meridional_flux_adjust: back to jki * set_merid_BT_cont: pull out meridional_flux_adjst * set_meridional_BT_cont: optimize cpu * rm remaining merid/zonal_flux_layer * zonal_flux_layere: improve naming a bit * zonal_flux_layere_OBC: improve naming a bit * improve flux_elem line wraps * optimize data transfers a bit * pass elem not arr * add some missing documentation * fix trailing spaces and missing var docs * add last param docs * meridional_flux_adjust: fix fpe err * fix another fpe * cleanup args for new helper subroutines * clean up enter/exit data * fix passing h twice * meridional_mass_flux: do concurrent * zonal_flux_adjust: use a few 3d tmp arrays to mirror meridional_flux_adjust * move target update out of continuity * initialize pbv%por_face_area[U/V] on GPU * cleanup some transfers from continuity_PPM * clean up a few minor things * zonal/meridional_flux_adjust: use scalar u/v_new * declare vp,up,h_tmp on gpu * remove h update * continuity_PPM: minimize mapping stmts * zonal_flux_adjust: minimize mapping stmts * set_zonal_bt_cont: minimize mapping stmts * merional_flux_adjust: minimize mapping stmts * set_merid_bt_cont: minimize mapping stmts * zonal/meridional_flux_adjust: tmp vars duhdu/dvhdv 3d -> scalar * separate alloc of private variables for gcc * u/vh_aux: 3d->2d * target teams loop recognizable by amdflang * Continuity CS outside of init This moves the Continuity CS to the dycore init function. For some reason, this avoids an answer change with CPU. (Possibly because alloc inside of a function doesn't quite match the CS outside of it?) A few minor data transfers are also added to fix up differences in the chksum log output. * zonal_mass_flux: isolate zonal_flux_layer * zonal_mass_flux: seperate local_specified_BC block loop * zonal_mass_flux: add j dim to tmp vars * zonal_mass_flux: seperate visc_rem_max init loops * zonal_mass_flux: seperate du_min/max_CFL init loops * zonal_mass_flux: separate duhdu/uh_tot_0 init loop * zonal_mass_flux: separate du_min/max_CFL aggress_adjust update case (untested) * zonal_mass_flux: separate du_min/max_CFL non-aggress_adjust update case * zonal_mass_flux: separate du_min/max_CFL non-use_visc_rem update case (untested) * zonal_mass_flux: separate du_min/max_CFL 0-clamp loop * zonal_mass_flux: separate do_I local_specified_BC init loop (untested) * copy zonal_flux_adjust that accepts 2d args * zonal_mass_flux move j loop into zonal_flux_adjust copy * zonal_flux_adjust_fused: use 2d internal arr * zonal_flux_adjusted_fused: separate all loops * zonal_mass_flux: separate u/du_cor update (untested) * zonal_mass_flux: replicate former present(uhbt) control flow * copy set_zonal_BT_cont with 2d args * zonal_mass_flux: move j-loop into set_zonal_BT_cont_fused * set_zonal_BT_cont_fused: use zonal_flux_adjust_fused * set_zonal_BT_cont_fused: separate init loop * set_zonal_BT_cont_fused: separate duR/L init loop * set_zonal_BT_cont_fused: separate u_0/L/R init loop * set_zonal_BT_cont_fused: use zonal_flux_layer_fused * set_zonal_BT_cont_fused: separate last 2 loops * copy merid_flux_layer that accepts entire arrs * merid_flux_layer_fused: separate loops * merididional_mass_flux: separate local_specified_BC loop (untested) * merididional_mass_flux: separate tmp variable init loops * merididional_mass_flux: separate dv_min/max_CFL calc loops * merididional_mass_flux: separate dv_min/max_CFL 0 clamp loop * merididional_mass_flux: separate simple_OBC_pt init loop (untested) * copy meridional_flux_adjust that accepts entire arrs * meridional_mass_flux: separate meridional_flux_adjust_fused * meridional_mass_flux: move j-loop into meridional_flux_adjust_fused * meridional_flux_adjust_fused: separate vh_aux,dvhdv init loops * meridional_flux_adjust_fused: 1d->2d tmp arrs * meridional_flux_adjust_fused: separate all arrs * meridional_mass_flux: separate (d)v_cor assignment loops * copy set_merid_BT_cont that accepts entire arrs * meridional_mass_flux: move j-loop into set_merid_BT_cont_fused * set_merid_BT_cont_fused: use meridional_flux_adjust_fused * set_merid_BT_cont_fused: separate tmp var init * set_merid_BT_cont_fused: separate short circuit loop * set_merid_BT_cont_fused: rm redundant k loop in short circuit * set_merid_BT_cont_fused: separate dvL/R init loop * set_merid_BT_cont_fused: make remaining tmp arrs 3d * set_merid_BT_cont_fused: use merid_flux_layer_fused * set_merid_BT_cont_fused: separate last loop * meridional_mass_flux: separate any_simple_OBC loop * zonal_edge_thickness: move k loop->PPM_reconstruction_x * meridional_edge_thickness: move k loop->PPM_reconstruction_y * PPM_reconstruction_x/y: move k loop->PPM_limit_pos * PPM_reconstruction_x/y: move k loop->PPM_limit_cw84 (untested) * zonal_BT_mass_flux: separate all loops (untested) * meridional_BT_mass_flux: separate all loops (untested) * set_zonal_BT_cont_fused: clean up var defs * use SZJB_(G) in do_I dclrns in merid* subroutines * set_merid_BT_cont_fused: clean up var defs * set_merid/zonal_BT_cont_fused -> set_merid/zonal_BT_cont * zonal_flux_adjust_fused: clean up var defs * zonal_flux_adjust_fused -> zonal_flux_adjust * zonal_flux_layer_fused -> zonal_flux_layer * meridional_flux_adjust_fused: clean up var defs * meridional_flux_adjust_fused -> meridional_flux_adjust * merid_flux_layer_fused: clean up var defs * merid_flux_layer_fused -> merid_flux_layer * fix call merid/zonal_flux_layer line conts * zonal_mass_flux: use visc_rem_u_tmp * zonal_mass_flux: separate du_min/max_CFL non-use_visc_rem update case (untested) * copy set_zonal_BT_cont with 2d args * meridional_mass_flux: use visc_rem_v_tmp * copy meridional_flux_adjust that accepts entire arrs * copy set_merid_BT_cont that accepts entire arrs * set_merid_BT_cont_fused: rm redundant k loop in short circuit * set_merid_BT_cont_fused: rm redundant k loop in short circuit * zonal_flux_adjust_fused -> zonal_flux_adjust * meridional_flux_adjust_fused -> meridional_flux_adjust * zonal_flux_layer: move GV in var dclrn * zonal/meridional_mass_flux: remove old visc_rem vars * remove redundant kloop * remove problematic omp directives * port PPM_limit_pos/CW84, PPM_reconstruction_x, zonal_edge_thickness * PPM_reconstruction_x: add enter/exit data stmts * continuity_zonal_convergence: save loop range for porting convenience * port continuity_zonal_convergence * zonal_mass_flux: array init -> loop init * zonal_mass_flux: port init loops * port zonal_flux_layer * zonal_flux_layer: add enter/exit data stmts * port zonal_flux_adjust * port set_zonal_bt_cont * zonal_flux_adjust: add more arrs in enter/exit data * set_zonal_BT_cont: add enter/exit data stmts * port zonal_flux_thickness loops * zonal_flux_thickness add enter/exit data stmts * zonal_mass_flux prepare obc for porting * zonal_mass_flux: merge any_somple_OBC loops * zonal_mass_flux port main loops * zonal_mass_flux port OBC loops * zonal_mass_flux: add enter/exit data stmts * continuity_ppm: add initial enter/exit data stmts * port PPM_reconstruction_y, meridional_edge_thickness * port continuity_merdional_convergence * ppm_reconstruction_y: add enter/exit data stmts * port merid_flux_layer * meridional_flux_adjust: port loops * meridional_flux_adjust: add enter/exit data stmts * meridional_mass_flux: port core loops * meridional_mass_flux: port OBC loops (untested) * meridional_flux_thickness: port core loops * meridional_flux_thickness: attempt port OBC loops (untested) * meridional_flux_thickness: add enter/exit data stmts * port set_merid_BT_cont loops * set_merid_bt_cont: add enter/exit data stmts * meridional_mass_flux: add enter/exit data stmts * zonal/meridional_mass_flux: add missing vars in enter/exit stmts * continuity_PPM: complete enter/exit data stmts * *_edge_thickness: do concurrent * zonal_mass_flux: do concurrent-ify * set_zonal_BT_cont: do concurrent-ify * zonal_flux_adjust: do concurrent-ify * zonal_flux_layer: do concurrent * zonal_flux_thickness: do concurrent * continuity_zonal_convergence: do concurrent * meridional_mass_flux: do concurrent * set_merid_bt_cont: do concurrent * meridional_flux_adjust: do concurrent * merid_flux_layer: do concurrent * meridional_flux_thickness: do concurrent * continuity_merdional_convergence: do concurrent * formatting * continuity_PPM: update LB * meridional/zonal_flux_adjust: separate ij-reduction * zonal_flux_adjust: a couple jki loops * set_zonal_bt_cont: some jki loops * zonal_mass_flux: some jki loops * optimise loops by using scalar zonal_flux_layer * zonal_flux_adjust: duhdu -> scalar * elemental zonal_flux_layere and separated OBC * zonal_flux_adjust: back to jki * zonal_flux_thickness: improve cpu perf a bit * zonal_flux_adjust: add some comments and guard early exit * zonal_flux_adjust: use omp target loop for private arrs * set_zonal_BT_cont: use omp target loop for private arrs * set_zonal_BT_cont: do conc jki loop * zonal_mass_flux: rm useless do_i init * zonal_flux_layere: precalc g_dy_Cu*por_face_areaU * zonal_flux_layere: precalc dh * zonal_flux_thickness: precalc dh * zonal_flux_adjust: make uhbt optional * zonal_flux_thickness: assign to outputs directly * mv zonal_flux_adjust from set_zonal_bt_cont->zonal_mass_flux * zonal_mass_flux: reuse visc_rem_u_tmp more * zonal_mass_flux: remove redundant if stmt * zonal_mass_flux: force inline zonal_flux_layere * zonal_flux_layere: inlinable gfortran -O3 but slower for ifort * zonal_flux_layere: make a bit "smaller" * zonal_mass_flux: move present(uhbt) or set_BT_cont to new subroutine * Revert "zonal_flux_layere: make a bit "smaller"" This reverts commit 316152eb3cf515c4179f8ceb738f9d259233915c. * Revert "zonal_flux_layere: inlinable gfortran -O3 but slower for ifort" This reverts commit 51ee896c54b4a349acb2ffe70489dc76558864b1. * fix gcc line truncation * pass doxygen tests * remove forceinline dirs * attempt document vars * a couple long lines * last trailing space * meridional_mass_flux: use zonal_flux_layere * zonal_flux_layere_OBC: make elemental * meridional_mass_flux: move dv_min/max_CFL calc into j-loop * meridional_mass_flux: move big chunk into separate subroutine * present_vhbt_or_set_bt_cont: merge couple of loops * meridional_flux_thickness: cpu optimize a bit * meridional_flux_adjust: back to jki * set_merid_BT_cont: pull out meridional_flux_adjst * set_meridional_BT_cont: optimize cpu * rm remaining merid/zonal_flux_layer * zonal_flux_layere: improve naming a bit * zonal_flux_layere_OBC: improve naming a bit * improve flux_elem line wraps * optimize data transfers a bit * pass elem not arr * add some missing documentation * fix trailing spaces and missing var docs * add last param docs * meridional_flux_adjust: fix fpe err * fix another fpe * cleanup args for new helper subroutines * clean up enter/exit data * fix passing h twice * meridional_mass_flux: do concurrent * zonal_flux_adjust: use a few 3d tmp arrays to mirror meridional_flux_adjust * move target update out of continuity * initialize pbv%por_face_area[U/V] on GPU * cleanup some transfers from continuity_PPM * clean up a few minor things * zonal/meridional_flux_adjust: use scalar u/v_new * declare vp,up,h_tmp on gpu * remove h update * continuity_PPM: minimize mapping stmts * zonal_flux_adjust: minimize mapping stmts * set_zonal_bt_cont: minimize mapping stmts * merional_flux_adjust: minimize mapping stmts * set_merid_bt_cont: minimize mapping stmts * zonal/meridional_flux_adjust: tmp vars duhdu/dvhdv 3d -> scalar * separate alloc of private variables for gcc * u/vh_aux: 3d->2d * target teams loop recognizable by amdflang * Continuity CS outside of init This moves the Continuity CS to the dycore init function. For some reason, this avoids an answer change with CPU. (Possibly because alloc inside of a function doesn't quite match the CS outside of it?) A few minor data transfers are also added to fix up differences in the chksum log output. * Continuity: Add locality to do concurrent Do concurrent inside of !$omp target teams loop seems to fail standard openmp tests if locality is not correctly set. This patch adds the correct locality to the four `!$omp target teams loop` directives. The domore argument has also been removed, and replaced with a `.not.any(do_i(:))` test. --------- Co-authored-by: Marshall Ward <marshall.ward@noaa.gov> * +Correct halo update sizes and reduce halo updates (#969) Added the new argument dyn_h_stencil to initialize_dyn_split_RK2 and the other 3 dynamic core initialization routines to return the size of the stencil for thicknesses as used by the dynamic core, depending on the options that are being used for the Coriolis and continuity schemes, and then used this in a set of halo updates in step_MOM_dynamics. With this change some additional halo updates that have recently been added inside of step_MOM_dyn_split_RK2 and the other 3 dynamic core time stepping routines could be eliminated. All answers are bitwise identical, but there is a new argument to 4 public interfaces. Co-authored-by: Marshall Ward <marshall.ward@gmail.com> * Correct multi-PE velocity truncation counts (#955) Modified the conditions that determine when to increment the count of velocity truncations to use the thickness as interpolated to velocity points to determine whether layers are thick enough to be counted, rather than the arithmetic mean thickness, and only count truncations that occur in the non-symmetric computational domain to avoid double counting. The filtering thicknesses should be very similar in the ocean interior, but they will differ at open boundary condition points. The corrected counting was verified by running the sloshing/layer test case with a maximum CFL set to 0.01 to create lots of truncations, and then verifying that the truncation count is now the same on 1 and 10 PEs, whereas before it was not. All solutions are bitwise identical, but the reported truncation counts in the ocean.stats files can change in multiple-PE cases with velocity truncations. * ANN parameterization of horizontal momentum eddy fluxes (#944) * Add MOM_ANN module * Mesoscale momentum parameterization with ANN - Computes subgrid stress using ANN in MOM_Zanna_Bolton - Uses MOM_ANN module for ANN inference Equivalent MOM_override for defaults ``` USE_ZB2020 = True ZB2020_USE_ANN = True USE_CIRCULATION_IN_HORVISC = True ZB2020_ANN_FILE_TALL = /path/to/ocean3d/subfilter/FGR3/EXP1/model/Tall.nc ``` * Mesoscale momentum parameterization with ANN (#2) Blank commit after squash/rebase was handled on command line * Moved MOM_ANN.F90 to src/framework/ * Minor refactor of MOM_ANN - Removed unused modules - Removed unused MOM_memory.h - Added input and output means which default to 0 and do not need to be present in the weights file - Gave defaults to means, norms, tests so that they do no need to be present in file - Added missing array notation "(:)" - Minor formatting * Adds unit tests and timing test to MOM_ANN - Added ANN_allocate, set_layer, set_input_normalization, and set_output_normalization methods to allow reconfiguration during unit tests - Added ANN_unit_tests with some simple constructed-by-code networks with known solutions - Added config_src/drivers/unit_tests/test_MOM_ANN.F90 to drive unit tests - Added config_src/drivers/timing_tests/time_MOM_ANN.F90 as rudimentary for timing inference * Adding multiple forms of inference - Adds inference operating on array (instead of single vector of features) - Implements several different versions of inference with various loop orders - Involves storing the transpose of A in the type - Tested by checking inference on same inputs is identical between variants - Added randomizers to assist in unit testing - Adds timing of variants to config_src/drivers/timing/time_MOM_ANN.F90 - Adds an interface (MOM_apply) to select preferred version of inference subroutine - Added command line args to time_MOM_ANN.F90 to allow more rapid evaluation of performance Variants explored, timed with gfortran (13.2) -O3 on Xeon: - vector_v1: - original inference from Pavel - vector_v2: - allocate work arrays just once, using widest layer - loop over layers in 2's to avoid pointer calculations and copies - speed up, x0.8 relative to v1 - vector_v3: - transpose loops - slow down, x1.54 relative to v1 - vector_v4: - transpose weights with same loop order as v1 - slow down, x1.03 relative to v1 - array_v1: - same structure as v2, working on x(space,feature) input/outputs - speed up, x0.41 relative to v1 - array_v2: - as for array_v1 but with transposed loop order - apply activation function on vector of first index while in cache - speed up, x0.35 relative to v1 - array_v3: - same structure as v2, working on x(feature,space) input/outputs - speed up, x0.58 relative to v1 * Renamed ANN variants and added some module documentation - Added module dox - Renamed _v1, _v2 etc to labels - Added ANN_apply_array_sio to ANN_apply interface - Replaced "flops" with "MBps" in timing output * Removed alternative variants of ANN in favor of optimized - Deleted variants of ANN that did not perform as well as the two versions that remain. * Apply array_sio function in ANN inference for momentum fluxes (#5) * Apply array_sio ANN inference for computation of momentum fluxes * remove trailing space * Initial commit * address Robert Hallberg code review * Restore deafult value of ZB_SCALING coefficient --------- Co-authored-by: Alistair Adcroft <Alistair.Adcroft@noaa.gov> Co-authored-by: Alistair Adcroft <adcroft@users.noreply.github.com> * Part one of thickness reservoirs - fixed the restart trouble * Better Kelvin wave results in layer mode * Response to comments on h_reservoirs PR. * Inline Harmonic Analysis Update 1: The accumulator of FtF is now associated with each FtSSH, instead of shared between all FtSSH, such that both are updated at the same model time steps, eliminating the potential inconsistency when different FtSSH are called at different time steps. * Inline Harmonic Analysis Update 2: HA_register are called inside HA_init. This should be done for all variables/fields. * Inline Harmonic Analysis Update 3: A logical flag is set to control whether harmonic analysis is to be enabled. * Inline Harmonic Analysis Update 4: time_ref and const_name are defined in HA_init, instead of being copied from MOM_tidal_forcing. This commit prepares for the separation of MOM_harmonic_analysis and MOM_tidal_forcing. * Inline Harmonic Analysis Update 5: MOM_harmonic_analysis is now independent of MOM_tidal_forcing, providing more flexibility for performing harmonic analyses on tidal constituents not available in MOM_tidal_forcing (e.g., MK3, M4). * Inline Harmonic Analysis Update 6: The frequencies of 8 overtides/compound tides (MK3, MN4, M4, MS4, S4, M6, S6, M8) have been added for the harmonic analysis. * Bug fix in MOM_open_boundary Fixed the inconsistency for defining the reference time of tides in MOM_tidal_forcing and MOM_open_boundary. * +Use I0 format to simplify integer output Use the I0 format that was introduced with Fortran 95 in 155 lines scattered across 40 files to simplify or shorten some error messages. In 21 cases, adjustl() calls that are no longer necessary for intended formatting were also eliminated. These changes have the effect of ensuring that there are still appropriate messages if there are, for example, more than 99 vertical layers or 9999 points (total) in a horizontal directions or more than 9999 PEs. In 15 cases, this change allowed for the elimination or reduction of if tests that formatted output based on the size of an integer. All answers are bitwise identical but there there may be some minor formatting changes in some error messages. * Vert friction: Fix index errors This patch fixes two class of index errors in multiple functions of `MOM_vert_friction.F90`: * `j=G%isc,G%jec` had been incorrectly applied to multiple loops. This went undetected because we almost exclusively use local indexing where `G%isc == G%jsc`, but is nonetheless a serious error. Thanks to Jorge Luis Gálvez Vallejo for reporting. * One errant loop in the shelf code had `i=is,je`. This was undetected due to poor ice shelf coverage testing. Thanks to Claire Yung for reporting. * Vert friction: Column loops moved in layers This patch moves the k-column loops inside of ji-layer loops, rather than outer-k loops of layers. The primary motivation is to restore performance at high-bandwidth runs, which were insufficiently tested during development of the k-j-i form. The inner-column loops show improved performance for both low and high-bandwidth runs. The high-bandwidth benchmark case: (128-core, 256x128 x 75 layer) ``` Profile Reference (Ocean vertical viscosity): 7.158s, 15.047s (-52.4%) ``` The low bandwidth case: (1-core, 32x32 x 75 layer) ``` Profile Reference (Ocean vertical viscosity): 3.911s, 4.788s (-18.3%) ``` For the GFDL OM5 production configuration at 503, runtimes of the slowest ranks were reduced in proportion to the high-bandwidth case above. For the reference dev/gfdl, ``` hits tmin tmax tavg (Ocean vertical viscosity) 288 4.303819 21.483670 14.452196 ``` After apply this patch, times reduce ~40% ``` hits tmin tmax tavg (Ocean vertical viscosity) 288 0.976130 13.398768 8.689331 ``` * Moving to columns allowed for removal of many `do_i` tests, since the test is applied before starting the loop. * The `touch_ij` dummy function was removed, since we're no longer trying to force an IPO optimization. * The shelf requires a re-calculation of the various thickness averages (h_arith, etc). These could be saved as 1D if it becomes a problem. * In addition to the usual regression testing, I also found no regressions in selected ice shelf configurations. * Bug fix for linear wave drag in MOM_barotropic * Linear wave drag is limited to be only applied to land points, using velocity point masks mask2dC[uv]. * Rayleigh_[uv] calculation and bt_rem_[uv] update from linear wave drag is limited for Htot>0 only. This patch eliminates potential NaN in Rayleigh_[uv] in an unusual scenario that Htot==0.0 and lin_drag_[uv]/=0. The changes do not change answers: bt_rem_[uv] is zero at land points regardless. Rayleigh_[uv] is added to [uv]_accel_bt which is masked before updating velocity. * Suppress warning message of negative eta in land In MOM_barotropic and non-Boussinesq mode, warning message on negative eta is now only issued at wet points, consistently with Boussinesq. * Flip the order of acceleration and velocity chksum In MOM_dynamics_split_RK2, now accleration chksum is printed before velocity with debug on, so that we could know which accleration term is responsible for a NaN in velocity. * vv gpu good (#44) * Vertvisc: First part of loop Moving gradually * add nvtx ranges and module warpper * port loop to openmp for now -> will move to do concurrent. Focused on porting to GPU now * more gpu * port the vertical velocity to GPUs * port the following vertical velocity loops to GPUs * more loops using omp target teams, validated * fix compute sanitizer death * meridional velcoity component ported * finish porting vertvisc * offload another loop * revert the last thing * add collapse(2) to the expensive loops * fix * fix validation * updates to porting: find coupling coef * almost done, just missingt he killer loop * lable the block of death * duct taped version of all funcs ported to GPUs * back to format and delete block constructs * remove the maps ins remnant yay! * remove need for local vertvisc_u/v variables * death loop of death has been ported * vertvisc coef does not have explicit mappings anymore * the maps are gone! * limit vel using single data transfer * do concurrent * some memory cleanup and code cleanup - profiling time * update and lots of nvtx markers for opt * limit vel betterments and notes for optimization * make limit vel good * h_ml transfer * Vertfrict: Move CS allocations to init Allocations of the vert frictions and visc control structures were moved from the dycore subroutine loops to the model initialization. Redundant grid (G) transfers were removed (although they were not triggering any transfers). Most (but not all) visc arrays have also been moved into the initialization. Kv_shear and Ray_[uv] need to be added. Trailing whitespace also (inadvertently) got removed, hopefully not a distraction. * coupling coeff fixes This appears to improve the CPU/GPU repro of the latest merge of dev/gpu into dev_gpu_vertvisc. * Resolves some of the u/v/h/dz field states between the barotropic and vertvisc functions. (Development of each assumed that the other was not complete). * This seems to fix some issues with a_[uv], Kv_tot, and dz inside of find_coupling_coef and vertvisc coef * chksum transfers were added to ensure consistent * The diff has moved to frhat[uv] in btcalc. The fields are all zero on GPU, so they are probably not transferred. Work in progress... * Remove redundant BT_cont transfer The BT_cont fields are now on GPU, so they don't need to be transferred before the btcalc call. Checksum transfers for h_ml were also added. Still assuming that it's computed on the CPU. * delete nvtcx * test to see if guarding and reintroducing the check helps * an ugly fix for a simple problem: just want to be sure before I dive * spaces * Vertvisc: formatting/style cleanup This patch reduces the formatting changes to this branch and brings it closer to both dev/gfdl and the MOM6 style guide. One minor non-cosmetic change is the reversion of ADp%dv_dt_str(i,J,1) calculation into one of the main Thomas loop. This returns it to a separate loop. * Vert friction: Explicit CS alloc in dycore init This patch removes the NVIDIA compiler preprocessing and explicitly compiles the `CS` in the dynamic core initilization. `vertvisc_init` is no longer responsible for allocation. The allocation test has also been removed. We just have to trust each other now (or segfault when it refs an unallocated field). The previous error in the CI was the absense of `CS` allocation in the other timestep methods (unsplit, unsplit RK, split RK2b). * Vert friction: Data transfer reduction Several modifications to the data transfer statements * Some redundant transfers were removed * Input/output transfers were moved up to the dycore loop * Vert friction: Cosmetic cleanup Remove some redundant comments, and an unused array loop. * Vert friction: nk_in_ml bugfix --------- Co-authored-by: Marshall Ward <marshall.ward@noaa.gov> * +Add the new parameter RESOLN_FUNCTION_OBC_BUG Added the new runtime parameter RESOLN_FUNCTION_OBC_BUG that can be set to false to take open boundary conditions into account when calculating the resolution functions at u-, v- or q-points. By default the wave speeds used to calculate resolution functions do not take OBCs into account and all answers are bitwise identical. * port main rk loop in a single commit (#47) * port main rk loop * typo and formatting * delete nan inducing copy * Dycore: NaN bugfix + minor diffs * hp upload after zero-initialization was causing random errors, likely in halo values. Now that zero-initialization is GPU-side, no need for upload. * Some nested do-concurrents were consolidated * Lots of whitespace fixes * Dycore: Mem reduction, first pass * Dycore: Mem cleanup second pass * Dycore: memcheck 3 * dycore: memcheck next * Dycore: memcheck, visc_rem and halos * dycore: Keep dz on GPU * Dycore: Remove [uv]_inst transfers * dycore: cs%eta fully on gpu * Dycore: Remove pbv copies A small change, but seems correct. pbv is input-only in step RK2, and is set outside in the main dycore loop (on CPU), followed by an upload. * dycore: memfix accel trim * Dycore: move up, vp down (and apparently an h download was redundant. looks ok, but havent 100% verified...) * dycore: Shift accels later in loop * dycore: tau[xy]_bot on GPU only * Dycore: [uv]_inst moved outside loop * Dycore: [uv]p, visc_rem_[uv] downshift * Dycore: [uvh]_av, uh, vh downshift Also some halo padding * Dycore: Remove multiple forcing/hp updates * Dycore: memfix [uv]_bc_accel * Dycore: Remove u_accel_bt and eta_pred transfers * Dycore: mem cleanup visc_rem_[uv] --------- Co-authored-by: Marshall Ward <marshall.ward@noaa.gov> * +*Fix 3-equation ice-ocean flux iteration (#972) 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 https://github.com/NOAA-GFDL/MOM6/issues/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> * Fix divide by zero fpe in apply_oda_incupd (#164) Move computation of tmp_h inside the if G%mask2dT test to avoid divide by zero error * Dycore: h/h_av mem cleanup * Dycore: Shift h, [uv]_av, [uv]h updates to halo * Dycore: static [uv]htr allocations NVIDIA was correcting for our missing uhtr allocations on the GPU in the slowest way possible. This patch adds an explicit allocation, which reduced the number of transfers dramatically. * add flag for gpu2gpu do_group_update mpi transfers for latest fms * Dycore: Move halo updates to GPU ALso: * Declate [uv]int_cor in MOM_barotropic as loop-locals * Add h to data transfer after RK2 step * DO_LOCALITY() bugfix; use GPU FMS in .testing The default FMS is not compatible with the new GPU-based MPI methods, so we just change it in .testing/Makefile for now. Down the road, we should probably move this into the .github config. DO_LOCALITY() now returns a ; if do concurrent locality modifiers are unsupported. This prevents errors associated with line continuations. * btstep: pre-calculate problematic eta checks on GPU * problem_eta->submerged,any_problem_eta->eta_is_submerged * Fixes shelfwave failure in debug mode - rotated OBC%segment%num_fields needs to be set. * Make sure reversed segments get rotated. * Bugfix for default TIDES_ANSWER_DATE in SAL 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. * Added local versions of density_elem and density_derivs without "this" argument * Fixes a typo in Recon1d PPM limiter 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. * Fixes shelfwave failure in debug mode - rotated OBC%segment%num_fields needs to be set. * Make sure reversed segments get rotated. * Frequency-dependent drag in tensor form 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. * Adds a PLM reconstruction scheme using least squares for the slope 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). * Vert friction: Force FMAs in tridiag solvers Switching the vertical friction loops from k/j/i to j/i/k replaced the evaluation of `b1` by FMA with a simpler version, causing an answer change when FMAs are enabled. Although less efficient, this patch adds an always-false loop to trick the compiler and force it to always execute `b1` by FMA. Specifically, loops of the following form execute `b1` by FMA. ``` do k=2,nz if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,k) c1(k) = dt * CS%a_v(i,J,K) * b1 b_denom_1 = CS%h_v(i,J,k) + dt * (Ray + CS%a_v(i,J,K) * d1) ---> b1 = 1.0 / (b_denom_1 + dt * CS%a_v(i,J,K+1)) d1 = b_denom_1 * b1 visc_rem_v(i,J,k) = (CS%h_v(i,J,k) + dt * CS%a_v(i,J,K) * visc_rem_v(i,J,k-1)) * b1 enddo ``` Switching to j/i/k ordering allows the Intel compiler to cache `a_[uv](K)` for use in the next iteration of `k` and evaluate `b1` by a single multiplication. If we insert an impossible branch, such as the following: ``` do k=2,nz if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,k) c1(k) = dt * CS%a_v(i,J,K) * b1 b_denom_1 = CS%h_v(i,J,k) + dt * (Ray + CS%a_v(i,J,K) * d1) b1 = 1.0 / (b_denom_1 + dt * CS%a_v(i,J,K+1)) d1 = b_denom_1 * b1 visc_rem_v(i,J,k) = (CS%h_v(i,J,k) + dt * CS%a_v(i,J,K) * visc_rem_v(i,J,k-1)) * b1 ---> if (dt < 0) exit enddo ``` then it blocks the lookahead logic of the compiler and forces the FMA execution as in the k/j/i version. There is a moderate impact on performance. ``` Before: hits tmin tmax tavg tstd tfrac grain pemin pemax (Ocean vertical viscosity) 300 2.717543 3.805039 3.523935 0.174203 0.064 31 0 511 ``` ``` After: hits tmin tmax tavg tstd tfrac grain pemin pemax (Ocean vertical viscosity) 300 2.780148 3.999669 3.761651 0.210061 0.069 31 0 511 ``` so this should only be considered a temporary fix until FMA answer changes are permitted. * fix segfault * send wt_* to gpu for btstep_timeloop * move uv_old download into if block * port find_ustar_mech_forcing interface of find_ustar * init I_Hbbl on GPU instead of CPU * clean up some loops * kji -> jki some loops * remove redundant CS%CA[u,v]_pref downloads * fuse some loops in vertical viscosity * move visc%nkml_visc_[u,v] to out of j loop * Vert frict: nk_in_ml reduce locality in v This patch adds a missing reduce(max: max_nk) locality modifier to nk_in_ml in the v-direction. This fixes an openmp regression that was detected in ifx 2025.2. * port find_uv_at_h * Corrected unit descriptions in 64 comments 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. * *Update TC testing parameters for late 2025 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 (#985) * 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 * Extend the PGF reconstruction to allow PLM-WLS 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 (#983) * 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. * Correction on total column thickness for wetting 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 * fix repeated symbol eror that ocurs in benchmark * Barotropic: frhat[uv] HYBRID repro fix This patch fixes a minor bit reproducibility regression in the calculation of `hat[uv]tot` and, consequently, `frhat[uv]` in the ifort compiler. Currently, the `hat[uv]tot` loops were moved to separate loops. This patch moves `hat[uv]tot` back into the principal `frhat[uv]` (or in this case, `hat[uv]`) loop. There may be performance consequences to this, I have not yet investigated. This was notably subtle, since many tests only call `btcalc()` once without the BTCONT `h_[uv]` inputs, was only observed in `frhatu`, and did not actually change solution answers. Nonetheless, this is a genuine answer change in an actively used solver, so we want to preserve bit repro until discussed and approved by the consortium. Although I can't see exactly the reason for the diff, it seems possible that Intel has added a reduction-like optimization, even at `-O0`. It also may not have occurred in production compilations with `-fp-model` flags. * Allow overshoot for for grounding test In commit b8c807be327c0, 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. * Initialize and integer only set on root_PE() 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. * Fix an uninitialized float in set_viscous_ML() `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. * Add floor to "h_marg" in continuity_PPM 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. * *+Fix CHANNEL_DRAG with bathymetry above sea level 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 66 unit descriptions in comments 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 (#995) * 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 * offload pass_uta_uhbta * (+) Decouple FMS infra from framework 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`. * Add mom_cap_outputlog.F90 that enables output logging diagnostics at a given hourly output frequency Author: Denise Worthen This feature is required for UFS operational configurations and is used to determine when MOM6 output (diagnostics and restart) have been completed. The log files created by this feature can be queried by the Global Workflow to either trigger downstream jobs or to ensure that if a run fails and a restart is required, model output is available consistent with a given restart file. * Revert vertical viscosity subroutines to JIK loop order used in dev/gfdl (#101) * Rewrote vertvisc_remnant to follow dev/gfdl jik structure. Uses omp directives and private vars to offload to gpu * Rewrote vertvisc to follow dev/gfdl jik structure and use omp directives + private vars in tridiagonal solver * Vertvisc with jik loops based on Marshals vv_coef branch, no transfers except ustar and hml * Switch to distribute parallel do and omp declare target to allow threading with find_coupling call inside of vertvisc_coef * Cleanup commented out code in vertvisc_coef * Revert non-k find_couling functions back to their dev/gfdl version, reintroduce conditionals to vertvisc_coef * Explicitly make coupling routines gpu routines to avoid crash when compiling with O2. Make transfers explicit * Harmonization with dev/gfdl * Several minor changes to harmonize this patch (and dev/gpu) with dev/gfdl. * touch_ij is removed since it is no longer used by the CPU. * Various whitespace changes to OMP directives and loop indices. * Unused OMP directives have been removed. * An `associated(ADp%dv_dt_visc)` typo has been fixed. * Added latent heat flux from ice shelf to ocean fluxes * Fixes wrong number of levels in z-coord diags 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. * Fix bug in registration of ALE sponge diagnostics for generic tracers (#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 * Fix reference of (rarely) unassociated pointer 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. * Adds 5 CMIP7 diagnostics for vertically integrated heat/salt content 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 * Optimized the ice-shelf CG scheme by reducing the number of times reproducing_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. * Check that frazil is allocated before adding it to ice-shelf water flux calculation. Needed for runs without frazil. * Added melt_mask for ice shelves * Added melt_mask to ice-shelf restart * comments and units * subroutine ice_shelf_solve_inner: Completed variable descriptions and units; converted cg_halo and max_cg_halo from real to integer * +Add trim_trailing_commas and ints_to_string Copied the function i2s from MOM_diag_mediator into the function ints_to_string in MOM_string_functions, and moved the code removing trailing commas from two places in MOM_diag_mediator into the new function trim_trailing_commas in MOM_string_functions. Because of the duplication of code between MOM6, SIS2 and the MOM6 ice shelf code, these functions would need to be replicated 3 or 6 times without these changes. Also added unit tests of both new functions to string_functions_unit_tests. All answers are bitwise identical but there are two new public functions in MOM_string_functions. * Call trim_trailing_commas from register_diag_field Call trim_trailing_commas from register_diag_field and register_static_field and ints_to_string from trim_trailing_commas and eliminated the now redundant routine i2s. All code functions exactly as before but there is less duplicative code. * Refactor nsten_halo in routine advect_tracer Move nsten_halo out of iteration loop * Fix OBC indexing bug in MOM_tracer_advect Fix a bug that tracers in domain outside of the OBC is falsely updated when then the OBC is in the interior. The bug was due to an indexing error in routine advect_x. * MOM_interpolate: use get_axis_size() The prior version of `get_external_field_info` incorrectly relied on the `size` output of `get_external_field_info_infra` to determine the size of an external field's axes, since all external fields are assumed to be domain-decomposed. Since axis metadata is generally opaque, we have introduced a new infra function, `get_axis_data`, which returns the size of an axis. * Adds the ability to read a CDEPS configuration file to provide in-line forcing. * Adds the ability to read a CDEPS configuration file to provide in-line forcing. Currently this is set up to read a non-climatological lrunoff data stream only. * ice-ocean-nolib: Fix SIS2 paths Patch to fix the SIS2 paths in the pipeline CI script. Explicitly excludes the icebergs stub, since we are using the actual icebergs model. * Correct the path to the Icepack interfaces The previous attempt to fix the automated no-library build of the ice-ocean model incorrectly specified the path to the Icepack_interfaces. This has now been corrected from `src/SIS2/config_src/external/Icepack_interfaces` to `src/SIS2/config_src/external/Icepack_interfaces` in pipeline-ci-tool.sh. The real mystery here is why the testing on the previous PR actually worked. * Delete unneeded masks args from 25 post_data calls Removed redundant mask arguments from 25 post_data() calls for 2-d arrays that were using masks that would have been set anyway based on the axes of these diagnostics. Explicit masks are only required for arrays that use unusual masks, pass atypically sized arrays (e.g., just the computational domain), or are static diagnostics that do not evolve in time. All answers and diagnostic output are bitwise identical. * vertvisc: add missing b_denom_1 map delete * vertvisc: remove scalar allocs * Add 2D meanSL field The spatially varying time mean sea level meanSL is used as a reference height to calculate, e.g., time mean ocean column thickness max(meanSL + bathyT, 0.0). This field allows the model run in a domain with spatically varying mean height, e.g. the Great Lakes system. This first commit insulates the changes from the rest of the model. It only adds the field to ocean_grid_type and dyn_horgrid_type, the transcription between the two types, and a routine to read it from a file. The field is not yet used by the rest of the code. * Use meanSL to calcualte mean column thickness This commit uses G%meanSL in 13 modules. The change is essentially replacing G%bathyT + G%Z_ref with G%meanSL + G%bathyT. Note that this does NOT mean parameter G%Z_ref is replaced by G%meanSL. G%Z_ref is factored in both G%meanSL and G%bathyT and it is kept as a useful consistency testing tool. Another cosmetic change is made by using G%meanSL + G%bathyT, instead of G%bathyT + G%meanSL, which (hopefully) can be easily interpreted as G%meanSL - (-G%bathyT). * Modify max_depth calculation using meanSL max_depth is really used as a maximum static thickness throughout the model, so meanSL needs to be considered. * +Fix how missing values are handled in post_data At no point does MOM6 code actually set arrays passed to the post_data() to have a missing value. Instead a missing value is set in output files entirely by masking. This commit eliminates the logic that would (inaccurately) try to reset fields that seem to match rescaled missing values to the output missing value. The previous code was inaccurate, in that a rescaled field could have taken on the unscaled missing value as a valid data point and still have been incorrectly marked is missing, although the odds of this happening are exceptionally small and it would only be cases with dimensional rescaling where this could have applied. For 2-d diagnostics, this commit eliminates a duplicative array syntax math expression that did exactly what the code now does. All solutions are identical, and because the missing value was not being explicitly it is unlikely that any diagnostics will change. * Remove ice-sheet melting/freezing contribution to fluxes%latent because it is already accounted for in fluxes%sens * Add tracing instrumentation to nuopc driver (#162) * adds calls to ufs tracing routines that will create a trace file which can then be visualized, which is found to be useful in identifying various performance issues. * +Add G%IdxCu_OBCmask and G%IdyCv_OBCmask Added the new elements `IdxCu_OBCmask` and `IdyCv_OBCmask` to the `ocean_grid_type` and `dyn_horgrid_type` to facilitate the application of no-gradient open boundary conditions at faces with essentially no added overhead. These new arrays are set initially in `set_derived_metrics()` and `set_derived_dyn_horgrid()`, but may be reset in `initialize_masks()` and `open_boundary_impose_land_mask()`. All answers are bitwise identical but there are a pair of new 2-d arrays in two transparent grid types. * Use G%IdxCu_OBCmask in 7 places Modified the code to use `G%IdxCu_OBCmask` and `G%IdyCv_OBCmask` in 7 places each in 6 modules. They are used instead of `G%OBCmaskCu*G%IdxCu` and `G%OBCmaskCv*G%IdyCv`, to which they are equivalent. This change should slightly speed up the model, and as expected all answers are bitwise identical. * Add option to scale tidal amplitude for bottom ustar. (#1016) * Add option to scale tidal amplitude for bottom ustar. - previously we used the tidal amplitude to compute ustar. - The additional factor translates between amplitude and time mean tidal current. - Setting the factor TIDEAMP_FACTOR<0 preserves old answers. * Update tideamp factor implementation for efficiency - factor out the negative "default" value to automatically set to multiply by 1.0 instead of using an if-block. - factor in the c-grid averaging 0.5 to further reduce extra operations, but clearly label the parameter to reflect this. --------- Co-authored-by: brandon.reichl <brandon.reichl@noaa.gov> * Add vertical tracer flux diagnostic for dye tracers (#1022) * Add vertical tracer flux diagnostic for dye tracers - Register vertical flux diagnostic in initialize_dye_tracer - Calculate net vertical flux from entrainment (positive upward) - Post flux diagnostic in dye_tracer_column_physics * changed diagnostic registration to be at interface, made sure boundary fluxes are zero * changed lines 338 and 354 as needed. Fixed accidental space on Line 1. * Regroup MOM_initialize_fixed params in param_doc This commit is meant to fix the issue that all parameters in MOM_initialize_fixed after OBC are logged under module MOM_open_boundary in MOM_parameter_doc. By moving log_version call after OBC, parameters from MOM_initialize_fixed are now logged under three "modules" in MOM_parameter_doc: 1. Parameters before OBC are under module MOM_grid_init, which also (incorrectly) includes topography relatd parameters. 2. module MOM_open_boundary 3. Parameters after OBC are under module MOM_initialize_fixed. The change makes sure OBC parameters are well separated from the other parameters. This is a hack rather than a fix. * Minor open_boundary_config refactor * Make OBC related calls in MOM_initialize_fixed explicitly conditional for readibility. * Early return in open_boundary_config if there is no se…

MJHarrison-GFDL reviewed Sep 6, 2019

View reviewed changes

marshallward mentioned this pull request Sep 9, 2019

Low CodeCov coverage threshold #996

Merged

kshedstrom and others added 3 commits September 16, 2019 11:04

More OBC cleanup.

f4a4aae

Deleted unwanted commented OBC code.

51fda4d

Merge branch 'dev/gfdl' into dev/esmg

53b8be5

adcroft merged commit 53b8be5 into mom-ocean:dev/gfdl Sep 17, 2019

claireyung added a commit to ACCESS-NRI/MOM6 that referenced this pull request Dec 17, 2025

Revert "Fix for ice-shelf friction velocity bugs (mom-ocean#995)"

8c993a6

This reverts commit 401ec5b.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added option to turn off baroclinic OBC time filtering#995

Added option to turn off baroclinic OBC time filtering#995
adcroft merged 3 commits into
mom-ocean:dev/gfdlfrom
ESMG:dev/esmg

kshedstrom commented Sep 6, 2019

Uh oh!

codecov-io commented Sep 6, 2019 •

edited

Loading

Uh oh!

MJHarrison-GFDL left a comment

Uh oh!

kshedstrom commented Sep 6, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Conversation

kshedstrom commented Sep 6, 2019

Uh oh!

codecov-io commented Sep 6, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

MJHarrison-GFDL left a comment

Choose a reason for hiding this comment

Uh oh!

kshedstrom commented Sep 6, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov-io commented Sep 6, 2019 •

edited

Loading