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 2022-10-27 #1586

Merged

Conversation

marshallward
Copy link
Collaborator

This pull request includes several major new features, along with the usual
batch of bugfixes and code improvements.

Major features include:

  • Support for tidal self-attraction and loading

  • Layer interface filtering

  • Improved support for wind stress in mixed layers, especially with respect to
    vertical viscosity.

  • Control of channel drag thickness layers

Along with these features are many new configuration parameters.

Bug fixes include the following: tracers in the ALE sponge, ice shelf
initialization, OBC tracer nudging, predictor timestep restart field handling,
split-timstep diagnostics, and oblique OBCs.

Note to reviewers:

Several parameters have new default values, which can potentially change
answers. These include parameters which may have been previously unconfigured.
Please review the parameters listed below, and compare to your existing
MOM_parameter_doc.all files.

Detailed descriptons for each change are available in the pull request
comments.


Features

Bugfixes

Misc

Contributors:

New Parameter Defaults

  • KAPPA_SHEAR_VERTEX_PSURF_BUG: True -> False
  • BAROTROPIC_TIDAL_SAL_BUG: True -> False
  • LAYER_Z_INIT_IC_EXTRAP_BUG: True -> False

Added to MOM_parameter_doc.all:

  • APPLY_INTERFACE_FILTER
    • INTERFACE_FILTER_TIME
    • INTERFACE_FILTER_MAX_CFL
    • INTERFACE_FILTER_ORDER
    • INTERFACE_FILTER_ISOTROPIC
  • KV_ML_INVZ2
  • TIDAL_SAL_SHT
  • SUBGRID_TOPO_AT_VEL
  • KD_ML_TOT
  • CHANNEL_DRAG_MAX_BBL_THICK
  • USE_POROUS_BARRIER
  • PORBAR_ANSWER_DATE
  • FIXED_DEPTH_LOTW_ML
  • PORBAR_ETA_INTERP
  • VON_KARMAN_CONST
  • PORBAR_MASKING_DEPTH
  • LOTW_VISCOUS_ML_FLOOR
  • KV_EXTRA_BBL
  • VON_KARMAN_BBL
  • STORE_CORIOLIS_ACCEL

Removed from MOM_parameter_doc.all:

  • KDML
  • HENYEY_IGW_BACKGROUND_NEW
  • KVBBL
  • KVML

Added to MOM_parameter_doc.short:

  • KV_ML_INVZ2

Removed from MOM_parameter_doc.short:

  • KAPPA_SHEAR_VERTEX_PSURF_BUG
  • BAROTROPIC_TIDAL_SAL_BUG
  • LAYER_Z_INIT_IC_EXTRAP_BUG
  • KVML

WenhaoChen89 and others added 30 commits August 19, 2022 14:50
* Change dumbbell initialization
* Change in Dumbbell Layer Mode
* Fix sponge diagnostics
* Fix ALE Sponge Diagnostics
* A minor style change removing spaces around = in optional. function arguments.
* Fix ALE sponge diagnostics
* character declaration fix
  - Ice shelf needs to be initialized prior to
    ocean state initialization. This fixes any cases with
    an ice shelf using the FMScap.
  Added the new runtime parameter CHANNEL_DRAG_MAX_BBL_THICK, to allow the
CHANNEL_DRAG parameterization to use the full dynamically determined depth of
the bottom boundary layer, or to use a fixed maximum thickness, regardless of
the settings of other parameters, although for now the default value is set
based on the settings of USE_JACKSON_PARAM and DRAG_AS_BODY_FORCE in order to
reproduce existing solutions by default.  There are new entries in some
MOM_parameter_doc files.  All answers in the MOM6-examples test suite are
bitwise identical, and this test suite has been tested and works with revised
values of this parameter.
* Porous barrier variables, which are grouped in porous_barrier_ptrs,
are changed from pointers to allocatables.
* Copies (aliases) of these variables in MOM_CS are removed.
* The name porous_barrier_ptrs is yet to be changed to minimize the
number of files needed to be edited.
* The interface porous_widths is deleted and subroutine por_widths is
renamed as porous_widths.
* porous_barrier_CS is added to control input and diagnostics of the
 module
* Diagnostics for both interface and layer averages weights are added in
subroutine porous_widths.
* An _init subroutine is added to facilitate reading parameters and
registering diagnostics
* checksum debugs are added within subroutine porous_widths.
This commit primarily fixes indexing bugs in subroutine porous_widths.

* Loop range is subroutine porous_widths is changed from data domain to
computation domain.
* find_eta call is now with a proper halo to cover all velocity points.
* Halo update for porous barrier fields is added in step_MOM_*.

Other changes:
* mask_depth, a component of porous_barrier_CS is now used to as
the criterion for masking. This removes the limit that the cell edge
depth has to be below sea surface.
* Output variable eta_cor was unused and is now changed to a local.
* Some unused indexing variables are removed.
parameters

* subroutine set_subgrid_topo_at_vel_from_file is added to read max, min
and avg depth at the edges from file.
* The subroutine is only called when a new runtime parameter
SUBGRID_TOPO_AT_VEL is True. Default is false.
* id_clock is added (as a private variable) to time porous barrier
calculations.
* Some small format fixes
* A runtime parameter PORBAR_ETA_INTERP is added to control different
methods for calculating interface height at the velocity points from
adjacent tracer cells.
* Two small thickness variables are added and scaled the unit of eta.
This is to assit the harmonic mean calculation, but eventually the
if-test eta_s - eta_prev>0 should be relaced by Angstrom.
* Runtime parameter POROUS_BARRIER_MASKING_DEPTH is renamed to make it
a liitlle bit shorter.
* log_version call is added to separate out porous_barriers module in
parameter file.
The main loops in porous_widths are simpified and cleaned up.
* calc_por_layer iss slightly modified to reduced the if-blocks in
the main loops.
* k-loops are moved out.
* To assist this new structure, two 2-D arrays are added to the stack.
* A new function eta_at_uv is added to treat different interpolations.
This is currently done at every point within the loop, and it is
probably adding too much an overhead. It is better pre-calculate
 eta_[uv] all together, which would increase the stack size.
The averaged weight over a layer (por_face_area[UV]) and the weight
at the interface (por_layer_width[UV]) are now calculated
separately, as the only two updates in step_mom are in fact
using only one of them. This reduces
some unnecessary calculations.

* Two new subroutines replaced the original `porous_widths`. They could
be combined with interface if we remove porous_barrier_ptrs type, and
use variables (of different nz) as output.
* There is a place holder switch do_next in the two calc_por subs. This
is used to further reduce calculations for all layers above D_max.
Will be tested and implemented in the next commit.
* A new subroutine calc_eta_at_uv is added to calculate eta at uv.
This simplifies the code, but also increases stack and some performance
overhead.
* A new runtime parameter USE_POROUS_BARRIER is added to control this
feature. The default is True at the moment to assist potential test. It
should be changed to false in the future.
* A boolean array is added to avoid unnecessary porous barrier
 calculations above the shallowest points.
* The layer-averaged weights are bounded by 1.0 explicitly, to avoid
the cases it goes beyond due to roundoff errors.
* For very thin layers (Angstrom), layer-averaged weights are set to
zeros for simplicity.
  * Perhaps it is more reasonable to use the inferace weight for these
  cases.
  * It might be useful to make the thin layer definition a runtime
  parameter.
* A bug related the size of eta_[uv] is fixed.

This commit is potentially answer-changing when the porouss barrier
module is used.
* The recently introduced ANSWER_DATE functionality is used to maintain a
version that should be bit-for-bit reproducible for previous porous barrier
related runs.
* Rename porous_barrier_ptrs to porous_barrier_type
* Some small documentation edits
* Fixed a bug that eta_[uv] was not properly initialized, which caused issues
with global chksums with DEBUG=True.
* Documents added to comply with Doxygen tests
* Fix nudged OBCs for tracers
Failed codecov.io uploads now report the stderr output.

This will allow us to start diagnosing the intermittent failures and,
possibly, find a solution.
Report error logs from codecov upload fail
  Added the new module MOM_interface_filter to allow for a biharmonic or
other order of filtering of the interface heights.  This new capability
is enabled with the new runtime parameter APPLY_INTERFACE_FILTER, and with
rates of filtering determined by the new runtime parameters
INTERFACE_FILTER_TIME, INTERFACE_FILTER_MAX_CFL, INTERFACE_FILTER_ORDER and
INTERFACE_FILTER_ISOTROPIC.  Set APPLY_INTERFACE_FILTER to True and provide a
positive timescale in INTERFACE_FILTER_TIME to enable the filtering.

  The comments describing THICKNESSDIFFUSE and THICKNESSDIFFUSE_FIRST were
also updated to clarify the distinctions with the new capabilities or to
clearly identify which parameters work together.

  By default, all answers are bitwise identical, but there are new entries or
modified comments in the MOM_parameter_doc files.
  Corrected the doxygen comments for the interface filtering module, to avoid
a section-name conflict with the GM module, and renamed some internal variables
for greater clarity when reading the code.  All answers are bitwise identical.
The sigsetjmp function is part of the POSIX, but is not required to be
defined as a symbol, and may be implemented as a macro.  Since Fortran
C bindings require a symbol, we cannot bind to macro implementations.

The prior implementation assumed a Linux glibc binding of __sigsetjmp
(accessed by a `sigsetjmp` macro), but this did not work on BSD and
MacOS builds, which have a dedicated `sigsetjmp` symbol.

Although the autoconf build included a macro to test and assign the
symbol to `SIGSETJMP_NAME`, this did not resolve builds based on mkmf or
similar build systems, and would fail to compile.

To resolve this, the SIGSETJMP_NAME points to a placeholder function,
`sigsetjmp_missing` which permits compilation but raises an error if
called.

Since this function is only used in our unit testing, and even then only
for tests which would otherwise raise FATAL, this change will not
disrupt any simulations.

However, it does mean that only "power" users who build with either
autoconf or `-DSIGSETJMP_NAME=\"...\"` will be able to run the unit
tests.  In practice, it should be sufficient to direct users to the
autoconf builds, and no actual disruptions are expected.
Disable sigsetjmp for default compilation
  Added a new argument to set_diffusivity_init and bkgnd_mixing_init to specify
when a physically based ocean surface boundary layer mixing scheme, like KPP,
ePBL or a bulk mixed layer is in use, and use this to restrict when a constant
diffusivity is used as a crude parameterization of the surface boundary layer.
This constant mixed layer diffusivity was formerly set with KDML, but is now set
with KD_ML_TOT, but KDML will still work for input but is no longer logged and
triggers a warning message if it is used.  Also removed KDML from the
description of the ADIABATIC option.  This leads to changes in the entries of
several MOM_parameter_doc files, and it could lead to answer changes that
set a non-default value of KDML and use ePBL or KPP, but it seems unlikely that
any such cases would exist.
  Replaced the runtime parameter KVBBL, which sets the total viscosity in the
bottom boundary layer when BOTTOMDRAGLAW is false, with KV_EXTRA_BBL, which sets
the viscosity anomaly relative to KV or other turbulent viscosities, and has a
default value of 0 m2/s, instead of the previous nonzero default for KVBBL.
Specifying KVBBL will still give equivalent results, but is no longer logged and
triggers a warning message if it is used.  This leads to changes in the entries
of several MOM_parameter_doc files, and while all existing solutions are bitwise
identical, by changing the definition of KVBBL there could be some changes at
the level of roundoff for cases that do not use BOTTOMDRAGLAW = True.
  Corrected a spelling error, from HARMOINIC to HARMONIC, in the get_param
description of one of the valid options for PORBAR_ETA_INTERP, and a few other
spelling errors in comments in MOM_porous_barriers.F90.  All answers are bitwise
identical, but an entry in some MOM_parameter_doc files will change.
  Renamed the runtime parameter KVML with KV_ML_INVZ2 and modified the comments
describing it to more clearly indicate what it does and to explicitly discourage
its use. Specifying KVML will still give identical results, but is no longer
logged and triggers a warning message if it is used.  The two descriptions of
the HBBL runtime parameter in MOM_vert_friction and MOM_set_viscosity were also
made consistent.  Also added units to several comments describing real variables
and corrected some spelling errors in comments in MOM_vert_friction.F90.  This
commit leads to changes in the entries of several MOM_parameter_doc files, but
all existing solutions are bitwise identical.
  Refactored the way that find_coupling_coef sets the viscous coupling in the
surface boundary layer, in preparation for adding new options.  All answers are
bitwise identical.
  Eliminated OBC arguments that are no longer used by three internal routines in
MOM_lateral_mixing_coeffs.  All answers are bitwise identical.
@codecov
Copy link

codecov bot commented Oct 28, 2022

Codecov Report

Merging #1586 (84682aa) into main (323ba0c) will decrease coverage by 0.18%.
The diff coverage is 33.92%.

❗ Current head 84682aa differs from pull request most recent head 7c4dfd4. Consider uploading reports for the commit 7c4dfd4 to get more accurate results

@@            Coverage Diff             @@
##             main    #1586      +/-   ##
==========================================
- Coverage   37.36%   37.18%   -0.19%     
==========================================
  Files         261      263       +2     
  Lines       71443    73035    +1592     
  Branches    13365    13608     +243     
==========================================
+ Hits        26698    27161     +463     
- Misses      39805    40859    +1054     
- Partials     4940     5015      +75     
Impacted Files Coverage Δ
src/core/MOM_PressureForce_Montgomery.F90 10.71% <0.00%> (ø)
src/core/MOM_continuity.F90 60.71% <ø> (ø)
src/core/MOM_dynamics_unsplit.F90 79.43% <ø> (ø)
src/core/MOM_dynamics_unsplit_RK2.F90 81.06% <ø> (ø)
src/core/MOM_variables.F90 46.96% <ø> (ø)
src/framework/posix.F90 50.00% <0.00%> (-5.27%) ⬇️
src/initialization/MOM_shared_initialization.F90 29.29% <0.00%> (-1.25%) ⬇️
src/parameterizations/lateral/MOM_MEKE.F90 41.92% <ø> (ø)
...parameterizations/lateral/MOM_interface_filter.F90 0.00% <0.00%> (ø)
...ameterizations/lateral/MOM_spherical_harmonics.F90 0.00% <0.00%> (ø)
... and 54 more

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@marshallward
Copy link
Collaborator Author

The regressions are due to changes in the default parameters, and can be overridden.

Copy link
Collaborator

@kshedstrom kshedstrom left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I approve this PR.

@jiandewang
Copy link
Collaborator

default value for KVML is changed from 1E-4 to 0, set KVML=1E-4 will be the same as KV_ML_INVZ2=1E-4. Just wan to confirm it.

@marshallward
Copy link
Collaborator Author

default value for KVML is changed from 1E-4 to 0, set KVML=1E-4 will be the same as KV_ML_INVZ2=1E-4. Just wan to confirm it.

Yes, that's right. If KV_ML_INVZ2 is unset (or set to a negative number for some reason), then KVML is enabled and used in its place.

I'm not sure how the original default value worked, but it may have been from KV.

Copy link
Collaborator

@abozec abozec left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

COAPS approves the changes

@jiandewang
Copy link
Collaborator

works fine with intel compiler but job crashed with gnu compiler at
https://github.com/NOAA-GFDL/MOM6/blob/dev-gfdl-candidate-main-2022-10-27/src/parameterizations/lateral/MOM_spherical_harmonics.F90#L298
with error: Fortran runtime error: Attempt to DEALLOCATE unallocated 'cs'

However, when compiled with extra DEBUG flags under gnu, job finished fine, see UFS DEBUG flag at
https://github.com/ufs-community/ufs-weather-model/blob/develop/cmake/GNU.cmake#L13

@marshallward I put the run log on GAEA at /lustre/f2/dev/ncep/Jiande.Wang/For-Marshall, see line 160

@marshallward
Copy link
Collaborator Author

marshallward commented Nov 2, 2022

Thanks @jiandewang I've looked at the code and I think it makes sense.

These lines appear in tidal_forcing_init():

call get_param(param_file, mdl, "TIDES", tides, &
"If true, apply tidal momentum forcing.", default=.false.)
if (.not.tides) return

If TIDES is false then the fields in sht_CS are never allocated.

(BTW this check is perhaps redundant, since there is another TIDES check in the RK2 solver. But either way, sht_CS is not initialized if TIDES is false.)

However, there is no similar check around tidal_forcing_end() or spherical_harmonics_end().


There are a few ways to fix it, probably the simplest is to just add if (allocated()) checks to spherical_harmonics_end(). If I were doing it, I would probably add the use_tides logical to either the RK2 CS or the tidal forcing CS and use that to decide whether or not to call spherical_harmonics_end().

Either fix is fine with me. Since we're in the middle of a PR review, maybe the simplest approach of if (allocated()) deallocate() is best.

This code is from @herrwang0 so he may have an opinion on how to fix this.

@herrwang0
Copy link
Contributor

It does look like there is a missing use_tides check in subroutine end_dyn_split_RK2, compared with initialize_dyn_split_RK2

call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
CS%tides_CSp)
call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
ntrunc, CS%vertvisc_CSp)

call hor_visc_end(CS%hor_visc)
call tidal_forcing_end(CS%tides_CSp)
call CoriolisAdv_end(CS%CoriolisAdv)

I agree with Marshall, adding use_tides check to end_dyn_split_RK2 should fix the bug and make it consistent with the init call. Adding if allocated() checks to the variables in spherical_harmonics_end()is also good to me (and again it is consistent with tidal_forcing_end()). So either (or both) is fine with me.


Side points:

  1. Related to this, I find that there are probably a bunch of end calls missing in the unsplit dynamical cores, say e.g.

    subroutine end_dyn_unsplit(CS)
    type(MOM_dyn_unsplit_CS), pointer :: CS !< unsplit dynamics control structure that
    !! will be deallocated in this subroutine.
    DEALLOC_(CS%diffu) ; DEALLOC_(CS%diffv)
    DEALLOC_(CS%CAu) ; DEALLOC_(CS%CAv)
    DEALLOC_(CS%PFu) ; DEALLOC_(CS%PFv)
    deallocate(CS)
    end subroutine end_dyn_unsplit

  2. Also, I wonder what is the general option on the initialized flag is nowadays? I suppose we could use it in the _end() routines to be more defensive and able to expose bugs like this (when end routines are erroneously called).

If (.not. CS%initialized) call MOM_error(FATAL, ...)

@marshallward
Copy link
Collaborator Author

@jiandewang I have a PR which hopefully fixes this issue. NOAA-GFDL#230
@herrwang0 if you want to test or review that would be helpful.

As for the questions above:

  • I think there are still significant gaps in the "finalize" usage. I did a valgrind check some time ago and caught many of them, but no doubt there are many more unresolved. I would also like to see these handled by final procedure language feature, but unfortunately it's not well-supported at the moment.

  • I don't personally support the initialized flag checks, and would rather this state be determined this the configuration logic of the run, but I think there's still quite a few around being used for important checks. Some have replaced older associated() checks from when these were pointers.

@jiandewang
Copy link
Collaborator

@jiandewang I have a PR which hopefully fixes this issue. NOAA-GFDL#230 @herrwang0 if you want to test or review that would be helpful.

As for the questions above:

  • I think there are still significant gaps in the "finalize" usage. I did a valgrind check some time ago and caught many of them, but no doubt there are many more unresolved. I would also like to see these handled by final procedure language feature, but unfortunately it's not well-supported at the moment.
  • I don't personally support the initialized flag checks, and would rather this state be determined this the configuration logic of the run, but I think there's still quite a few around being used for important checks. Some have replaced older associated() checks from when these were pointers.

@marshallward the fixing works fine in UFS

@jiandewang jiandewang self-requested a review November 7, 2022 18:02
This patch moves the local `use_tides` flag of the split RK2 solver,
which tracks the TIDES parameter, into the RK2 control structure (CS),
and uses it to conditionally call `tidal_forcing_end().

The `tidal_forcing_init()` function conditionally allocates the
spherical harmonic fields, but `tidal_forcing_end()` is always called,
which will in turn always attempt to deallocate the fields.  In some
compilers, under certain levels of debugging, this can raise a runtime
error.
Copy link
Collaborator

@gustavo-marques gustavo-marques left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

NCAR approves this PR.

@marshallward
Copy link
Collaborator Author

This looks ready to merge. Thank you as always for the verification.

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.