Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GFDL to main (2024-11-27) #1646

Closed
wants to merge 125 commits into from
Closed

Conversation

marshallward
Copy link
Collaborator

Celebrating American Thanksgiving with a new PR from GFDL 🦃. This patch includes new schemes to support eddy backscatter parameterization, tidal modeling, internal tide propagation, and shear mixing near shelves. It also includes many improvements and bugfixes to existing solvers.

Features and new capabilities

Updates with minor improvements or new diagnostics

Refactoring

Bugfixes (including self-consistency corrections)

Infrastructure

Contributors

kshedstrom and others added 30 commits June 2, 2024 05:55
* Spatially variable fields for MLE%Bodner

 - Allows reading in 2d fields for Cr and for MLD_decaying_Tfilt.

* Finish the job?

* Renamed one variable, fixed some units
* Added option to convert flux through static ice front into icebergs

Added variables for the accumulated iceberg mass and heat flux due to
calving from ice shelves (flux through the static ice front). These
will be passed to the coupler and SIS2/iceberg module to initialize
bergs. Also fixed the ice-shelf SMB override and reorganized
ice-shelf post data calls so that they do not strictly have to be
called at multiples of the ice velocity time step.

* Added ice-shelf scalar diagnostics

Added ice-shelf scalar diagnostics related to volume-above-floatation
and surface/basal mass balance. Had to modify ice-shelf diag mediator
to allow scalar diagnostics.

* Fixed units for volume-above-floatation and Cp_ice. Renamed volume-above-floatation variable from 'vab' to 'vaf'

* Fixed write_ice_shelf_energy call within subroutine solo_step_ice_shelf so that it is passing the correct arguments

* Fixed syntax of calving units

* Fixed units for ice shelf calving and scalar diagnostics
  Refactored the Idealized_Hurricane module to clean up strange or unscaled
variable units and eliminate dimensional scaling factors with the latest answer
dates.  This includes the introduction of 18 new runtime variables to replace
hard-coded dimensional constants and two runtime logical flags to reproduce
existing bugs.  An inconsistency in the sign convention for the distance from
the hurricane center with idealized_hurricane_wind_forcing in (the probably not
yet used) single column mode was also corrected. Also added descriptions with
units for all the variables in this module.  By default or with appropriate
parameter settings all answers are bitwise identical.
* Add rescaling paramter to KD Shear

Add a parameted (beta) to rescale the distance to the nearest solid
boundary that is used within calculation of kappa shear. While rescaling
this value is unphysical, this adjustment can be thought of as not
allowing shear instabilities to grow up to the full distance to the
nearest boundary because of other turbulent processes which would disrupt
their growth.

* Rename parameter and swtich to multiplication

Rename the parameter beta to lz_rescale to be more descriptive.
To avoid two divisions in a single line, the inverse of lz_rescale
is calculated and stored as I_lz_rescale_sqr. Where the calculation
is done logic has been added to check that lz_rescale is greater
than zero to avoid dividing by zero or making lz_rescale negative.

* Add Comment where lz_rescale is used

Based on Wei's suggestion, add a comment where lz_rescale is used
and put the multiplication by I_lz_rescale to the next line.

* Remove Trailing whitespace

---------

Co-authored-by: Theresa Morrison <[email protected]>
Co-authored-by: Theresa Morrison <[email protected]>
Co-authored-by: Theresa Morrison <[email protected]>
Co-authored-by: Theresa Morrison <[email protected]>
fix style errors

Modify axes_data to make it allocatable

fix style
  Added the new runtime parameters ROBUST_STOKES_PGF and LA_MISALIGNMENT_BUG to
allow for the selection of Stokes pressure gradient force calculations that work
properly in the limit of vanishingly thin layers and to correct a sign error in
the calculations of the misalignment between the waves and shears in the
Langmuir number calculations when LA_MISALIGNMENT is true.  By default, all
answers are bitwise identical, but as these options are not yet widely used it
might make sense to move aggressively to obsolete the previous code once more
extensive testing has take place.  There will be new entries in the
MOM_parameter_doc files for some cases, but there are not yet any such cases
in the MOM-examples regression suite.
-Account for re-entry boundary conditions and non-symmetric grids
-Optimized SSA CG scheme to eliminate unneeded loops and pass_var calls
-Added a missing pass_var call that should fix a reproducibility issue
 on different PE layouts
-Automatically adjust ice-shelf velocity data when it is read in from
 file so that it agrees with assigned BCs
-Changed some ice-shelf mass units that were affecting partially-fully
 cells at the ice front, so that they are always mass per ice-shelf
 area (not per grid area)
  This commit improves the documentation of the remaining real variables with
previously undocumented units in the framework directory, with some other
similar unit documentation improvements in other files.  Units were added to
comments describing about 79 real variables in 4 framework modules,
atmos_ocean_fluxes.F90 for 2 drivers, MOM_oda_incupd.F90, MOM_oda_driver.F90
and MOM_variables.F90.

  The variable nhours_incupd in initialize_oda_incupd was renamed to
incupd_timescale and the factor of 3600.0 that converts ODA_INCUPD_NHOURS
from hours into seconds was moved from where this variable is used into the
scale argument where it is set to help document the meaning and units of this
variable as simply and clearly as possible.

  The unused variable smb was removed from Shelf_main in ice_shelf_driver.F90.
The internal variable tmp in RGC_initialize_sponges was renamed rho to reflect
its contents, and its contents and units are now described in a comment.

  Although the units have been added to the description in MEKE_vec as though it
is being used in a dimensionally consistent way, I suspect that this might
actually be in MKS units, in which case a unit conversion factor might be needed,
and a comment has been added to note this.  However, the machine learning code
that is used to set this array comes from an external package that is not being
used yet at GFDL, so it is not clear what code should be examined to address
this question.

  A comment was added in set_up_global_tgrid noting a unit rescaling factor
that appears to be missing from a dimensional constant.

  All answers are bitwise identical, and for the most part only comments are
changed, although two internal variables were renamed.
  Corrected the sign convention used for the neutral_slope_x and neutral_slope_y
diagnostics in cases when there is not an equation of state being used to match
the usual sign convention (interfaces sloping up to the east or north is a
positive slope) and to match what is done when an equation of state is being
used to calculate neutral density slopes rather than just basing the slopes of
the interfaces alone.  The same correction was made to the slope_x and slope_y
arguments returned from calc_isoneutral_slopes, and self-consistent changes were
made to the overturning streamfunction calculations in thickness_diffuse_full.

  In addition, there were some corrections made to the documentation at the end
of MOM_thickness_diffuse.F90, and the calculate_slopes argument to
calc_slope_functions_using_just_e() that was always hard-coded to .true. was
eliminated to avoid any confusion about whether these changes might propagate
beyond the interface height diffusion calculations.  This commit addresses most
aspects of the issue discussed at github.com//issues/359, although
it does not change the sign convention for the overturning streamfunction.

  All solutions are bitwise identical, but there is a change in the sign
convention for two diagnostics and two arguments to a publicly visible interface
in pure layer mode when no equation of state is used.
  Added the new runtime parameter DETERMINE_TEMP_CONVERGENCE_BUG to avoid using
layout-dependent tests for temperature and salinity tolerances to determine when
to stop iterating in determine_temperature.  Also added logic that correctly
determines when iterations can be stopped because no further changes occur, and
rearranged to loops to avoid revisiting layers that have converged.

  Because the default tolerances for the temperature, salinity and density
changes are set to be the same and the partial derivatives of density with
temperature and salinity are less than 1 in the MKS units used here for a
realistic equation of state of seawater, this bug does not impact any of the
existing regression test cases, but this bug could lead to non-reproducibility
with non-default values.

  This commit changes the entries in the MOM_parameter_doc files that have
Z_INIT_ALE_REMAPPING = .false. and FIT_TO_TARGET_DENSITY_IC = .true., by adding
a new runtime parameter and by not logging DETERMINE_TEMP_T_TOLERANCE and
DETERMINE_TEMP_T_TOLERANCE if they are not used because
DETERMINE_TEMP_CONVERGENCE_BUG is false.  By default, all answers are bitwise
identical.
  Corrected the field being posted for the Kd_interface diagnostic inside of
diabatic_ALE() to be the minimum of Kd_heat and Kd_salt, which includes the
contributions from ePBL, but excluding any extra double-diffusive mixing, making
it analogous to the diagnostic that is currently being posted from
layered_diabatic() or diabatic_ALE_legacy().  Diabatic_ALE() is used when
USE_REGRIDDING=True and USE_LEGACY_DIABATIC_DRIVER=False, in which case a
diagnostic will change, correcting the issue noted at
/issues/398.  All solutions are bitwise identical.
  This commit adds the necessary information about the turns of the model grid
relative to the (unrotated) coupler_type data fields that are inside of the
forcing type and surface_state type and are used with passive tracers, so that
certain tracer packages can now be used with rotated grids.  The previous
version had problems where the model would not run when ROTATE_INTEX = True and
the CFCs (and probably other passive tracers) were being used, as noted at
/issues/621.  These problems have now been fixed.

  There are new calls to coupler_type_spawn() in allocate_forcing_by_ref() and
allocate_surface_state() that manage the creation of the coupler_2d_bt_types for
unrotated types based on information from their rotated counterparts.

  There is a new optional turns arguments to allocate_forcing_by_ref() and new
optional sfc_state_in and turns arguments to allocate_surface_state(), and these
are now being used in at least 6 places where unrotated temporary surface_state
or forcing types are being set up.

  There are also new optional turns argument to extract_coupler_type_data() and
set_coupler_type_data() that are used to deal with the fact that the internal
data arrays in the coupler_bc_types are never rotated, unlike the other MOM6
arrays, because they have to be passed directly into the generic tracer code.
These new turns arguments are used in 14 calls from various tracer packages,
including in 6 calls from the OCMIP2_CFC code.

  There are 4 new calls to deallocate_surface_state() or
deallocate_forcing_type() that were added to avoid (presumably minor) memory
leaks when grid rotation is enabled.  These new calls are in
initialize_ice_shelf_fluxes(), shelf_calc_flux() and in the surface flux
diagnostics section of step_MOM().

  All answers are bitwise identical in any cases that worked previously, but
some cases with grid rotation that previously were failing with cryptic error
messages are now running successfully.  There are new optional arguments to
several publicly visible routines.
The "eta has dropped below bathT" warning message has indices off by 1.
This is because G%isd_global = G%isd + HI%idg_offset and G%isd is
(most of the time) 1.

G%isd_global/G%jsd_global are now replaced by
G%HI%idg_offset/G%HI%jdg_offset.
Add an option to normalize wt_[uv] in the barotropic solver. This is
fix step 1 of 2 to address the issue that visc_rem is applied
incorrectly to the barotropic solver.

A runtime flag BAROTROPIC_VERTICAL_WEIGHT_FIX is added to control the
behavior of this commit. The current default is false (not applying the
fix).
Add an option to use `dt` instead of `dt_pred` in the calculation of
vertvisc_remnant() at the end of predictor stage. The change affects
the following `continuity`() call and `btstep`() call in corrector step.

Combined with the changes from the previous commit, the new options
should solve the issue of inconsistent barotropic velocities from
barotropic solver and baroclinic step.

A runtime flag VISC_REM_TIMESTEP_FIX is added to control the behavior of
 this commit. The current default is false (not applying the fix).
This commit is the third fix related to visc_rem in split mode. h_[uv]
from `zonal_flux_thickness`() and `meridional_flux_thickness`() is
currently multiplied by `visc_rem_[uv]`. This seems to be
double-counting as `visc_rem_[uv]` is accounted for in the barotropic
solver vertical weights. Moreover, for options other than BT_cont, the
velocity point thickness does not include visc_rem. The fix removes the
visc_rem factor from h_[uv] when BT_cont is used.

A runtime flag VISC_REM_CONT_HVEL_FIX is added to control the behavior
 of this commit. The current default is false (not applying the fix).

Other small changes:
* Write out explicitly all keyword arguments in `continuity`() calls
in MOM_dynamics_split_RK2.
* Rename a runtime parameter BAROTROPIC_VERTICAL_WEIGHT_FIX from a
previous commit to  VISC_REM_BT_WEIGHT_FIX, to be consistent with the
style of the other two flags.
* Add VISC_REM_TIMESTEP_FIX parameter to MOM_dynamics_split_RK2b
* Correct a bug using dimensional h_neglect in calculating
non-dimensional Iwt_[uv]_tot, by replacing the division with
an Adcroft_reciprocal style division.

* Add a parameter VISC_REM_BUG to control the defaults of all three
viscous remnant bug related parameters. The individual parameters should
 be removed in the future.
  Refactored PressureForce_FV_nonBouss and PressureForce_FV_Bouss to use more
3-dimensional arrays in preparation for the changes that we are expecting from
Claire Yung to reduce the pressure gradient errors in ice shelf cavities.  These
anticipated changes will involve selecting an internal interface at which to
define the shape of the interfaces with depth when there is an ice shelf, rather
than assuming that the top surface pressure varies linearly with depth between
the two top corners.  In PressureForce_FV_nonBouss, this refactoring includes
turning za, intx_za and inty_za into 3-d arrays and moving the code setting
these arrays outside of the k-loop setting the pressure gradient forces.
Changes in PressureForce_FV_Bouss are analogous but involve turning 7 arrays
(pa, dpa, intz_dpa, intx_pa, intx_dpa, inty_pa and inty_dpa) into 3-d arrays.
In both cases, the code for a reduced gravity model is consolidated into a
single block toward the end of the routines.  These changes to use more 3-d
arrays could increase the computational time for the pressure gradient
calculations, but in short regression tests any such increase in computational
time appears to be small.  No external interfaces are changed, and all answers
are bitwise identical.
  It was noted in /issues/491 that the call to
initialize_regridding() for diagnostics gives a segmentation fault when the
maximum depth of the ocean is greater than 9250 m unless DIAG_COORD_DEF_Z is
explicitly set.  This commit fixes this problem by (1) adding an extra layer to
the bottom of the hard-coded WOA09 list of diagnostic depths when the ocean is
deeper than the 9250 m maximum depth in that file, and (2) changing the default
behavior for diagnostic grids to be be uniform when the maximum ocean depth is
deeper than 9250 m, for which WOA09 is no longer likely to be a convenient
choice.  With this change, MOM6 no longer has segmentation faults for the
configuration described in issues/491.  All solutions are bitwise identical, but
some z-space diagnostic grids will be modified in cases that previously had
segmentation faults.
  Added the new runtime parameter BACKSCATTER_UNDERBOUND that can be set to
false to avoid a potential source of numerical instabilities from excessively
large biharmonic viscosities due to overly generous bounds where there are
negative Laplacian viscosities due to backscatter parameterizations.  Setting
this to true reproduces the previous solutions, while setting it to false treats
negative and zero Laplacian viscosities the same way when bounding the
biharmonic viscosity.  When the Laplacian viscosities are all non-negative, the
value of BACKSCATTER_UNDERBOUND makes no difference.  The default is true for
historical reasons, but this option probably should be set to false to avoid
certain numerical instabilities from the horizontal viscosities.  By default,
all answers are bitwise identical, but there are new entries in
MOM_parameter_doc files for cases that set both BETTER_BOUND_KH and
BETTER_BOUND_AH to true.
Add FRICTWORK_BUG parameter (default is true). If FRICTWORK_BUG is false, use the new way to compute FrictWork using the thickness flux. The previous version (realized by setting FRICTWORK_BUG as true) uses the velocity to compute FrictWork, which overestimates the total energy dissipation rate compared with the domain integral of KE_horvisc computed in MOM_diagnostics.
Bracket counter for <..> was correct but (..) was broken, it was
counting ( but not ).  I can only presume that such cases were very
rare.

This patch fixes the closed parenthesis count.
- Addressed compiler warnings about unused variables and
  argument intent.
- Note there is one dummy argument that was unused and the best way
  to address the warnings was to remove the argument, i.e. and API change.
  This was only in a debugging/utility function so does not impact the
  rest of the model.
- Address warnings about unused/uninitialized variables in MOM_remapping.F90
  - Sets values that appear to the compiler to be potentially unset, but
    always are in practice
- Added new driver, time_MOM_remapping, to config_src/timing_tests/
  that exercises the remapping_core_h() function with PCM and PLM.
- We should add other reconstruction schemes too, but the main reason for
  adding this now is to monitor the impact of refactoring the function
  remapping_core_h() which is mostly independent of the reconstruction
  schemes.
- Fixed .testing/Makefile to not fail with `make build.timing -j`
Added a simple driver for the remapping unit tests.
- Added a simple class within MOM_remapping.F90 to greatly abbreviate
  the lines comparing values of arrays.
- Translated the existing remapping unit tests to use the testing
  class.
- The first block of remap_via_sub_cells() computes the intersection between
  two columns (grids). I've split out this block into a stand alone function
  so that we can re-use it in a to-be-written analog of remap_via_sub_cells()
  that will remap integrated quantities.
- Added some in-line documentation of algorithm used by remap_via_sub_cells()
- Added new column-wise function intersect_src_tgt_grids()
- Added tests of intersect_src_tgt_grids() that check against known results
  - Had to add a new helper function to MOM_remapping to report state of
    integer arrays.
- This second phase of remap_via_sub_cells() maps the values from the source
  grid to the sub-grid using the pre-calculated arrays/indices from
  intersect_src_tgt_grids().
  This phase as-is is only used for remapping of intensive variables but
  it does implicitly calculate intermediate extensive quantities which are
  analogous to what we need when remapping momentum. Splitting this phase
  out primarily allows us to test this bit of code and it has indeed revealed
  an edge case that fails.
- Added new column-wise function remap_src_to_sub_grid()
- Added tests of remap_src_to_sub_grid() that check against known results
  - One test is commented out corresponding to an edge case that reveals
    the current code is wrong. Will enable this test after checking how
    many experiments are impacted but the associated bug fix.
marshallward and others added 28 commits September 17, 2024 15:20
Adding null read operations so that compilers will not warn about unused
input variables in dummy functions.
This is a minor patch which allows the code to be compiled under F2023
standardization.

* Arguments are declared in the order that they are used.  In this case,
  it means that array lengths are declared before the arrays using them.

* The syntax of various IO statements has been cleaned up.
Some functions in the FMS cap used arguments which depended on other
arguments, which were declared out of order.

* ocean_model_data2D_get
* ocean_model_get_UV_surf

This patch moves the array index size definitions before the array
definitions.

This is required for language standard compliance.
This patch enables Fortran 2018 standard compliance testing.  We could
do 2023, but our current CI of choice doesn't yet have a compiler which
can do this.

The actual content of this PR is a decoupling of FCFLAGS_DEBUG and
FCFLAGS_FMS.  There is now a default FCFLAGS macro which is used by the
other two macros.  One can now optionally configure FCFLAGS_DEBUG
without worrying about the impact on FCFLAGS_FMS.

The motivation here is that we don't want to test for F2018/2023
compliance in FMS.
Temporarily cut off the control of VISC_REM_CONT_HVEL_FIX from
VISC_REM_BUG and change the default of VISC_REM_CONT_HVEL_FIX to False.
…l ice-sheet velocities (#714)

* Updated with dev/gfdl. Then, added parameters for min ice thickness, min basal traction, max surf slope, and min ice viscosity to use for ice dynamics

* noted how the MIN_BASAL_TRACTION parameter input is in units [Pa m-1 yr], but converts automatically to [Pa m-1 s] in the code

* Added separate scalar ice-shelf diagnostics for Antarctica and Greenland

* Added sx_shelf and sy_shelf as ice shelf diagnostics to save the surface slope fields used in the shallow shelf approximation (SSA). This is particularly helpful for determining whether unrealistic velocities are caused by unrealistically steep surface slopes, which can sometimes arise for example, on coarse grid cells that cover both a steep mountainous region and a realively flat ice shelf. Then, the MAX_SURFACE_SLOPE parameter can be tuned to set an upper bound on the SSA surface slope to avoid these steep-slope-induced problematic velocities.

* Fix doxygen errors for ice-sheet process_and_post_scalar_data routine

* FMA fix associated with enforcement of max allowed ice-shelf surface slope
* EBT Backscatter

Backscatter code using the equivalent barotropic (EBT) mode documented in Yankovsky et al. (2024).

* Modifications to MOM_hor_visc:
- Separate FrictWork and FrictWork_bh loops
- Simplified the computation for MEKE%mom_src and MEKE%mom_src_bh when MEKE%backscatter_Ro_c /= 0.

(cherry picked from commit 8ffc6a8)

---------

Co-authored-by: Wenda Zhang <[email protected]>
Co-authored-by: Marshall Ward <[email protected]>
bigthetao refers the the uppercase greek letter theta used for
conservative temperature. Salinity so uses the roman alphabet
and a bigso is ambiguous as many papers use S for salinity hence
absso seems like the most unambiguous.
This commit fixes a bug that target code cannot be built for regression
tests when BUILD is specified with command line arguments. The fix
explicitly specifies BUILD argument for target codebase, avoiding the
variable being passed from the parent make.
…es, and ice speed; and the magnitude of ice surface slope and driving stress. Fixed some misleading units for variables related to ice viscosity
…ics. Also fixed incorrect loop bounds for a seldom-used ice-shelf convergence criterion
* Fixes to preserve rotational symmetry when using Fused-multiply-adds

Tested in ice-ocean mode and ice-shelf-only mode using ISOMIP and
MISOMIP+, respectively

* Added more parentheses to preserve rotational symmetry when using FMAs.

Note that the changes here are to the MOM_ice_shelf_dynamics.F90 function 'quad_area', which is only called by the subroutine 'bilinear_shape_functions'; this subroutine is currently unused so that rotationally-consistent use of FMAs in both 'quad_area' and 'bilinear_shape_functions' has not been tested
Fix a bug that layer/interface weights may not be reset to 1.0 if no
calculations are needed. The bug has negligible impact for barotropic
applications but may affect baroclinic runs which are yet to performed.
Other than staying up to date, this purportedly fixes a "Node.js deprecation"
message we first noticed in the coverage workflow but has appeared elsewhere
intermittently.
- Changed gitlab runner to the mom6-account on gaea
- Added gitlab variable MOM6_RUN_JOB_DURATION to control the allowed run
  duration during bad days for the system. Defaults to 15:00 (15 mins)
- Added FC=ftn MPIFC=ftn CC=cc environment vars when invoking make
  in .testing
  Eliminated the previously optional h_neglect and h_neglect_edge arguments to
remapping_core_h and remapping_core_w, and added optional h_neglect and
h_neglect_edge arguments to initialize_remapping for storage in the remapping
control structures.  The answer_date, h_neglect and h_neglect_edge arguments
to ALE_remap_scalar were also eliminated, as they are no longer used.
The new routine open_boundary_setup_vert was added to manage these changes.
A new verticalGrid_type argument was added to wave_speed_init.

  The h_neglect argument to 29 routines was made non-optional, because there is
no way to provide a consistent hard-coded default for a dimensional parameter.
The (mostly low-level) routines where this change to a non-optional h_neglect
was made include  build_reconstructions_1d, P1M_interpolation,
P3M_interpolation, P3M_boundary_extrapolation, PLM_reconstruction,
PLM_boundary_extrapolation, PPM_reconstruction, PPM_limiter_standard,
PPM_boundary_extrapolation, PQM_reconstruction, PQM_limiter,
PQM_boundary_extrapolation_v1, build_hycom1_column, build_rho_column,
bound_edge_values, edge_values_explicit_h4, edge_values_explicit_h4cw,
edge_values_implicit_h4, edge_slopes_implicit_h3, edge_slopes_implicit_h5,
edge_values_implicit_h6, regridding_set_ppolys, build_and_interpolate_grid,
remapByProjection, remapByDeltaZ, and integrateReconOnInterval.

  In two modules (MOM_open_boundary and MOM_wave_speed), two separate copies of
remapping_CS variables were needed to accommodate different negligible
thicknesses or units in different remapping calls.

  In those cases that also have an optional h_neglect_edge argument the default
is now set to h_neglect.  Improperly hard-coded dimensional default values for
h_neglect or h_neglect_edge were eliminated in 16 places as a result of this
change, but they were added back in 2 unit testing routines (one of these is
archaic and unused) that do not use exercise dimensional rescaling tests.

  Because the previously optional arguments were already present in all calls
from the dev/gfdl or main branches of the MOM6 code, and because the places
where arguments have been removed they are balanced by equivalent arguments
added to initialize_remapping, all answers are bitwise identical in the
regression tests, but this change in interfaces could impact other code that is
not in the main branch of MOM6.
* subroutine initialize_MOM
* MOM_input in test cases
This patch replaces the `CS%[uv]_file < 0` checks with
`CS%[uv]_file == -1`.  FMS1 returns negative file handles for missing or
otherwise error-prone files, but the FMS2 IO framework relies on
`newunit=` to autogenerate handle IDs, which are always negative and
cannot be used with checks for negative values.

The check is replaced with equality with -1.  `newunit` is guaranteed to
never return -1 for a valid file, so this is a valid check for a missing
file.  It also lets us continue to use -1 as the initial (unopened)
value.

Behavior is compatible with `mpp_open()` output, so this can also be
used with the FMS1 API.

A better solution would be to introduce some validation function which
is defined by each API, but there is not yet any need for such
sophistication.
This commit contains two changes to gcov report generation and codecov
upload.

1. Separation of the file parser unit tests from the others was causing
   them to be exluded from report.cov.unit.

   This patch reworks the rules to replace
   MOM_file_parser_tests.F90.gcov with driver code,
   build/unit/test_%.F90.gcov.

   This does assume at all drivers will look the same (test_*.F90) but
   that part can be reworked if it ever becomes a problem in the future.

   Thanks to @adcroft for multiple suggestions in this PR.

2. Github appears to internally store all its repositories in another
   directory with the name as the repo (in this case MOM6/).

   Normally this is hidden to everyone, but it was causing some
   confusion with the codecov upload tool, and was unable to match the
   source code to the .gcov report.

   The .codecov.yml config file was modified to adjust for this path
   change, and should now correctly allow coverage to be reported
   alongside the file.

   (The GitHub Actions app likely makes this adjustment, but we need to
   do it manually since we upload directly to Codecov.io.)
This particular bugfix was proved to be unnecessary and therefore
incorrectly added. Runtime parameter VISC_REM_BUG now controls two
bugfixes, related to RK2 time stepping and barotropic vertical weight,
which are valid.
VISC_REM_TIMESTEP_FIX -> VISC_REM_TIMESTEP_BUG
VISC_REM_BT_WEIGHT_FIX -> VISC_REM_BT_WEIGHT_BUG

The "signs" of the flags are flipped accordingly.
…lt_fingers

- Double diffusion code was updated to include runtime parameter for maximum density ratio, but it was never implemented in the actual calculation of the double diffusion diffusivities.  It is now used in the calculation and Rrho0 is removed from the code.
- The default value of Max_Rrho_salt_fingers is changed from 2.55 to 1.9 to avoid changing behavior from the old code.  However, it should be updated to 2.55 and noted that it will change answers for any run with double diffusion enabled through MOM6.  It does not affect the behavior of double diffusion enabled through CVMix.
Most modules that use remapping have their own parameter e.g.
Z_INIT_REMAPPING_SCHEME.  However, OBC, ODA and SPONGE currently
sliently use the primary ALE REMAPPING_SCHEME parameter.

Added logged parameters OBC_REMAPPING_SCHEME, ODA_REMAPPING_SCHEME
and SPONGE_REMAPPING_SCHEME to allow customization.  All take
REMAPPING_SCHEME as their default to maintain compatibility.

For ODA, also added ODA_REMAPPING_USE_OM4_SUBCELLS that takes
REMAPPING_USE_OM4_SUBCELLS as its default to maintain compatibility,
and ODA_BOUNDARY_EXTRAP that takes BOUNDARY_EXTRAPOLATION as its
default to maintain compatibility.  Note that BOUNDARY_EXTRAPOLATION
is a bug, it should have been REMAP_BOUNDARY_EXTRAP, but the former
remains the default for compatibility.

For SPONGE, also added SPONGE_BOUNDARY_EXTRAP that takes
BOUNDARY_EXTRAPOLATION as its default to maintain compatibility.
Note that BOUNDARY_EXTRAPOLATION is a bug, it should have been
REMAP_BOUNDARY_EXTRAP, but the former remains the default for
compatibility.

Answers are unchanged, but there are new logged parameters.
Added SQG vertical structure in Varmix to provide vertical profile for diffusivities
- added function calc_sqg_struct in MOM_lateral_mixing_coeffs to compute sqg_struct
- added sqg_expo to set the exponent of sqg_struct
- to use sqg_struct for the backscatter, set BS_use_sqg=true, sqg_expo>0., and BS_EBT_power=0.
- if SQG_USE_MEKE=True, use the eddy length scale from MEKE to compute sqg_struct
- added eddy length scale Le in MEKE if SQG_USE_MEKE=True
- added MEKE%Le into restart file if SQG_USE_MEKE=True
- added MEKE in Varmix
- registered N2_u and N2_v diagnostics when SQG_EXPO>0

(cherry picked from commit 6d3df0541c33d6f6d1f9fcb695f1a1eb961ec1b3)

* Compute vertical structures for khth, khtr, backscatter, and kdgl90 all in VarMix

- Vertical structures including khth_struct, khtr_struct, BS_struct, and kdgl90_struct are now computed in VarMix
- Each diffusivity/viscosity have two vertical structure options, equivalent barotropic (EBT) and surface quasigeostrophic (SQG) mode structures
- KHTH_USE_EBT_STRUCT, KHTR_USE_EBT_STRUCT, KDGL90_USE_EBT_STRUCT and BS_EBT_POWER parameters, which already existed, still control whether to use the EBT structure for khth, khtr, kdgl90, and backscatter, respectively
- Added KHTH_USE_SQG_STRUCT, KHTR_USE_SQG_STRUCT, KDGL90_USE_SQG_STRUCT and BS_USE_SQG parameters to control whether to use the SQG structure for khth, khtr, kdgl90, and backscatter, respectively
- If neither EBT nor SQG is called, no vertical structure will be used for that diffusivity/viscosity
- An error will be called if both EBT and SQG structures are called for the same diffusivity/viscosity
- Added dt as an input of calc_resoln_function. dt is needed for calc_sqg_struct called in calc_resoln_function

Bug fixes
- Changed ()**0.5 to sqrt in sqg_struct, which solves the dimension error
- Added if statement for using BS_struct in MOM_hor_visc
- Added if statement for CS%sqg_struct<=0.
- Changed the subroundoff of Coriolis parameter to 1e-40 in calc_sqg_struct
- Changed parameter BS_USE_SQG to BS_USE_SQG_STRUCT

---------

Co-authored-by: Alistair Adcroft <[email protected]>
  When HOR_VISC_ANSWER_DATE was introduced to replace HOR_VISC_2018_ANSWERS on
August 4, 2022 as a part of github.com//pull/179, the logic was
incorrectly specified, using the older form with the newer answer date and vice
versa, because `(CS%answer_date > 20190101)` was used instead of
`(CS%answer_date < 20190101)`.  (Curiously, using exactly 20190101 actually
gives the intended result.)

  This commit modifies this logic so that the older (mildly dimensionally
inconsistent) version is now being used for answer dates between 20190102 and
20241201, but a very late answer date uses the corrected form.  The offending
block of code is only used when USE_MEKE is true and the Rossby number scaling
of the biharmonic energy source is enabled by setting MEKE_BACKSCAT_RO_C > 0,
which does not appear to be very common.  To avoid logging the description of
this ugly new logic in MOM_parameter_doc files where it does not impact the
solutions, new logic was added to limit the logging of  HOR_VISC_ANSWER_DATE.
Also, as there are no known non-Boussinesq cases with this combination of MEKE
parameters, the minimum value of HOR_VISC_ANSWER_DATE was changed to 20241201
when the model is in non-Boussinesq mode.   Because this bug went undetected
when it was first introduced, it probably is not widely used, and it might make
sense to obsolete HOR_VISC_ANSWER_DATE and eliminate the older, inconsistent
block of code.  This commit could change answers for answer dates that are above
20241201 with the MEKE Rossby number scaling enabled via MEKE_BACKSCAT_RO_C > 0.
@marshallward
Copy link
Collaborator Author

Oh darn, wrong branch. See you all soon...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.