EAMxx Vertical Coordinate

EAMxx physics are being designed to work with a very generic vertical coordinate system and we want to remove as many assumptions as possible and encapsulate the key operations in module subroutines.

We currently assume a finite difference / finite volume type coordinate, with most variables carried at cell centers (and considered to be cell averages through the additional assumption that quantities change linearly with height between their lower and upper boundaries), and a few variables given at cell interfaces ( interface level height z_int and vertical velocity w_int). The layer midpoint positions are taken as physical midpoints with z_mid(k)=(z_int(k+1)+z_int(k))/2.

In the dycore, we are using an “eta” coordinate. HOMME uses the common convention from Arakawa and Lamb (1977) where eta=A+B and pi=A(k)*p0 + B(k)*pi_surf, resulting in eta=1 at the surface and eta close to zero at the model top. To keep the physics as general as possible, we should, following EAM, avoid any knowledge of this coordinate system and work in either pressure or z coordinates.

Open

Interface variables

The height of all layers z_mid and z_int, as well as T_mid, p_mid and pseudo_density (mid), and surface pressure and ptop are given at the start of each physics step. The physics currently applies the thermodynamics at constant pressure, meaning the pressure is held fixed and the height (z_mid and z_int) has to be recomputed after each parameterization, if that parameterization changed the temperature.

p_int is a special case: this should probably also be provided by the dycore. In HOMME, when running in hydrostatic mode (or calling the physics with hydrostatic pressure), p_int is determined by the coordinate system. With nonhydrostatic pressure, p_int needs to be interpolated from surface pressure and p_mid. In HOMME nonhydrostatic mode, there is also a nontrivial calculation of surface pressure described here: https://acme-climate.atlassian.net/wiki/spaces/NGDNA/pages/2579530246

Because of historical use of a hydrostatic mass coordinate, instead of working with density “rho”, we use pseudo_density=g*rho*dz. But there are no (or very few) places in the code where a mass coordinate or hydrostatic approximation is assumed, and we should be able to remove these with consistent use of pseudo_density. One additional assumption is ptop=constant. This assumption could probably be removed in the physics, but it is an assumption made in the theta-l dycore. In the future, we should also consider recomputing the pressure to account for phase changes and water fluxes after each parameterization instead of at the very end of the physics step.

Vertical coordinate requirements:

Because the mid point values represent cell averages, averaging from interface quantities to mid point quantities is naturally defined as:

pmid(k) = (pint(k+1)+pint(k))/2

In addition, derivatives with respect to z of interface quantities are naturally defined by:

dz (mid points) = (zint(k+1) - zint(k))

Integrals of mid point quantities are naturally defined as the anti-derivative of the above equation. Meaning that if we know dz at mid points, and at the surface, we can sum dz to determine zint(k).

To be consistent with the above, mass is defined as the vertical sum of pseudo_density. mass of a quantity given by its specific humidity (referred to here as wet mixing ratio q) is q=sum(q*pseudo_density).

EAM’s geopotential_t() routine:

A key operation tied to this discretization is using the EOS to compute dz from p, pseudo_density and T_v:

p=rho*R*T_v = pseudo_density*R*T_v/(g*dz)

Which is then integrated to get z_int.

This form is generic and applies to any vertical coordinate system with the requirements above, and any equation of state which can define a virtual temperature.

Local Energy

local energy: c_p*T - R*T_v + PHI

Column Energy

Column energy includes the potential energy in the layer above the model top. Total column energy is the correct conserved quantity and takes into account the pressure work required to change the position of the model top:

Column energy: sum_k ( c_p*T - R*T_v + PHI ) + pi_top PHI_top

For processes which have no mass flux, it is advantages to use integration by parts to write the column energy in these forms (where the boundary term remains constant):

Column energy: sum_k ( c_p*T + (phydro/p-1 ) R*T_v ) + pi_surf PHI_surf (where phydro = hydrostatic pressure)

Column energy (Hydrostatic, used by EAM): sum_k ( c_p*T ) + pi_surf PHI_surf

Computing Interface Quantities

zint: Naturally defined by integrating dz at midpoints. zmid computed from zint via simple averaging.

T_int: RRTMG currently computes this using linear interpolation in pressure space. This is generic for any vertical coordinate.

pint: Currently we use pint() mostly as a way to carry around surface pressure (pint(nlev+1)) and ptop (pint(1)), and hence it would be good to populate the interior values as well. RRTMGP uses pint(), but only for interpolation weights in order to compute T_int. HOMME currently supports a get_nonhydro_pressure() routine for midpoint pressure and this interface needs to be extended to interface values. The best algorithm is to use:

This has the nice property that it inverts the INT->MID formula when applied to hydrostatic pressure, and thus close to an inverse of this relationship for NH pressure. (see https://acme-climate.atlassian.net/wiki/spaces/NGDNA/pages/1927512393 )

We considered several other algorithms which are less desirable (left here for reference):

1. Invert pmid(k)=(pint(k+1)+pint(k))/2

   1. This operation can be inverted if we know pint(nlev+1) or pint(1). But it is unstable and can result in a large 2dx oscillation.

2. Linear interpolation in physical space

   1. pint(k)=dz(k-1)*pmid(k) + dz(k)*pmid(k-1) ) / ( dz(k-1) + dz(k))

   2. Other variants: interpolation of log(p) in physical space, interpolate in pressure space using p_int and p_mid

3. Linear interpolate in eta space:

   1. pint(k)=(deta(k-1)*pmid(k) + deta(k)*pmid(k-1) ) / ( deta(k-1) + deta(k))

   2. This algorithm should be avoided in the physics ( no reason for physics to assume an “eta” coordinate).

4. Conservative interpolation

   1. A conservative interpolation method for a wet mixing ratio quantity “q”, would be one that preserves

      1. mass = sum_k [ pseudo_density(k) * (qint(k+1)+qint(k))/2 ]