ROMS KINETIC_ENRG is half volume-average KE

[m.hadfield]

I'm reporting this in the ROMS Bugs forum, though it may be a feature...

Does anyone other than me find it surprising that the quantity written to stdout at each time step under the heading KINETIC_ENRG is half the volume-average kinetic energy (for a 2D simulation)?

The relevant code is in diag.F. First a column-integrated kinetic energy is calculated for each (i,j). With SOLVE3D undefined the expression is

Code: Select all

            ke2d(i,j)=(zeta(i,j,krhs)+h(i,j))*                          &
     &                0.25_r8*(ubar(i  ,j,krhs)*ubar(i  ,j,krhs)+       &
     &                         ubar(i+1,j,krhs)*ubar(i+1,j,krhs)+       &
     &                         vbar(i,j  ,krhs)*vbar(i,j  ,krhs)+       &
     &                         vbar(i,j+1,krhs)*vbar(i,j+1,krhs))

This quantity ke2d has units m^3/s^2 and is a finite-difference form of 0.5*(h+zeta)*(ubar^2+vbar^2).

This array is summed over the domain in the variable avgke, then this variable is further modified by the expression

Code: Select all

          avgke=0.5_r8*avgke/volume

But we have already applied a factor of 0.5

The other oddity is the fact that for a 3D simulation KINETIC_ENRG is multiplied by density, ie it has units kg m^-1 s^-2. If we compare a 2D simulation with a 3D simulation in which the velocities are vertically uniform, KINETIC_ENRG is about 1025 times as large in the 3D case as in the 2D case.

[arango]

Actually, the problem is not with the half factor which is needed to compute the volume averaged kinetic energy per unit mass (in transport form and at interior RHO-points), but in the way that the velocity components are squared:

Code: Select all

            cff1=0.5_r8*(ubar(i,j,krhs)+ubar(i+1,j,krhs))
            cff2=0.5_r8*(vbar(i,j,krhs)+vbar(i,j+1,krhs))
            ke2d(i,j)=(zeta(i,j,krhs)+h(i,j))*rho0*                     & 
     &                (cff1*cff1+cff2*cff2)

A similar approach is needed for the 3D case. Notice that for consistency, I am using the rho0 factor as you mentioned.

[shchepet]

Actually what is proposed by Hernan is not consistent with Lilly discretization on a C-grid, even thought it may be viewed as approximatrely correct.

Hernan's procedure implies that velocity components are averaged from their U-,V-points to RHO-points first, and then squered there. This is a WRONG way of doing it, because this is NOT the discrete form of kinetic energy which is conserved EXACTLY if time is continuous, advection is discretized by Lilly scheme and Coriolis/curvilinear terms by Arakawa-Lamb, 1977.

The form reported by Mark is better: it is exact in the case of Cartezian coordinates (never mind the factor of 2), but is still wrong in curvilinear. Hernan's form also.

The correct form is:

KE = sum over all i,j of

0.25*(D(i,j)+D(i-1,j))*dn_u(i,j)*dm_u(i,j)*ubar(i,j)**2
+ 0.25*(D(i,j)+D(i,j-1))*dm_v(i,j)*dn_v(i,j)*vbar(i,j)**2

where D=h+zeta is total depth; dn_u is 1/pn computed at u-points, and dm_v is 1/pm computed at v-points.

This and only this form which is consistent with Lilly discretizaton (a second-order version of ROMS; as well as POM). The rule of thumb here is that, KE is a weighted sum of non-averaged squares ubar and vbar, weighted by control volumes of U- and V-points.

To verify it one can simpy multiply discrete momentum equations in conservation form, thus by ubar and vbar respectively and after summation by parts verify that all terms associated with advective fluxes reduce to side boundaries (hence no KE is generated/destroyed inside the domain), while terms associated with Coriolis/curvilinear terms cancel out exactly.

The form proposed by Mark is correct in the case of Cartezian coordinates because it is linear with respect to D=h+zeta, and therefore one can switch order of summation; but in the case of curvilinear grid it is not exact, since pm and pn's are now different.

Note that SCRUM-generation codes always did it the wrong way: averaging velocities before taking square. As the result, inhereted codes may be confusiong.

An obvious consequence is that if one considers a 2dx wave or a checker-board pattern, then Hernan's formula tells that the KE of that motion is zero, while Mark's formula yields a positive number. C-grided models have non-trivial KE associated with 2dx patterns.

[arango]

I agree. The total energy density present in a system is equal to the initial energy plus the energy supplied by the external forces minus the energy lost by mixing process. This kind of conservation energy is not trivial in a model like ROMS or any other terrain-following coordinates model. I always look at these diagnostics as approximated quantities and not as an energy conservation statement. They are useful to determine if the model its is blowing-up and to check the behavior of the total basin volume. They are also useful in spin-up problems.

I was trying to clarify above the 0.25 in the squaring of the velocity components. Now, there are very subtle numerical issues when computing such quantities in a model like ROMS. Recall, that we are computing a volume integral which requires special parallel considerations. This triple integral is prone to randoff. I guess that it can be coded differently to aliviate that. Notice also that in the current version the land sea masking is not included.