Preprocessing¶

This page explains what the preprocessing step does, what it expects as input, and what it produces.

The function you call is:

from flexprep.preprocessing import preprocess

dm_out = preprocess(dm_in)

What it does¶

Precipitation de-accumulation¶

Inputs	Units (in)	Output	Units (out)	Requirements / Notes
`cp`, `lsp`	m (accumulated from step 0)	`cp`, `lsp` (rates)	mm h⁻¹	Need ≥ 2 steps; first valid at step 1.

Convective (cp) and large-scale (lsp) precipitation are accumulated totals (in meters) from the start of the forecast. We convert them to an intensity (in mm h⁻¹) by differencing consecutive steps and normalizing by the elapsed time:

\[\text{rate}_i = \frac{P_i - P_{i-1}}{t_i - t_{i-1}} \times 1\ \text{hour}\]

t_i - t_{i-1} comes from the step coordinate.
Multiply by 1000 to convert m/h → mm/h.
The first valid value is at step 1, so step 0 is dropped.

Radiation / Heat Flux / Surface Stress de-accumulation¶

Inputs	Units (in)	Output (rate)	Units (out)	Requirements / Notes
`ssr` — surface net shortwave radiation	J m⁻² (accumulated)	`ssr`	W m⁻²	Need ≥ 2 steps
`sshf` — surface sensible heat flux	J m⁻² (accumulated)	`sshf`	W m⁻²	—
`ewss` — eastward turbulent surface stress	N·s m⁻² (accumulated)	`ewss`	N m⁻²	—
`nsss` — northward turbulent surface stress	N·s m⁻² (accumulated)	`nsss`	N m⁻²	—

These variables are provided as accumulated totals over time. During preprocessing, we turn them into instantaneous per-second rates by differencing consecutive forecast steps and dividing by the elapsed time.

ssr: J m⁻² → W m⁻²
sshf: J m⁻² → W m⁻²
ewss, nsss: N·s m⁻² → N m⁻²

Requires ≥ 2 steps; the first valid rate is at step 1.

Omega (vertical pressure velocity) — concept & implementation¶

Inputs (required)	Units	Output	Units (out)	Requirements / Notes
`sp` — surface pressure	Pa	`omega`	Pa s⁻¹	Uses hybrid coeffs; slice levels 40..136; drop step 0 if present
`etadot` — vertical velocity in η	s⁻¹
`ak`, `bk` — hybrid coefficients	—			Derived from GRIB `pv`; consistent grid required

This section explains what the vertical omega calculation does, why it is needed, and how the discrete operator works in the preprocessing pipeline.

Target variable: omega (Pa s⁻¹)
Source variable: etadot (s⁻¹)

1) Why we need `dp/dη`¶

IFS provides vertical velocity in η-coordinates:

etadot = dη/dt

FLEXPART requires pressure vertical velocity:

omega = dp/dt

Using the chain rule:

omega = (dp/dη) * etadot

2) Hybrid vertical coordinate (η)¶

ECMWF IFS uses a hybrid sigma–pressure coordinate:

p(η) = A(η) + B(η) * ps

A(η) = ak, B(η) = bk
ps is surface pressure (Pa)
Near surface: B ≈ 1
Aloft: B → 0

For two adjacent half-levels:

Δp_k = ΔA_k + ps * ΔB_k

3) Discrete `dp/dη`¶

Relative to a reference pressure pref = 101325 Pa:

scale = (ΔA + ps*ΔB) / (ΔA + pref*ΔB)

In code:

ps * (dak/ps + dbk) / (dak/pref + dbk)

4) Layer → level mapping¶

Using a centered scheme:

omega_k - omega_{k-1} ≈ 2 * F_k

This explains the 2.0 factor.

5) Vertical accumulation¶

The omega profile is reconstructed by vertical accumulation:

given: omega_k - omega_{k-1} ≈ 2 * F_k
solve: omega_k by accumulating over k

6) Putting it together¶

pref = 101325.0
dak = ak.diff("level")
dbk = bk.diff("level")

omega_layer = 2.0 * etadot * ps * (dak/ps + dbk) / (dak/pref + dbk)
omega = omega_layer.reduce(cumdiff, dim="level")

Inputs, outputs, assumptions¶

Inputs - sp — surface pressure (Pa) - etadot — vertical velocity in η (s⁻¹) - ak, bk — hybrid coefficients

Output - omega — pressure vertical velocity (Pa s⁻¹)

Assumptions - All inputs share grid and level indexing - Missing ak/bk → omega cannot be computed