Diagnostic expressions

EAMxx builds derived diagnostics from expressions:

field_names:
  - q_flux := (qc+qv)*p_mid
  - cold_T  := T_mid.where(T_mid<273.15).mean('lev')

Expressions are parsed by dexpr, the expression parser shared across E3SM components.

Expressions must be given a name

An expression is not a usable NetCDF variable name, and field_names entries are written to file verbatim. So any expression you want written must be given a name with :=:

field_names:
  - q_flux := (qc+qv)*p_mid    # good
  - (qc+qv)*p_mid              # error: "An expression must be given an output name"

The expression itself is recorded in the variable's alias_of attribute, so the provenance stays in the file.

This applies only to what is written. Expressions nested inside other expressions, and expressions in the aliases: section (intermediates that other diagnostics use but that are not written), need no name of their own.

Operators

Operator	Meaning
`a + b`, `a - b`, `a * b`, `a / b`	element-wise arithmetic
`-a`	negation, shorthand for `(-1)*a`
`>`, `>=`, `<`, `<=`, `==`, `!=`	comparison, only inside `where(..)`
`( )`	grouping

Operands may be fields, other expressions, named physical constants (Rgas, P0, ...), or numeric literals. Constant names are not case sensitive. An operand that names both a constant and a field (including a name introduced with :=) is rejected. Earth's rotation rate is omega_earth; omega is the vertical pressure velocity field.

field_names:
  - p_scaled   := p_mid/100
  - T_scaled   := T_mid*0.5
  - dry_air    := p_mid - qv*p_mid
  - p_ratio    := p_mid / P0

A numeric literal has no units. Multiplying or dividing by one keeps the units of the field, so p_mid/100 holds the pressure in hPa, but is written with the units of p_mid, which are Pa; nothing converts units for you. Adding a number to, or subtracting it from, a field with units is rejected when the stream is set up, since + and - need operands with the same units: T_mid - 273.15 is an error.

Precedence is the usual one, so qc+qv*p_mid is qc+(qv*p_mid). When in doubt, parenthesize. (That example only illustrates grouping: qc+(qv*p_mid) adds a mixing ratio to a pressure, and is rejected at run time for incompatible units, while (qc+qv)*p_mid is fine.)

Functions

Functions are written as methods on the field they act on, X.f(..), and are named after their xarray counterparts wherever there is an honest analogue.

Expression	Meaning
`X.isel(lev=10)`	`X` at level 10
`X.isel(lev=0)`	`X` at the model top
`X.isel(lev=-1)`	`X` at the lowest level
`X.interp(p_mid=500, units='hPa')`	`X` interpolated to 500 hPa
`X.interp(z_mid=10, reference='surface')`	`X` interpolated to 10 m above the surface
`X.interp(z_mid=1000, reference='sealevel')`	`X` interpolated to 1000 m above sea level
`X.mean('col')`	horizontal (area-weighted) mean
`X.mean('lev')`	vertical mean
`X.sum('lev')`	vertical sum
`X.mean('lev', weights='dp')`	vertical mean weighted by `dp`
`X.where(qv>0.01)`	`X` where `qv>0.01`, missing elsewhere
`X.derivative('p_mid')`	vertical derivative with respect to `p_mid`
`X.histogram(bins=[0,1,2])`	histogram of `X` over the given bin edges
`X.zonal_mean(bins=20)`	zonal mean in 20 latitude bins
`X.shift(time=1)`	`X` at the previous step
`X.over_dt()`	`X` divided by the timestep
`X.tend()`	tendency of `X` over the last step
`X.broadcast_like(Y)`	`X` broadcast to the layout of `Y`
`p_mid.searchsorted(500, units='hPa')`	index of the first level with pressure above 500 hPa
`z_mid.searchsorted(1000)`	index of the first level below 1000 m above sea level
`height_mid.searchsorted(100)`	index of the first level below 100 m above the surface

Notes:

isel indexes levels the way xarray does: lev=0 is the model top and lev=-1 the bottom. Other negative indices are rejected, since they cannot be resolved without knowing the level count.
interp takes exactly one of p_mid or z_mid, named after the coordinate as in xarray. units defaults to Pa and m respectively; reference ('surface' or 'sealevel', default 'sealevel') applies to z_mid only.
A target outside a column's data range is filled with the missing-value marker, except near the surface with reference='surface'. interp(p_mid=..) fills at either end: there is no meaningful way to extrapolate past the model top or below the surface in pressure space. interp(z_mid=.., reference='sealevel') does the same at both ends, since a target elevation can legitimately sit below a mountain's surface, where the quantity is not defined. But interp(z_mid=.., reference='surface') is expected to always be well formed close to the ground -- wind_speed.interp(z_mid=2,reference='surface') should mean something even when the lowest model level sits well above 2m -- so it extrapolates using the nearest level instead of filling when the target is below the bottom entry. That exception is for the near-surface end only: above the model top, reference='surface' fills too, exactly like the other two cases.
shift looks BACK in time, following xarray's sign convention: a positive shift moves data forward, so shift(time=1) is the value from the previous step. shift(time=-1) would be the value from a step ahead, which a running model cannot know, and is rejected -- do not read it as "one step back" by analogy with isel(lev=-1). Only time=1 is available, since only one step of history is kept.
A shift of a shift is rejected. X.shift(time=1).shift(time=1) would be one step back, not two: a shift copies its operand at the start of each step, and an operand that is itself a shift has already been updated by then. So it is rejected at setup, and so are the same through an alias (l1 := X.shift(time=1), then l1.shift(time=1)) and X.shift(time=1).tend(). For two steps back, put an operation between the shifts: (X.shift(time=1)*1.0).shift(time=1) (on the first step of a run, see the first step). The shift of any other expression is one step back, e.g. bt1.shift(time=1) with bt1 := X - X.shift(time=1).
shift(time=1) of a diagnostic needs a stream that evaluates every step. An instant stream evaluates diagnostics on output steps only, so the value one model step back does not exist there unless the stream writes every step (frequency_units: nsteps with frequency: 1). In any other instant stream, (qv*T_mid).shift(time=1) is rejected at setup, and so is (qv*T_mid).tend(); use an averaged stream. The shift of a model field (T_mid.shift(time=1)) is one model step in any stream.
mean('col') is always area weighted, so weights does not apply to it. weights for 'lev' is 'dp' or 'dz'.
broadcast_like repeats X along the dimensions of Y that X lacks, so that X can be combined with Y. E.g., dp.broadcast_like(horiz_winds) would produce a version of dp (call it dp_bcast) with layout (ncol,2,lev), where dp_bcast(i,:,k)==dp(i,k). The layout of X must be an ordered subset of that of Y. Notice that Y is only used for its layout, its values are never read.
where takes a single comparison. and/or are not supported; chain where(..) calls instead. A chain is the intersection of its conditions: X.where(A>0).where(B>0) is valid only where both hold.
X.sum('lev') of a masked X is masked where no level is valid. Invalid levels are skipped, and a column with none left has no sum; it is written as missing rather than as 0. mean('lev') behaves the same way. This is xarray's sum('lev', min_count=1), not its default sum('lev'), which would give 0.
A scalar threshold is in the field's own units, as the model stores them. T_mid.where(T_mid<273.15) is in kelvin, and p_mid.where(p_mid>85000) in pascals, not hectopascals. A comparison between two fields is checked when the stream is set up: T_mid.where(T_mid>qv) is rejected, since the units differ. (A field that never declared its units is not checked.)
mask and lev are placeholders, not fields. Neither exists in the field manager, and neither means anything on its own. mask is only legal as the receiver of .where(..): mask.where(cond) is the 0/1 indicator of cond itself, rather than some field sampled by it. lev is only legal as the left operand of the comparison inside where(..), where it means the level index. So mask.where(lev>5) is a request, while mask, mask*1.0, mask.mean('lev') and lev>5 are not: each of those leaves mask or lev as an operand to be resolved on its own, and the run aborts with The key 'mask' is not associated to any registered product, since it is neither a field nor a diagnostic. Neither name comes from xarray.
A negative number is written as it reads, as a threshold or as a factor: omega.where(omega<=-0.1), -2*T_mid, T_mid*(-2). In a where condition, the field goes on the left of the comparison.
over_dt and tend cannot go in an instant stream. An instant stream takes a snapshot at t0, before any step has run, and dt is zero there, so the diagnostic aborts with "dt must be positive". Write them to an averaged stream (averaging_type: average, max, min) instead. This is a property of the diagnostic, not of the syntax.
A mask.where(..) result cannot be written to file as is. Its output has int data type and every variable is written as a real, so listing mask.where(..) in field_names is rejected when the stream is set up. Multiplying by 1 promotes it to Real, and mask.where(ps>100000)*1.0 writes fine. Note the *1.0 applies to the field that where produced, which is a perfectly ordinary field; it is not an operation on mask, which per the bullet above cannot be multiplied.
X.tend() is shorthand for (X - X.shift(time=1)).over_dt(), and inherits the over_dt restriction above.
searchsorted gives, per column, the index of the first level whose pressure is above the target pressure (for a height, which decreases with the level index, the first level below the target); the index is nlevs if there is none. That is numpy's searchsorted(side='right') along the level dimension. Its receiver picks the coordinate: p_mid/p_int for a pressure (units defaults to Pa), z_mid/z_int for a height above sea level, and height_mid/height_int for a height above the surface (units defaults to m). interp(..) is built on it, so fields interpolated to the same target share one search. Its output has int data type, so like mask.where(..) it cannot be written as is.
Histogram bin edges are numbers in increasing order. They may be negative, and may be written with an exponent (bins=[-1,0,1e-5,1e-4]).

Named quantities are not operations, and stay plain names in an expression: LiqWaterPath, RelativeHumidity, z_mid, dz, precip_liq_surf_mass_flux, and so on.

Where this departs from xarray

isel, interp, mean, sum, where and shift mean what they mean in xarray, on the same arguments. The rest have no xarray spelling we could adopt:


`histogram`	not in xarray; the keyword follows xhistogram
`zonal_mean`	xarray would be `groupby_bins('lat', 20).mean()`, a call chain we cannot represent
`over_dt`, `tend`	no xarray equivalent
`mean('lev', weights=..)`	xarray would be `weighted(dp).mean('lev')`, again a chain
`searchsorted`	named after numpy's (`side='right'`), but takes `units`, and for heights, which decrease with the level index, it looks for the first level below the target

X.mean('lev', weights='dp') is (X*dp).sum('lev')/dp.sum('lev'), expressed as one diagnostic rather than as a graph of three.

Two operand names have no xarray analogue at all, because they are not fields: mask and lev, both legal only inside where(..). See the note above for what they mean and where they may appear.

Errors

An expression that does not parse, calls a function that does not exist, or uses an operator no diagnostic implements is rejected with a message pointing at the offending part. Unknown function names come back with the list of available ones.

Error! Unsupported expression: no diagnostic implements the operator '**'.
 - subexpression: (T_mid**2)

Worked example

fields:
  physics_pg2:
    aliases:
      # An intermediate: used below, but not written to file. No name needed
      # for the expression itself beyond the alias.
      - warm_T := T_mid.where(T_mid>273.15)

    field_names:
      # Column-mean temperature, but only where it is above freezing
      - warm_T_colavg := warm_T.mean('lev')
      # Surface temperature where the air is moist
      - sfc_moist_T   := T_mid.where(qv>0.0).isel(lev=-1)
      # Moisture flux
      - q_flux        := (qc+qv)*p_mid

Adding a function

The set of callable functions is registered by EAMxx, in components/eamxx/src/share/io/eamxx_dexpr_diags.cpp, not by dexpr itself. Adding one means adding a FunctionSpec to the registry there and a case to the translator that maps it onto a diagnostic; nothing under share/dexpr changes.

Each FunctionSpec carries an example of how the call is written. That is not just documentation: two checks use it, so a function that does not hang together fails fast rather than at the moment a user first writes it.

dexpr::validate_registry, called where the registry is built, proves the example parses, matches the spec that declared it (arity, keyword names, required keywords), and really calls that function. Declare mean as taking no positional arguments while writing T_mid.mean('lev') and this catches it.
The dexpr_every_registered_function_is_buildable case in src/share/io/tests/create_diag.cpp runs every example through create_diagnostic, so a function registered without a translator case is caught too.

So adding a function means: register it with an example, translate it, and both checks come along for free. The dexpr command line tool does the same for the generic vocabulary -- dexpr check "<expr>" and dexpr check-registry.

Those two checks stop at "a diagnostic was built". Every expression on this page is also run end to end -- computed and written to NetCDF -- by the standalone test in components/eamxx/tests/multi-process/dynamics_physics/homme_shoc_cld_p3_rrtmgp_pg2, which carries two extra output streams for the purpose: output_diags_ins.yaml and, for the ones an instant stream cannot hold, output_diags_avg.yaml. Keep them in step with this page: an example added here belongs in one of those.