The bldfm package

This package contains the BLDFM model. It includes core modules for the PBL model, spectral solver, and utility functions, as well as higher-level modules for configuration, interface, caching, and synthetic data generation. Plotting and I/O are provided by abltk.

Core modules

bldfm.pbl_model module

bldfm.pbl_model.phi(x)

Stability correction function for eddy diffusivity in Monin-Obukhov Similarity Theory (MOST).

Parameters:

x (float or numpy.ndarray) – Dimensionless stability parameter (z / mol).

Returns:

Stability correction factor for eddy diffusivity.

Return type:

numpy.ndarray

Notes

• For stable conditions (x > 0), a linear relationship is used.

• For unstable conditions (x < 0), the Businger-Dyer formulation is applied.

bldfm.pbl_model.psi(x)

Stability correction function for Monin-Obukhov Similarity Theory (MOST).

Parameters:

x (float or numpy.ndarray) – Dimensionless stability parameter (z / mol).

Returns:

Stability correction factor for momentum.

Return type:

numpy.ndarray

Notes

• For stable conditions (x > 0), a linear relationship is used.

• For unstable conditions (x < 0), the Businger-Dyer formulation is applied.

bldfm.pbl_model.vertical_profiles(n, meas_height, wind, ustar=None, z0=None, mol=1000000000.0, prsc=1.0, closure='MOST', domain_height=None, stretch=None, z0_min=0.001, z0_max=2.0, tke=None)

Computes vertical profiles of horizontal wind components and eddy diffusivity in the planetary boundary layer (PBL) based on Monin-Obukhov Similarity Theory (MOST) or other closure models.

Parameters:

• n (int) – Number of vertical grid points between z0 and meas_height.

• meas_height (float) – Measurement height above the ground.

• wind (tuple of floats) – Zonal (u) and meridional (v) wind components at the measurement height.

• ustar (float or numpy.ndarray) – Friction velocity [m/s].

• mol (float or numpy.ndarray, optional) – Monin-Obukhov length. Default is 1e9 (neutral conditions).

• prsc (float, optional) – Prandtl or Schmidt number. Default is 1.0.

• closure (str, optional) – Closure model to use. Options are “MOST”, “CONSTANT”, or “OAAHOC”. Default is “MOST”.

• z0 (float or numpy.ndarray, optional) – Roughness length. Default is -1e9 (auto-calculated).

• z0_min (float, optional) – Minimum allowable roughness length. Default is 0.001.

• z0_max (float, optional) – Maximum allowable roughness length. Default is 2.0.

• tke (float or numpy.ndarray, optional) – Turbulent kinetic energy [m²/s²] for the “OAAHOC” closure. If not provided, it will default to 1.0.

Returns:

• z (numpy.ndarray): 1D array of vertical grid points.

• profiles (tuple of numpy.ndarray): 1D arrays of horizontal wind components (u, v) and eddy diffusivities (Kx, Ky, Kz) at each vertical grid point.

Return type:

tuple

Raises:

ValueError – If invalid closure type is provided.

Notes

• The “OAAHOC” closure uses a one-and-a-half order closure model based on Schumann-Lilly.

• The “MOST” closure uses Monin-Obukhov Similarity Theory.

References

• Kormann, R., & Meixner, F. X. (2001). An analytical footprint model for non-neutral stratification. Boundary-Layer Meteorology, 99(2), 207–224. https://doi.org/10.1023/A:1018991015119

• Schumann, U. (1991). Subgrid length-scales for large-eddy simulation of stratified turbulence. Theoretical and Computational Fluid Dynamics, 2(5), 279–290. https://doi.org/10.1007/BF00271468

bldfm.solver module

bldfm.solver.steady_state_transport_solver(srf_flx, z, profiles, domain, levels, modes=(512, 512), meas_pt=(0.0, 0.0), srf_bg_conc=0.0, footprint=False, analytic=False, halo=None, precision='single', cache=None)

Solves the steady-state advection-diffusion equation for a concentration with flux boundary condition given vertical profiles of wind and eddy diffusivity using the Fourier, linear shooting, and semi-implicit Euler methods.

Parameters:

• srf_flx (ndarray of float) – 2D field of surface kinematic flux at z=z0 [m/s].

• z (ndarray of float) – 1D array of vertical grid points from z0 to zm [m].

• profiles (tuple of ndarray) – Tuple containing 1D arrays of vertical profiles of zonal wind, meridional wind [m/s], and eddy diffusivities [m²/s].

• domain (tuple of float) – Tuple containing domain sizes (xmax, ymax) [m].

• levels (float or ndarray of float) – Vertical level for output or optionally 1D array of output levels.

• modes (tuple of int, optional) – Tuple containing the number of zonal and meridional Fourier modes (nlx, nly). Default is (512, 512).

• meas_pt (tuple of float, optional) – Coordinates of the measurement point (xm, ym) [m] relative to srf_flx, where the origin is at srf_flx[0, 0]. Default is (0.0, 0.0).

• srf_bg_conc (float, optional) – Surface background concentration at z=z0 [scalar_unit]. Default is 0.0.

• footprint (bool, optional) – If True, activates footprints (Green’s function) for output. Default is False.

• analytic (bool, optional) – If True, uses the analytic solution for constant wind and eddy diffusivity. Default is False.

• halo (float, optional) – Width of the zero-flux halo around the domain [m]. Default is -1e9, which sets the halo to max(xmax, ymax).

Returns:

• z (ndarray of float) – Heights [m] at levels.

• conc (ndarray of float) – 2D or 3D field of concentration at levels or Green’s function.

• flx (ndarray of float) – 2D or 3D field of kinematic flux at levels or footprint.

bldfm.utils module

bldfm.utils.compute_wind_fields(u_rot, wind_dir)

Computes the zonal (u) and meridional (v) wind components from a rotated wind speed and direction using the meteorological convention.

Parameters:

• u_rot (float) – Rotated wind speed.

• wind_dir (float) – Wind direction in degrees (meteorological convention: direction the wind is coming FROM, clockwise from north). 0 = from north, 90 = from east, 180 = from south, 270 = from west.

Returns:

A tuple (u, v) where:

• u (float): Zonal wind component (east-west, positive = eastward).

• v (float): Meridional wind component (north-south, positive = northward).

Return type:

tuple

bldfm.utils.ideal_source(nxy, domain, src_loc=None, shape='diamond')

Creates a synthetic source field in the shape of a circle or diamond. Useful for testing purposes.

Parameters:

• nxy (tuple) – Number of grid points in the x and y directions (nx, ny).

• domain (tuple) – Physical dimensions of the domain (xmax, ymax).

• shape (str) – Shape of the source field. Options are “circle” or “diamond”. Default is “diamond”.

Returns:

A 2D array representing the source field.

Return type:

numpy.ndarray

bldfm.utils.parallelize(func)

bldfm.utils.point_measurement(f, g)

Computes the convolution of two 2D arrays evaluated at a specific point.

Parameters:

• f (numpy.ndarray) – First 2D array.

• g (numpy.ndarray) – Second 2D array.

Returns:

The result of the convolution at the specified point.

Return type:

float

bldfm.fft_manager module

class bldfm.fft_manager.FFTManager(wisdom_file='fftw_wisdom.pkl', num_threads=1, cache_keepalive=30)

Bases: object

Manages pyFFTW with wisdom and caching for optimal performance in Dask environments.

This class addresses memory issues in Dask parallelized environments by: - Using pyfftw numpy interface for compatibility - Loading/saving FFTW wisdom for optimal planning - Enabling pyfftw caching to prevent repeated object allocation - Providing proper memory management and cleanup

clear_cache()

Clear pyfftw cache.

fft2(input_data, norm='backward')

Perform 2D forward FFT using pyfftw numpy interface.

Parameters:

• input_data – Input array for FFT

• norm – Normalization mode (“forward”, “backward”, “ortho”)

Returns:

Complex array with FFT result

ifft2(input_data, norm='backward')

Perform 2D inverse FFT using pyfftw numpy interface.

Parameters:

• input_data – Input array for inverse FFT

• norm – Normalization mode (“forward”, “backward”, “ortho”)

Returns:

Complex array with inverse FFT result

bldfm.fft_manager.fft2(input_data, norm='backward')

Global function for 2D FFT using the FFT manager.

bldfm.fft_manager.get_fft_manager(num_threads=1, cache_keepalive=30)

Get or create the global FFT manager instance.

If the manager already exists with different num_threads, it is re-initialised with the new thread count.

bldfm.fft_manager.ifft2(input_data, norm='backward')

Global function for 2D inverse FFT using the FFT manager.

bldfm.fft_manager.reset_fft_manager()

Reset FFTManager singleton (call in forked worker processes).

Configuration and interface

bldfm.config_parser module

YAML configuration parser and dataclass schema for BLDFM.

Defines the configuration hierarchy:

BLDFMConfig ├── DomainConfig ├── List[TowerConfig] ├── MetConfig ├── SolverConfig ├── OutputConfig └── ParallelConfig

class bldfm.config_parser.BLDFMConfig(domain: ~bldfm.config_parser.DomainConfig, towers: ~typing.List[~bldfm.config_parser.TowerConfig], met: ~bldfm.config_parser.MetConfig, solver: ~bldfm.config_parser.SolverConfig = <factory>, output: ~bldfm.config_parser.OutputConfig = <factory>, parallel: ~bldfm.config_parser.ParallelConfig = <factory>)

Bases: object

Top-level BLDFM configuration.

domain: DomainConfig

met: MetConfig

output: OutputConfig

parallel: ParallelConfig

solver: SolverConfig

towers: List[TowerConfig]

class bldfm.config_parser.DomainConfig(nx: int, ny: int, xmax: float, ymax: float, nz: int, modes: Tuple[int, int] = (512, 512), halo: float | None = None, ref_lat: float | None = None, ref_lon: float | None = None, output_levels: List[int] | None = None, full_output: bool = False)

Bases: object

Configuration for the computational domain.

full_output: bool = False

halo: float | None = None

modes: Tuple[int, int] = (512, 512)

nx: int

ny: int

nz: int

output_levels: List[int] | None = None

ref_lat: float | None = None

ref_lon: float | None = None

xmax: float

ymax: float

class bldfm.config_parser.MetConfig(ustar: float | List[float] | None = None, mol: float | List[float] = 1000000000.0, wind_speed: float | List[float] = 5.0, wind_dir: float | List[float] = 270.0, z0: float | None = None, timestamps: List[str] | None = None)

Bases: object

Meteorological forcing data (scalar for single-time, list for timeseries).

get_step(i: int) → dict

Get met parameters for a single timestep.

Return type:

dict with scalar values for ustar, mol, wind_speed, wind_dir, timestamp.

mol: float | List[float] = 1000000000.0

property n_timesteps: int

Number of timesteps in the timeseries.

timestamps: List[str] | None = None

ustar: float | List[float] | None = None

validate()

Validate that at least one of ustar/z0 is provided, and list lengths match.

wind_dir: float | List[float] = 270.0

wind_speed: float | List[float] = 5.0

z0: float | None = None

class bldfm.config_parser.OutputConfig(format: str = 'netcdf', directory: str = './output')

Bases: object

Output configuration.

directory: str = './output'

format: str = 'netcdf'

class bldfm.config_parser.ParallelConfig(num_threads: int = 1, max_workers: int = 1, use_cache: bool = False)

Bases: object

Parallelism configuration.

max_workers: int = 1

num_threads: int = 1

use_cache: bool = False

class bldfm.config_parser.SolverConfig(closure: str = 'MOST', precision: str = 'single', footprint: bool = False, surface_flux_shape: str = 'diamond', analytic: bool = False, src_loc: Tuple[float, float] | None = None)

Bases: object

Solver configuration.

analytic: bool = False

closure: str = 'MOST'

footprint: bool = False

precision: str = 'single'

src_loc: Tuple[float, float] | None = None

surface_flux_shape: str = 'diamond'

class bldfm.config_parser.TowerConfig(name: str, lat: float, lon: float, z_m: float, x: float = 0.0, y: float = 0.0)

Bases: object

Configuration for a single measurement tower.

compute_local_xy(ref_lat: float, ref_lon: float)

Compute local x/y from lat/lon and reference origin.

lat: float

lon: float

name: str

x: float = 0.0

y: float = 0.0

z_m: float

bldfm.config_parser.load_config(path: str | Path) → BLDFMConfig

Load a BLDFM configuration from a YAML file.

Parameters:

path (str or Path) – Path to the YAML configuration file.

Returns:

Parsed and validated configuration.

Return type:

BLDFMConfig

bldfm.config_parser.parse_config_dict(raw: dict) → BLDFMConfig

Parse a BLDFM configuration from a dictionary.

Parameters:

raw (dict) – Raw configuration dictionary (e.g. from YAML).

Returns:

Parsed and validated configuration.

Return type:

BLDFMConfig

bldfm.interface module

High-level interface for running BLDFM simulations.

Provides convenience functions that wrap the manual workflow of compute_wind_fields -> vertical_profiles -> steady_state_transport_solver into single function calls driven by configuration objects.

bldfm.interface.run_bldfm_multitower(config: BLDFMConfig, surface_flux: ndarray | None = None) → dict

Run BLDFM for all towers and all timesteps.

Parameters:

• config (BLDFMConfig) – Full BLDFM configuration with towers and met data.

• surface_flux (ndarray, optional) – 2D surface flux field (reused for all towers/timesteps).

Returns:

Mapping of tower_name -> list of result dicts (one per timestep).

Return type:

dict

bldfm.interface.run_bldfm_parallel(config: BLDFMConfig, max_workers: int | None = None, parallel_over: str = 'towers', surface_flux: ndarray | None = None) → dict

Run BLDFM in parallel using ProcessPoolExecutor.

Parameters:

• config (BLDFMConfig) – Full BLDFM configuration with towers and met data.

• max_workers (int, optional) – Number of worker processes. Defaults to config.parallel.max_workers.

• parallel_over (str) – Parallelization strategy: - “towers”: distribute towers across workers, each runs full timeseries - “time”: for each tower, distribute timesteps across workers - “both”: flatten all (tower, timestep) pairs across workers

• surface_flux (ndarray, optional) – 2D surface flux field (not passed to subprocesses to avoid large serialization; each worker generates its own ideal source).

Returns:

Mapping of tower_name -> list of result dicts (one per timestep). Same format as run_bldfm_multitower.

Return type:

dict

bldfm.interface.run_bldfm_single(config: BLDFMConfig, tower: TowerConfig, met_index: int = 0, surface_flux: ndarray | None = None, cache=None) → dict

Run a single BLDFM solve for one tower at one timestep.

Encapsulates the 3-step workflow:

1. compute_wind_fields(wind_speed, wind_dir) -> (u, v)

2. vertical_profiles(nz, z_m, (u, v), …) -> (z, profiles)

3. steady_state_transport_solver(srf_flx, z, profiles, …) -> (grid, conc, flx)

Parameters:

• config (BLDFMConfig) – Full BLDFM configuration.

• tower (TowerConfig) – Tower to compute footprint/concentration for.

• met_index (int) – Index into the met timeseries (0 for single-timestep configs).

• surface_flux (ndarray, optional) – 2D surface flux field. If None, generates an ideal source.

• cache (GreensFunctionCache, optional) – Cache instance for footprint reuse.

Returns:

Result dictionary with keys: - grid: (X, Y, Z) coordinate arrays - conc: concentration field - flx: flux field - tower_name: name of the tower - timestamp: timestamp or index - params: dict of met parameters used

Return type:

dict

bldfm.interface.run_bldfm_timeseries(config: BLDFMConfig, tower: TowerConfig, surface_flux: ndarray | None = None) → list

Run BLDFM for all timesteps in the met config for a single tower.

Parameters:

• config (BLDFMConfig) – Full BLDFM configuration with timeseries met data.

• tower (TowerConfig) – Tower to compute footprint/concentration for.

• surface_flux (ndarray, optional) – 2D surface flux field (reused for all timesteps).

Returns:

One result dict per timestep (same format as run_bldfm_single).

Return type:

list of dict

bldfm.cli module

Command-line interface for BLDFM.

Usage:

bldfm run config.yaml bldfm run config.yaml –dry-run bldfm run config.yaml –plot

bldfm.cli.cmd_run(args)

Execute a BLDFM run from a YAML config file.

bldfm.cli.main()

Main CLI entry point.

Data and I/O

bldfm.synthetic module

Synthetic data generators for testing and demonstration.

Provides functions to generate realistic-looking meteorological timeseries and tower grid configurations without requiring real observational data.

bldfm.synthetic.generate_synthetic_timeseries(n_timesteps: int = 48, dt_minutes: int = 30, start_time: str = '2024-01-01T00:00', ustar_range: Tuple[float, float] = (0.1, 0.8), mol_range: Tuple[float, float] = (-500.0, 500.0), wind_speed_range: Tuple[float, float] = (1.0, 8.0), wind_dir_mean: float = 270.0, wind_dir_std: float = 30.0, seed: int | None = None) → dict

Generate a synthetic meteorological timeseries with diurnal cycle.

Produces ustar, Monin-Obukhov length, wind speed, and wind direction arrays that follow a realistic diurnal pattern (unstable daytime, stable nighttime) with added noise.

Parameters:

• n_timesteps (int) – Number of time steps to generate.

• dt_minutes (int) – Time step interval in minutes.

• start_time (str) – ISO-format start time string.

• ustar_range (tuple of float) – (min, max) friction velocity [m/s].

• mol_range (tuple of float) – (min_negative, max_positive) Monin-Obukhov length [m]. Negative = unstable, positive = stable.

• wind_speed_range (tuple of float) – (min, max) wind speed [m/s].

• wind_dir_mean (float) – Mean wind direction [degrees].

• wind_dir_std (float) – Standard deviation of wind direction [degrees].

• seed (int, optional) – Random seed for reproducibility.

Returns:

Dictionary with keys matching MetConfig schema: ustar, mol, wind_speed, wind_dir, timestamps (all lists).

Return type:

dict

bldfm.synthetic.generate_towers_grid(n_towers: int = 4, center_lat: float = 50.95, center_lon: float = 11.586, spacing_m: float = 500.0, z_m: float = 10.0, layout: str = 'grid', seed: int | None = None) → List[dict]

Generate tower configurations in various spatial layouts.

Parameters:

• n_towers (int) – Number of towers to generate.

• center_lat (float) – Center point of the tower array [decimal degrees].

• center_lon (float) – Center point of the tower array [decimal degrees].

• spacing_m (float) – Approximate spacing between towers [meters].

• z_m (float) – Measurement height [meters] (same for all towers).

• layout (str) – Spatial layout: “grid”, “random”, or “transect”.

• seed (int, optional) – Random seed for reproducibility (used by “random” layout).

Returns:

Tower configurations matching TowerConfig schema.

Return type:

list of dict

bldfm.cache module

Disk cache for Green’s function results.

When footprint=True, the solver output depends only on the vertical profiles, domain geometry, and measurement point — not on the surface flux field. Caching these results avoids redundant solves when re-running with the same PBL/domain configuration.

class bldfm.cache.GreensFunctionCache(cache_dir='.bldfm_cache')

Bases: object

Disk-based cache for Green’s function solver outputs.

Cache key is a SHA-256 hash of the solver inputs that determine the Green’s function: vertical grid, profiles, domain, modes, measurement point, halo, and precision.

Parameters:

cache_dir (str or Path) – Directory for cache files. Created if it does not exist.

clear()

Remove all cached files.

get(z, profiles, domain, modes, meas_pt, halo, precision)

Look up cached result.

Returns:

(grid, conc, flx) if cached, None on miss.

Return type:

tuple or None

put(z, profiles, domain, modes, meas_pt, halo, precision, grid, conc, flx)

Store a result in the cache.

Plotting

Plotting functions have moved to abltk.plotting. See the abl-tk documentation for the full API reference. Common imports:

from abltk.plotting import plot_footprint_field, plot_footprint_on_map