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
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. list-table::
   :header-rows: 1

   * - 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:

.. math::

   \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
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. list-table::
   :header-rows: 1

   * - 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
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. list-table::
   :header-rows: 1

   * - 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
----------------------

.. code-block:: python

   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