## ROMS KINETIC_ENRG is half volume-average KE

Bug reports, work arounds and fixes

Moderators: arango, robertson

Message
Author
Posts: 521
Joined: Tue Jul 01, 2003 4:12 am
Location: NIWA

### ROMS KINETIC_ENRG is half volume-average KE

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
Posts: 1114
Joined: Wed Feb 26, 2003 4:41 pm
Location: IMCS, Rutgers University
Contact:
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.
Last edited by arango on Mon Apr 10, 2006 5:29 pm, edited 1 time in total.

shchepet
Posts: 185
Joined: Fri Nov 14, 2003 4:57 pm
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