GLS & MY25 Advection and Constraints. Posted from PECS2024

[harcourt]

Attached are revisions of 4 code modules for the predictor and corrector steps of GLS and MY25 second moment closure options. They correct an incongruity in the default cpp selection of predictor and corrector advection schemes, one that seems to go back at least 2 decades to version 2.0 per archived copies of these files.

Currently, in the default turbulence advection mode (i.e. NOT defining either K_C4ADVECTION or K_C2ADVECTION), the corrector steps under both GLS and MY25 closures default to third-order upstream but the two predictor steps of MY25 and GLS each default to fourth-order centered.

The result of this current (V4.1-V4.2) incongruity in default turbulence advection is that the turbulence length scale is wildly noisy for GLS under strong advection by tides or coastal currents, and the resulting tendency to instability has been found to lock in the high default setting of gls_kmin=7.6E-6 and spuriously high diffusivity in quiescent low-stratification. The problem of spuriously high K's for low N is discussed in posts by M. Scully https://myroms.org/forum/viewtopic.php?p=8433#p8433 and by J. Pringle in https://myroms.org/forum/viewtopic.php?p=18566#p18566.

The modified code modules ...

(1) ... extend the cpp conditional statement in the predictor steps to make K_C4ADVECTION an explicit directive and third order upwind the default, as it is in the corrector modules.

In addition, there are two other changes to the GLS corrector module:
(2) The turbulence length scale is constrained, via a floor on dissipation, to always be less than 0.4 times the total water depth, and ...
(3) ... the potentially noisy error in the Gh limiter I pointed out 6 years ago in viewtopic.php?p=19637#p19637 has also been corrected.

Not sure if (2) and (3) overlap in prophylactic effect or if are really still necessary in the interior domain using the new code with harmonized advection steps, but (2) was instrumental in identifying the problem corrected by (1). That is because it stabilized tidally driven ROMS in LiveOcean (of Parker MacCready & co., Univ. Washington Oceanography) while still producing values of turbulence length scale much larger than water depth. This led to looking closely at the advective schemes, because the corrector step of advection occurs *after* application of (2) limiting turbulence length scale via a floor on dissipation.

The GLS length scale still appears noisy at the domain boundaries and (i think) at river inflow points in default mode with these corrections, though this does not seem to cause LiveOcean to crash. I've been working on a fix for that but have not yet succeeded.

The correction of the cpp if-then-else decision tree is also implemented in the attached for the MY25 choice because although MY25 already constrains the length scale more strongly via the ad-hoc wall function in the dissipation term of the q2l equation, the incongruity of the advection schemes and the negative impact on the GLS closure indicates that MY25 is 'on the ropes' of this ad-hoc length scale constraint probably much more than optimal to convey the skill of its tuning in 1D column models into a 3D model with advection.

An alternative to the changes harmonizing the default selection of turbulence advection scheme is to define either K_C4ADVECTION or K_C2ADVECTION in the current code. Note that these cpp options do not appear to be discussed in current ROMS documentation, and also that doing so might cause advection of turbulence to differ materially from the default third order upstream advection of density/scalars.

KPP has no turbulence scalars to advect, so users selecting that parameterization need not be concerned. Or at least not about this!

Ramsey Harcourt,
With colleagues K. Ravi Prakash & J. B. Mickett at APL/UW Seattle

[harcourt]

Sorry, but those two prestep routines have errors in the new code changes still.

In each of the prestep routines there is a pair of lines, one for each horizontal direction of third-order upstream advection, that says:

IF (cff.gt.0.0_r8) THEN

The first of these should read

IF (XF(i,j).gt.0.0_r8) THEN

and the second should read

IF (EF(i,j).gt.0.0_r8) THEN

These lines are changed in the 4 files in the updated set attached to this post. I'll see if I can delete the erroneous ones above.

However, we are now still testing the new prestep routines in light of the error, pointed out by Guangyu Xu.

Another thing to note is that these changes are only trying to bring the tke and length scale advection routines under gls and my25 options up into consistency with each other. It looks like the Third-order upstream advection of other tracers T,S,etc.. may be using a slightly different formulation now in gls and MY25 corstep routines (which takes the upstream direction from the interpolated velocity, not separately on each as for T,S in pre_step3d.F).

-Ramsey