-
Notifications
You must be signed in to change notification settings - Fork 60
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
Add GL90 parameterization for stacked shallow water #268
Conversation
This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990), Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient nu is computed from kappa_GM via thermal wind balance, and the following relation: nu = kappa_GM * f^2 / N^2. The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free boundary conditions at the top and bottom. In the current implementation, kappa_GM is assumed either (a) constant or as (b) having an EBT structure. A third possible formulation of nu is depth-independent: nu = f^2 * alpha The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth. Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode, in which case we have 1/N^2 = h/g'. More specifically, this commit adds a new subroutine that computes the couping coefficient associated with GL90 via a_cpl_gl90 = nu / h = kappa_GM * f^2 / g' or a_cpl_gl90 = nu / h = f^2 * alpha / h. Further, a_cpl_gl90 is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary layer, whose depth is set by Hbbl_gl90, and 1 otherwise. This modification is necessary to avlid fluxing momentum into vanished layers that ride over steep topography. Finally, a_cpl_gl90 is added to a_cpl, where the latter is the coupling coefficient associated with the remaining vertical stresses, used in the vertical viscosity solver. More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1), Appendix B.
Codecov Report
@@ Coverage Diff @@
## dev/gfdl #268 +/- ##
============================================
- Coverage 37.11% 37.09% -0.03%
============================================
Files 263 263
Lines 73493 73582 +89
Branches 13700 13718 +18
============================================
+ Hits 27279 27296 +17
- Misses 41185 41251 +66
- Partials 5029 5035 +6
📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
@marshallward any ideas why the style / doxygen test fails? I didn't even touch line 90 in |
The log seems to not understand the
I'm not sure why it's getting confused here, where there are other working cases, but I will try to have a look soon. (I'm out of action today and possibly a bit longer, but will try to get to it.) |
Oh, I remember a similar issue happening before. Sometimes a docstring of a previous variable will break, and the error won't manifest until a later line. Have a look at one of the previous lines. |
It's something about the apostrophe No idea why, but it doesn't like that apostrophe. |
Change g' --> g^prime because doxygen and style test does not like the apostrophe.
Nice catch, that fixed it! Thanks. It looks like there is a separate issue with codecov. |
The codecov one is unfortunately server-side (current theory is they limit the number of uploads within a certain timeframe). Resubmitting it usually works, though sometimes you have to wait a bit. |
Ok, we can re-run the github action in a bit? In any case, I'll leave this PR in draft mode until #263 is reviewed and merged because this PR depends on it. |
This variable is analogous to KHTH_USE_EBT_STRUCT, but is specifically for the GL90 scheme. If the user sets KD_GL90_USE_EBT_STRUCT = True, an EBT structure will be applied to KD_GL90. Before this commit, KHTH_USE_EBT_STRUCT took the two-fold role of imposing an EBT structure on KHTH (used in the GM scheme) and KD_GL90 (used in the GL scheme), which may have been confusing.
This PR is ready for review now. @Hallberg-NOAA or @adcroft would be great reviewers as they have seen this code before, while it was being developed. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Overall this excellent PR looks like it is close to being ready to go, but before we merge it into dev/gfdl, I have a few minor specific comments (see below), that I think should be considered. Thank you for this contribution, @NoraLoose!
e253c9b
to
6ff428f
Compare
All get_param calls that are for variables that are only used when GL90 is in use get the argument do_not_log = .not.CS%use_GL90_in_SSW. Prior, these variables have been wrapped by an if-statement.
6da281a
to
2712b20
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am now satisfied that this PR is correct, and I am approving it conditional upon it passing the pipeline testing.
Because of all of the corrections and the merging onto this branch, I am recommending that when this is accepted, it should be done via a squash merge. @NoraLoose would you object to this? |
This PR is undergoing pipeline testing at https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/17774. |
@Hallberg-NOAA sounds good. |
Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/17774 ✔️ 🟡 |
* Fix biharmonic Leith Biharmonic Leith uses Del omega at is-1 and js-1. This unavoidably requires u at js-3 and v at is-3, which are unavailable. It also requires Del omega at Ieq+1 and Jeq+1, which requires v at Ieq+3 and u at Jeq+3, which are unavailable. This necessitates a halo update. Fixes several bugs in Leith+E. - Fixes indexing when computing smoothed vorticity and its gradient - Crucially, computes `vert_vort_mag` when using Leith+E - Fixes some logic in the smoothing code - Other minor indexing fixes * Leith+E Logic Update Ah is required at h and q points. The original code computed Ah at h points, then packed into Ah_h, then applied upper bounds to Ah. If Ah_h is in the diag_table or if debug is true, then the value of Ah with upper bounds get packed into Ah_h. Then, at q points the code unpacks Ah_h. This update makes sure that the upper bound gets applied to q points, not just h points. * Leith+E halo updates The main thing that this commit does is to perform smoothing of u and v outside of the loop over layers. This swaps nz 2D blocking halo updates for a single blocking 3D halo update. * Leith+E smoothing This commit adds a runtime flag, SMOOTH_AH. If True (default) then `m_leithy` and `Ah` are both smoothed, which leads to many blocking communications. If False then these fields are rougher, but there is less communication. * Leith+E eliminate pass-var This commit removes one halo update in Leith+E. To achieve this requires re-indexing two assignments. The value of Ah and Kh are computed at h points, then re-used at q points. Without the halo update it is necessary to offset the assignment at h and q points, e.g. Kh(I,J) = Kh_h(i+1,j+1,k), to avoid accessing values that have not been computed. * Leith+E OBC Adds code so that Leith+E works with OBC. * Leith+E eliminate halo update This commit eliminates one more halo update in Leith+E. * *Correct rotational symmetry with USE_LEITHY This commit revises the smoothing code used when USE_LEITHY = True to give answers that respect rotational symmetry and it also corrects some horizontal indexing bugs and problems with the staggering in some halo update and smooth_x9 calls and reduces some loop ranges to their minimal required values. The specific changes include: 1. Corrected a horizontal indexing bug when interpolating Kh_h and Ah_h to corner (q) points when USE_LEITHY = True. These had previously been inappropriately copied from the thickness point to the southwest of the corner point. This required symmetric-memory-mode calculations of the thickness point viscosities whenever USE_LEITHY is true, but to avoid adding complicated logic, the symmetric-memory loop bounds are used for the calculation of Kh. 2. Revised smooth_x9 to give rotationally symmetric answers and split it into the two routines smooth_x9_h and smooth_x9_uv to reduce the memory used by this routine and reduce the use of optional arguments. 3. Eliminated 4 unneeded halo update calls, and added error handling for the case where Leith options are used with insufficiently wide halos. 4. Added new integers to indicate the loop ranges over which the viscosities and related variables should be calculated, depending on which options are active, and then adjusted 91 do-loop extents horizontal_viscosity code to reflect the loop ranges over which arrays are actually used. 5. Added a new 2-d variable for the squared viscosity for smoothing that can be used for halo updates and to avoid having a variable with confusingly inconsistent dimensions at various points in the code. 6. Corrected the position arguments on 2 smooth_x9 calls and 4 pass_var calls that are used when USE_LEITHY=.true. and SMOOTH_AH=.true. As previously written, these smooth_x9 and pass_var calls would work when in non-symmetric memory mode but would give incorrect answers when in symmetric memory mode. These revisions change answers when USE_LEITHY is true, but answers are bitwise identical in all other cases. --------- Co-authored-by: Robert Hallberg <[email protected]>
GL90 parameterization
This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990), Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient
nu
is computed fromKD_GL90
(corresponding tokappa_GM
) via thermal wind balance, and the following relation:Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode, in which case we have
1/N^2 = h/g'
.The vertical viscosity
del_z ( nu del_z u)
is applied to the momentum equation with stress-free boundary conditions at the top and bottom.GL90 viscosity coefficient: Current options
In the current implementation,
KD_GL90
is assumed either (a) constant or (b) horizontally varying. In both cases, (a) and (b), one can additionally impose an EBT structure in the vertical forKD_GL90
. Another possible formulation of nu is depth-independent:The latter formulation would be equivalent to a
kappa_GM
that varies asN^2
with depth.Figure: Vertical viscosity coefficient
Kv_v
for(upper left) baseline with GL90 parameterization turned off;
(upper right) GL90 parameterization turned on and spatially constant
KD_GL90
;(lower left) GL90 parameterization turned on and
KD_GL90
which is horizontally constant but varies in the vertical with and EBT structure;(lower right) GL90 parameterization turned on with vertically constant GL90 viscosity, achieved by setting
alpha_GL90 = 3 * 10^7 m^2 s
.Note that in the last panel, the diagnostic
Kv_v
is not exactly constant in the vertical becausewhere
h_v
is the layer thickness used at a v-velocity points in the vertical viscosity scheme. The third factor is the composition of a harmonic and arithmetic mean, and is not equal to the identity.Implementation done via coupling coefficients
More specifically, this commit adds a new subroutine that computes the coupling coefficient associated with GL90 via
or
Further,
a_cpl_gl90
is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary layer, whose depth is set byHbbl_gl90
, and 1 otherwise. This modification is necessary to avoid fluxing momentum into vanished layers that ride over steep topography.Finally,
a_cpl_gl90
is added toa_cpl
, where the latter is the coupling coefficient associated with the remaining vertical stresses, used in the vertical viscosity solver.More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1), Appendix B.
Follow-up PR
There will be a follow-up PR which implements diagnostics related to the GL90 parameterization.