Hi.

I've been working on putting different turbulence closures into the COAWST version of ROMS and into GOTM, and I think I have found a longstanding problem in the published account of an aspect of the GLS turbulence closure model, and in the resulting implementation in ROMS.

The problem involves a smoothing of the transition to a limiter on the nondimensional buoyancy forcing Gh in the ROMS/Nonlinear/gls_corstep.F subroutine. In a current or recent version of this routine with svn line

!svn $Id: gls_corstep.F 830 2017-01-24 21:21:11Z arango $

the code with the indicated line numbers reads:

**Code:**

1058 Gh=MIN(Gh,Gh-(Gh-gls_Ghcri)**2/ &

1059 & (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

and this is indeed consistent with the published equations for this aspect of the application of limiters to Gh with the original statement in Eq. (19) of Burchard and Peterson (1999) and its reiteration in the Burchard et al. (1999) GOTM report to the EC. This equation was altered erroneously in Eq. (33) of Warner et al (2005) but this is a typographical error (placing 'Gh - ' into the numerator rather than in front of the fraction). That typo has no direct bearing here on its longstanding implementation in gls_corstep.F, which is faithful to the preceding mathematical statements of how this smoothed limiter should be applied to Gh.

The problem instead is that the original mathematical statement of the smoothed Gh-limiter in Burchard and Peterson (1999) or in Burchard et al. (1999) appear to be incorrect, given the effect it is reported to have in the accompanying Figure 1 of Burchard and Peterson (1999) or Figure 2.3 of Burchard et al (1999). In order to obtain the behavior reported in those figures -- and which makes much more sense than the currently coded limiter -- the lines of code in gls_corstep.F need to be changed to read

**Code:**

1058 Gh=MIN(Gh,Gh-MAX(0.0_r8,Gh-gls_Ghcri)**2/ &

1059 & (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

and the Equation (33) in Warner et al (2005) should read something like:

Gh = Gh_{unlimit} - MAX(0,Gh_{unlimit} - Gh_{crit})^2 / (Gh_{unlimit} + Gh0 + 2 Gh_{crit})

or

Gh = Gh_{unlimit} - (Gh_{unlimit} - Gh_{crit})^2 / (Gh_{unlimit} + Gh0 + 2 Gh_{crit}), for Gh_{unlimit} > Gh_{crit}

Because Burchard et al (1999) and Burchard and Peterson (1999) apply this limiter to a differently-defined nondimensional buoyancy forcing \alpha_N with the opposite sign, the equations (2.90 & 19, resp.) in those publications should have read something like

\alpha_N = MAX{ \tilde{\alpha}_N , \tilde{\alpha}_N - MIN(0,\tilde{\alpha}_N - \alpha_c)^2/(\tilde{\alpha}_N + \alpha_{min} - 2 \alpha_c)

or as is but appended:

\alpha_N = MAX{ \tilde{\alpha}_N , \tilde{\alpha}_N - (\tilde{\alpha}_N - \alpha_c)^2/(\tilde{\alpha}_N + \alpha_{min} - 2 \alpha_c), for \tilde{\alpha}_N < \alpha_c

.

That this is indeed the case can be verified by plotting the different resulting behaviors of Gh or \alpha_N and comparing it to those figures.

The effect of this model error is that when unstable buoyancy forcing is very strong, vertical diffusivity gets unintentionally shut down over the range of Gh approaching the critical value Ghcri, before rising back up to unlimited Gh at Ghcri and transitioning smoothly beyond that as intended, to a maximum value. The smoother is working right for Gh>Ghcri, but puts a notch in Gh down to zero just below Gh=Ghcri, and that notch is wider for some gls closure options than others.

In evaluating the impact of this, bear in mind this is also the regime where local single-point closures are expected to perform poorly due to their typical lack of a nonlocal buoyancy flux. Having the vertical diffusivity turn off when even the correctly limited diffusivity underpredicts buoyancy flux exacerbates this poor performance. Dysfunctionality due to the limiter error has probably been attributed to this fundamental model feature, so that may explain why the erroneous limiter has stuck around so long.

Other implementations such as current/recent versions of GOTM appear to have unceremoniously abandoned this limiter smoothing in favor of a simpler hard upper limit on Gh -- I could not find corresponding code in GOTM. A current or recent implementation in NEMO does show vestiges of an earlier presence in the definition of an unused parameter rghcri in their NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 corresponding to Ghcri.

If anyone is aware of other models where this issue of the incorrectly determined smooth Gh limiter has propagated, please add a comment or direct attention here.

Ramsey Harcourt

University of Washington, APL

harcourt@uw.eduReferences

Burchard, H., Bolding, K., Villarreal, M.R., 1999. GOTM––a general ocean turbulence model. Theory, applications

and test cases. Technical Report EUR 18745 EN, European Commission.

Burchard, H., Petersen, O., 1999. Models of turbulence in the marine environment. A comparative study of twoequation

turbulence models. J. Marine Syst. 21, 29–53.

Warner, J. C., C. R. Sherwood, H. G. Arango and R. P. Signell, 2005: Performance of four turbulence closure models implemented using a generic length scale method, Ocean Modelling, 8, 81-113.