bug in calculating the Richardson number in zdfkpp.F90

[pmyers]

Dear NEMO team: I was looking at how the Richardson number is computed, and noticed it is done differently in two routines, zdftke.F90 and zdfkpp.F90. After examining the two, I think there may be a bug in the version of zdfkpp.F90 - In zdfkpp, the averaging factor of 0.5 is multiplied in after the shear is squared, and so I don't believe in gets correctly squared. Additionally, when the actual Richardson number is computed, epsilon is not added to the square of the shear term, so it would be possible to get division by zero errors when the shear is zero (first time step when starting model from rest, for exampe). I also note that the two computations are based on different velocities (ub,vb) and (un,vn).

Paul

[gm]

Replying to pmyers@ualberta.ca:

Dear Paul:

I discover with you that the Richardson number is compyted differently in zdfkpp and zdftke. The former choose to averaged the square of the velocity shear, while the latter compute the square of the averaged shear. In other word the shear is larger in KPP than in TKE. I don't know for KPP, but for TKE, the choice was first made a long time ago for zdfric (Pacanoski & Philander formulation) and just apply in TKE. I completely agree on the error in zdfkpp (there is an additional factor of 0.5 which is wrong).

For the epsilon, on line 352 of zdfkpp and on line 369 of zdftke you have a 1.e-20 added so non division by zero (N.B. if not on all land points you must have a floating point exception!).

Last point: the time used on velocity. KPP is a diagnostic computation, un or ub can be used (I guess). TKE is a prognostic equation using a forward time stepping: ub have to be used ! Note that with energetic consideration it can be shown that the used of dz(ub)*dz(un) instead of dz( ub2) is better (in this case, the KE dissipated by vertical momentum mixing is used to create TKE). This should be change in forthcoming issue.

Gurvan

Dear NEMO team: I was looking at how the Richardson number is computed, and noticed it is done differently in two routines, zdftke.F90 and zdfkpp.F90. After examining the two, I think there may be a bug in the version of zdfkpp.F90 - In zdfkpp, the averaging factor of 0.5 is multiplied in after the shear is squared, and so I don't believe in gets correctly squared. Additionally, when the actual Richardson number is computed, epsilon is not added to the square of the shear term, so it would be possible to get division by zero errors when the shear is zero (first time step when starting model from rest, for exampe). I also note that the two computations are based on different velocities (ub,vb) and (un,vn).

Paul

[rblod]

The additional factor 0.5 has been removed.

Thanks a lot for your comments

Rachid