Skip to content

Meteorology

Gridded meteorology (MetField) on the les_fvm analysis grid, the wrfout loader, and the adapter that feeds it to the Eulerian solver as wind. Tier 0.

met

Meteorological forcing for the transport tiers (Tier 0 prerequisites).

The transport models are driven by an external meteorology; this package owns the container those models consume and the loaders that populate it from numerical-weather-prediction output.

  • :class:MetField — gridded wind, temperature, pressure and PBL height on a les_fvm analysis grid, on the les_fvm C-grid stagger, with a time axis.
  • :func:load_wrf — read a wrfout file onto an analysis grid.
  • :func:to_prescribed_wind — expose a :class:MetField as the :class:~plumax.les_fvm.wind.PrescribedWindField the Eulerian solver takes, with piecewise-linear interpolation between met times.

MetField

Bases: Module

Meteorological forcing on an analysis grid with a time axis.

Attributes:

Name Type Description
plume_grid PlumeGrid3D

Analysis grid every field lives on.
times Float[Array, ' n_time']

Seconds since t0, strictly increasing, shape (n_time,).
u Float[Array, 'n_time nz ny nx']

Eastward wind [m/s] at U-points (T-point x + dx/2), interior shape (n_time, nz, ny, nx).
v Float[Array, 'n_time nz ny nx']

Northward wind [m/s] at V-points (T-point y + dy/2).
w Float[Array, 'n_time nz ny nx']

Vertical wind [m/s] at T-points.
temperature Float[Array, 'n_time nz ny nx']

Air temperature [K] at T-points.
pressure Float[Array, 'n_time nz ny nx']

Air pressure [Pa] at T-points.
pbl_height Float[Array, 'n_time ny nx']

Boundary-layer height above ground [m] at horizontal T-points, shape (n_time, ny, nx).
t0 str

ISO-8601 timestamp that times counts from (informational; static, so it never enters a trace).

n_time property 

n_time: int

Number of met times.

wind_at 

wind_at(t: Float[Array, '']) -> tuple[Float[Array, 'nz ny nx'], Float[Array, 'nz ny nx'], Float[Array, 'nz ny nx']]

Interior (u, v, w) at time t, piecewise-linear in time.

Met is typically hourly while the transport step is seconds, so the wind is interpolated linearly between the two bracketing met times (policy (b) of the Tier 0 design page). Outside [times[0], times[-1]] the nearest end value is held constant.

Source code in src/plumax/met/field.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def wind_at(
    self, t: Float[Array, ""]
) -> tuple[
    Float[Array, "nz ny nx"],
    Float[Array, "nz ny nx"],
    Float[Array, "nz ny nx"],
]:
    """Interior ``(u, v, w)`` at time ``t``, piecewise-linear in time.

    Met is typically hourly while the transport step is seconds, so the
    wind is interpolated linearly between the two bracketing met times
    (policy (b) of the Tier 0 design page).  Outside ``[times[0],
    times[-1]]`` the nearest end value is held constant.
    """
    return (
        interpolate_in_time(self.times, self.u, t),
        interpolate_in_time(self.times, self.v, t),
        interpolate_in_time(self.times, self.w, t),
    )

to_prescribed_wind 

to_prescribed_wind(field: MetField) -> PrescribedWindField

Expose a :class:MetField as the wind object the Eulerian solver takes.

The returned :class:~plumax.les_fvm.wind.PrescribedWindField samples :meth:MetField.wind_at — so simulate_eulerian_dispersion( wind_field=to_prescribed_wind(field), ...) runs the FV transport on loaded meteorology with no other change. The caller's domain_* triples must describe the same grid field.plume_grid was built on.

Source code in src/plumax/met/field.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def to_prescribed_wind(field: MetField) -> PrescribedWindField:
    """Expose a :class:`MetField` as the wind object the Eulerian solver takes.

    The returned :class:`~plumax.les_fvm.wind.PrescribedWindField` samples
    :meth:`MetField.wind_at` — so ``simulate_eulerian_dispersion(
    wind_field=to_prescribed_wind(field), ...)`` runs the FV transport on
    loaded meteorology with no other change.  The caller's ``domain_*``
    triples must describe the same grid ``field.plume_grid`` was built on.
    """

    def interior_fn(t):
        return field.wind_at(t)

    return PrescribedWindField(plume_grid=field.plume_grid, interior_fn=interior_fn)

load_wrf 

load_wrf(path: str | Path, *, plume_grid: PlumeGrid3D, origin: tuple[float, float] = (0.0, 0.0), times: slice | None = None) -> MetField

Load a wrfout file onto an analysis grid.

Parameters:

Name Type Description Default
path str | Path

NetCDF wrfout file.

 required
plume_grid PlumeGrid3D

Analysis grid (from :func:plumax.les_fvm.make_grid).

 required
origin tuple[float, float]

WRF-frame metres of the analysis grid's (x=0, y=0); see the module docstring.

 (0.0, 0.0)
times slice | None

Optional slice over the file's Time axis.

 None

Returns:

Name Type Description
A MetField

class:MetField on plume_grid.

Raises:

Type Description
ValueError

If any analysis T-, U- or V-point falls outside the WRF mass-point domain (the loader never extrapolates horizontally), or if times selects nothing.
Source code in src/plumax/met/wrf.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def load_wrf(
    path: str | Path,
    *,
    plume_grid: PlumeGrid3D,
    origin: tuple[float, float] = (0.0, 0.0),
    times: slice | None = None,
) -> MetField:
    """Load a ``wrfout`` file onto an analysis grid.

    Args:
        path: NetCDF ``wrfout`` file.
        plume_grid: Analysis grid (from :func:`plumax.les_fvm.make_grid`).
        origin: WRF-frame metres of the analysis grid's ``(x=0, y=0)``;
            see the module docstring.
        times: Optional ``slice`` over the file's ``Time`` axis.

    Returns:
        A :class:`MetField` on ``plume_grid``.

    Raises:
        ValueError: If any analysis T-, U- or V-point falls outside the
            WRF mass-point domain (the loader never extrapolates
            horizontally), or if ``times`` selects nothing.
    """
    # Time decoding off: WRF gives XTIME CF units (``minutes since ...``),
    # which xarray would otherwise turn into datetimes; we want raw minutes.
    with xr.open_dataset(path, decode_times=False) as ds:
        if times is not None:
            ds = ds.isel(Time=times)
        if ds.sizes["Time"] == 0:
            raise ValueError("load_wrf: `times` selects no time steps")
        dx_w = float(ds.attrs["DX"])
        dy_w = float(ds.attrs["DY"])
        u_stag = ds["U"].values
        v_stag = ds["V"].values
        w_stag = ds["W"].values
        theta_pert = ds["T"].values
        pressure = ds["P"].values + ds["PB"].values
        geopotential = ds["PH"].values + ds["PHB"].values
        terrain = ds["HGT"].values
        pblh = ds["PBLH"].values
        seconds, t0 = _time_axis(ds)

    # Destagger to mass points (WRF stores U/V/W on the faces).
    u_m = 0.5 * (u_stag[..., :-1] + u_stag[..., 1:])
    v_m = 0.5 * (v_stag[..., :-1, :] + v_stag[..., 1:, :])
    w_m = 0.5 * (w_stag[:, :-1] + w_stag[:, 1:])
    z_stag = geopotential / GRAVITY
    z_mass = 0.5 * (z_stag[:, :-1] + z_stag[:, 1:])
    if terrain.ndim == 3:  # WRF writes HGT with a Time axis
        terrain = terrain[:, None]
    z_agl = z_mass - terrain
    temperature = (theta_pert + THETA_BASE) * (pressure / P_REFERENCE) ** KAPPA

    n_y_w, n_x_w = u_m.shape[-2:]
    x_t = np.asarray(plume_grid.x, dtype=np.float64) + origin[0]
    y_t = np.asarray(plume_grid.y, dtype=np.float64) + origin[1]
    x_u = x_t + 0.5 * plume_grid.dx
    y_v = y_t + 0.5 * plume_grid.dy
    # WRF heights are above ground; the analysis grid's lower boundary is its
    # surface (``make_grid`` semantics), so interpolate at heights above
    # ``z_min = z[0] - dz/2`` rather than at the raw ``z`` coordinates.
    z_abs = np.asarray(plume_grid.z, dtype=np.float64)
    z_t = z_abs - (z_abs[0] - 0.5 * plume_grid.dz)
    _check_inside(x_u, y_v, dx_w=dx_w, dy_w=dy_w, n_x_w=n_x_w, n_y_w=n_y_w)
    _check_inside(x_t, y_t, dx_w=dx_w, dy_w=dy_w, n_x_w=n_x_w, n_y_w=n_y_w)

    def onto(field: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
        column = _interpolate_vertical(z_agl, field, z_t)
        return _bilinear(column, x / dx_w, y / dy_w)

    dtype = plume_grid.x.dtype
    as_jax = lambda a: jnp.asarray(a, dtype=dtype)
    return MetField(
        plume_grid=plume_grid,
        times=as_jax(seconds),
        u=as_jax(onto(u_m, x_u, y_t)),
        v=as_jax(onto(v_m, x_t, y_v)),
        w=as_jax(onto(w_m, x_t, y_t)),
        temperature=as_jax(onto(temperature, x_t, y_t)),
        pressure=as_jax(onto(pressure, x_t, y_t)),
        pbl_height=as_jax(_bilinear(pblh, x_t / dx_w, y_t / dy_w)),
        t0=t0,
    )