Possible Coding Error in gls_corstep.F: application of gls_Kmin and gls_Pmin to control advection

[harcourt]

Not 100% sure, but I think there is a coding error in the application of gls_Kmin and gls_Pmin to limit values produced by horizontal and vertical advection of tke and gls for models using any of the GLS-wrapped versions of second moment closure (and hence entering ROMS/Nonlinear/gls_prestep.F & ROMS/Nonlinear/gls_corstep.F).

A note down in the middle of the file mentions that...

! At entry, it is assumed that the turbulent kinetic energy fields
! "tke" and "gls", at time level "nnew", are set to its values at
! time level "nstp" multiplied by the grid box thicknesses Hz
! (from old time step and at W-points).

... and if you look in gls_prestep.F this does indeed appear to be the last thing done to gls and tke at 4th index nnew in gls_prestep.F...

cff=0.5_r8*(Hz(i,j,k)+Hz(i,j,k+1))
...

tke(i,j,k,nnew)=cff*tke(i,j,k,nstp)
gls(i,j,k,nnew)=cff*gls(i,j,k,nstp)

... but when advection is applied in the horizontal immediately thereafter upon entering gls_corstep.F, the limits of gls_Kmin and gls_Pmin are therefore applied to the product of (interpolated) Hz and tke; gls by ...

tke(i,j,k,nnew)=MAX(tke(i,j,k,nnew),gls_Kmin(ng))
...
gls(i,j,k,nnew)=MAX(gls(i,j,k,nnew),gls_Pmin(ng))

... in three different loops over horizontal indices, (once) after adding vertical and (twice) for horizontal advection terms using FXK, FEK,FCK, FXP, FEP, FCP (and these Fxx's do appear to have the application of a factor of thickness from interpolated Huon, Hvom & Hz).

So this means that any benefit intended to model stability by applying kmin to oscillations produced by these advection terms is largely disengaged when thickness Hz, Hvom, Huon is larger than O(1) and is applied using a lower limit when it is smaller than O(1).

I attached a possible fix alt_gls_corstep.F relative to a recent distributed ROMS version.

Looking for feedback if this change looks incorrect to someone with more experience working under the hood with this code.

Likely not a critical repair impacting results, but may conceivably add to other complications associated with and numerical stability when trying to reduce gls_Kmin from its default excessively high setting (see axe grinding, PS below).

Ramsey Harcourt
APL/UW Seattle

P.S. I noticed this possible problem while looking for why vertical mixing was excessively high in deep ocean locations with low stratification, a problem pointed out to me by Guangyu Xu at UW/APL, but the source of high levels of AKs, AKt, AKv was not obvious since the high diffusivity was far from uniform as would reflect a background setting choice. After some spelunking through gls_corstep, I eventually (re-)identified the problem and immediately found my way to Malcolm Scully's 2011 identification of the problem with the default input file setting of gls_Kmin=7.6e-6 m^2/s resulting in a stealth minimum on diffusivity and viscosity scaling as gls_Kmin/N at https://myroms.org/forum/viewtopic.php?p=8433#p8433 , and a more recent 2018 re-identification and exploration of this pitfall by Jamie Pringle at https://myroms.org/forum/viewtopic.php?p=18566#p18566. I found this high-diffusivity/viscosity defect of default gls_Kmin settings particularly hard to find on the board if you don't know you are actually looking for a problem rooted in the 'gls_Kmin' setting, but easy to pull up on the board once you do. So besides the possible coding error above that I noted along the way, this is here is my plug added to other voices that the default suggested gls_Kmin=7.6e-6 setting for GLS closures should be either lowered or marked clearly on example input files or documentation as potentially leading to excessively high deep ocean vertical mixing. Like for N~1e-3/s you get AKt, AKs, Akv ~ 1e-3 m^2/s !!!

[jcwarner]

ok thanks for digging in that deep. it helps a lot.
let me look into it. I do remember those min values were only suggestions, from the literature. But if there is an issue with the min being applied to tke*Hz and not to tke (and gls), then that needs to be resolved.
-j

[sflin]

Hi Harcourt,

Thank you for posting your finding and sharing the code. We recently also noticed that vertical mixing is excessively high in the deep ocean. The high values of Akv closely follow the bathymetry.
If I did not take it wrong, in your revised code, gls_Kmin and gls_Pmin were multiplied by the grid-box height (cff1=0.5_r8*(Hz(i,j,k)+Hz(i,j,k+1))).
The value of cff1 is usually much larger than 1 in deep water. That means tke would possibly become larger in the revised code, and then lead to a larger Akv. How could that address the issue of too strong mixing?


Thanks,
sflin



cff1*gls_Kmin(ng)

Code: Select all

!  Time-step horizontal advection.
!
        DO j=Jstr,Jend
          DO i=Istr,Iend
!RRHvvv
!see "At entry,..." comment below and execution of it in gls_prestep.F
!If minimum needs to be applied it needs a factor of Hz, which could be O(100)m
            cff1=0.5_r8*(Hz(i,j,k)+Hz(i,j,k+1))
!RRH^^^
            cff=dt(ng)*pm(i,j)*pn(i,j)
            tke(i,j,k,nnew)=tke(i,j,k,nnew)-                            &
     &                      cff*(FXK(i+1,j)-FXK(i,j)+                   &
     &                           FEK(i,j+1)-FEK(i,j))
!RRHvvv
!           tke(i,j,k,nnew)=MAX(tke(i,j,k,nnew),gls_Kmin(ng))
            tke(i,j,k,nnew)=MAX(tke(i,j,k,nnew),cff1*gls_Kmin(ng))
!RRH^^^
            gls(i,j,k,nnew)=gls(i,j,k,nnew)-                            &
     &                      cff*(FXP(i+1,j)-FXP(i,j)+                   &
     &                           FEP(i,j+1)-FEP(i,j))
!RRHvvv
!           gls(i,j,k,nnew)=MAX(gls(i,j,k,nnew),gls_Pmin(ng))
            gls(i,j,k,nnew)=MAX(gls(i,j,k,nnew),cff1*gls_Pmin(ng))
!RRH^^^
          END DO
        END DO

[arango]

I haven't had the time to check this issue. We need to analyze the caping of the values. It is probably in the wrong place and needs to be postponed after all the terms are computed. ROMS equations and conservation statements are written in flux form, so things change between horizontal and vertical time-stepping. If you follow the units, they are terms of transport. Please, wait until we look into it.

[xmuliuy]

Hello, it appears there may be a small coding error in the gls_corstep.F file. I recently read the paper by Warner published in Ocean Modelling in 2005, the code provides more details. However, I noticed a discrepancy in line 1011 where the computation of Gh differs slightly from the formulation presented in Equation (33) of the paper.
The formulation for Equation (33) is provided as a figure attached below.
the origin code:

Code: Select all

            Gh=MIN(Gh,Gh-(Gh-gls_Ghcri)**2/                             &
     &                    (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

In the code, it is necessary to add a parenthesis in the numerator, the code may be change to:

Code: Select all

            Gh=MIN(Gh,(Gh-(Gh-gls_Ghcri)**2)/                             &
     &                    (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

Warner, J. C., et al. (2005). "Performance of four turbulence closure models implemented using a generic length scale method." Ocean Modelling 8(1): 81-113.