Inconsistency in definition of energy fluxes

[illorenzo7]

@feathern @BWHindman I am noticing what seems to be a subtle issue in the definition of the energy flux diagnostics, i.e., in

Diagnostics_Thermal_Equation.F90
and
Diagnostics_KE_Flux.F90

There are really two issues.

The first one only occurs when the equations are non-dimensional:

1. In forming the total energy equation to calculate the fluxes (after the equations are non-dimensionalized), the equations need to have the same units. E.g., we can multiply the KE equation (u dot [momentum equation]) by L^2/T^2 (where L and T are the units of length and time in the chosen non-dimensionalization), and multiply the heat equation by the unit of internal energy (e.g., Tbar [S] for anelastic, where [S] is the unit of entropy), and then add them. Or equivalently, we can multiply just the kinetic energy equation by (L^2/T^2) / (Tbar [S]), and then add it to the heat equation.

This latter option is currently what is done in the code, but only for the total enthalpy flux only. See line 608 of Diagnostics_Thermal_Equation.F90 and note c_8 = (L^2/T^2) / (Tbar [S]) (this = Di Pr / Ra if T = [viscous diffusion time]). The other enthalpy fluxes (pp, mm, etc.) are wrong, which I am fixing in pull request #480.

The second issue I believe affects all reference types:

2. The heat equation actually needs to be turned into a "potential energy equation" before being added to the KE equation, to correctly cancel the buoyancy work term. Basically, it needs to be multiplied by the integral of g(r). For anelastic (and dS/dr = 0), integral(g) ~ T_bar and the resultant fluxes (enthalpy flux, heat flux, viscous flux, conductive flux) wind up being correct. But for stable layers or Boussinesq, I don't think the fluxes are right.

Maybe someone has already thought about this and maybe I'm missing something?

Cheers,

Loren