"""
HYSPLIT's config: every setting that shapes the particles HYSPLIT makes.
:class:`HysplitConfig` holds them. Most are HYSPLIT's own, under HYSPLIT's
names; a few are used by PYSTILT itself. Which input file each goes to is
the HYSPLIT driver's business (:mod:`stilt.transport.hysplit.driver`).
"""
from __future__ import annotations
from pathlib import Path
from typing import Any, ClassVar, Literal, Self
from pydantic import ConfigDict, Field, field_validator, model_validator
from stilt.particles import HNF_PLUME_COLUMNS
from stilt.transport import TransportConfig
#: Most hourly ZICONTROL factors HYSPLIT can hold (``ZIPRESC(150)`` in
#: hymodelc.F); it reads more without a bounds check.
MAX_ZISCALE_HOURS = 150
#: The wind-error settings, in the order of the lines of HYSPLIT's WINDERR file.
WIND_ERROR_SETTINGS = ("siguverr", "tluverr", "zcoruverr", "horcoruverr")
#: The mixed-layer error settings, in the order of the lines of HYSPLIT's ZIERR file.
ZI_ERROR_SETTINGS = ("sigzierr", "tlzierr", "horcorzierr")
[docs]
class HysplitConfig(TransportConfig):
"""
HYSPLIT's config: every setting that shapes the particles it makes.
In ``config.yaml`` they are flat, top-level keys when HYSPLIT is the
project's model, as it is by default, and a variant may override any of
them. Most are HYSPLIT ``SETUP.CFG`` entries with
HYSPLIT's own names; see the HYSPLIT user guide for the full meaning of
each. The wind-error group (``siguverr``, ``tluverr``, ``zcoruverr``,
``horcoruverr``) perturbs the particles' winds, and the mixed-layer group
(``sigzierr``, ``tlzierr``, ``horcorzierr``) perturbs each particle's
footprint by a random mixed-layer height error. Each group is set in full
or not at all.
"""
model_config = ConfigDict(arbitrary_types_allowed=True)
#: Fields left out of a run's recorded settings. What they point at is
#: recorded instead: the build's version and the checksums of data files
#: that differ from the bundled ones.
UNRECORDED: ClassVar[frozenset[str]] = frozenset({"exe_dir", "data_dir"})
numpar: int = Field(
200,
description=(
"Number of particles released per simulation. More particles give a "
"less noisy footprint and take longer to run."
),
)
exe_dir: Path | None = Field(
None,
description=(
"Directory holding a custom ``hycs_std`` build to run in place of "
"the one bundled with PYSTILT, with a ``version`` file beside it. "
"The version is recorded with each run, so runs of different "
"builds are kept apart. A build that writes release-time (t = 0) rows to "
"``PARTICLE_STILT.DAT`` gives exact release heights for multipoint "
"and slant receptors."
),
)
data_dir: Path | None = Field(
None,
description=(
"Directory of HYSPLIT data tables (``ASCDATA.CFG``, ``LANDUSE.ASC``, "
"``ROUGLEN.ASC``, ``TERRAIN.ASC``) to use in place of the bundled "
"ones. A table it does not hold comes from the bundled set. The "
"tables change the particles, so each table that differs from the "
"bundled one is recorded with the run by its checksum."
),
)
varsiwant: list[
Literal[
"time",
"indx",
"long",
"lati",
"zagl",
"sigw",
"tlgr",
"zsfc",
"icdx",
"temp",
"samt",
"foot",
"shtf",
"tcld",
"dmas",
"dens",
"rhfr",
"sphu",
"lcld",
"zloc",
"dswf",
"wout",
"mlht",
"rain",
"crai",
"pres",
"whtf",
"temz",
"zcfx",
]
] = Field(
default_factory=lambda: [
"time",
"indx",
"long",
"lati",
"zagl",
"foot",
"mlht",
"pres",
"dens",
"samt",
"sigw",
"tlgr",
],
description=(
"Particle variables ``hycs_std`` writes to the trajectory output, "
"by the lowercase four-letter names HYSPLIT reads. Each is a "
"column of the particles under the same name, except ``indx``, "
"``time``, ``long``, and ``lati``. ``zcfx`` is the vertical "
"displacement by convective flux along the trajectory, in m. "
"The default is the set footprints need, plus ``pres`` for "
"pressure weighting."
),
)
capemin: float = Field(
-1.0,
description=(
"Convection option. -1 turns convection off, -2 uses the Grell "
"scheme, and a positive value mixes vertically when CAPE exceeds "
"it, in J/kg."
),
)
cmass: int = Field(
0,
description="Compute grid concentrations (0) or grid mass (1).",
)
conage: int = Field(
48,
description="Particle age at which particles and puffs convert, in hours.",
)
cpack: int = Field(
1,
description="Packing of the binary concentration grid.",
)
delt: float = Field(
1.0,
description=(
"Integration time step, in minutes. 0 lets HYSPLIT choose; a "
"negative value sets the minimum step."
),
)
dxf: float = Field(
1.0,
description="Horizontal x-grid offset factor for ensemble runs.",
)
dyf: float = Field(
1.0,
description="Horizontal y-grid offset factor for ensemble runs.",
)
dzf: float = Field(
0.01,
description="Vertical offset factor for ensemble runs (0.01 is about 250 m).",
)
efile: str = Field(
"",
description="Name of a time-varying emissions file. Blank uses none.",
)
emisshrs: float = Field(
0.01,
description="Duration of the particle release, in hours.",
)
frhmax: float = Field(
3.0,
description="Maximum horizontal puff-rounding parameter.",
)
frhs: float = Field(
1.0,
description="Horizontal puff-rounding fraction for merging.",
)
frme: float = Field(
0.1,
description="Mass-rounding fraction for enhanced merging.",
)
frmr: float = Field(
0.0,
description="Mass-removal fraction for enhanced merging.",
)
frts: float = Field(0.1, description="Temporal puff-rounding fraction.")
frvs: float = Field(0.01, description="Vertical puff-rounding fraction.")
hscale: float = Field(
10800.0,
description="Horizontal Lagrangian timescale, in seconds.",
)
ichem: int = Field(
8,
description="HYSPLIT chemistry and output mode. 8 is the STILT emulation mode.",
)
idsp: int = Field(
2,
description="Particle dispersion scheme: 1 for HYSPLIT, 2 for STILT.",
)
initd: int = Field(
0,
description="Initial distribution as particles, puffs, or a mix. 0 is 3D particles.",
)
k10m: int = Field(
1,
description=(
"Use the 10 m winds and 2 m temperature as the lowest meteorology "
"level (1) or skip them (0)."
),
)
kagl: int = Field(
1,
description="Write trajectory heights above ground (1) or above sea level (0).",
)
kbls: int = Field(
1,
description=(
"Derive boundary-layer stability from surface fluxes (1) or from "
"wind and temperature profiles (2)."
),
)
kblt: int = Field(
5,
description=(
"Boundary-layer turbulence scheme: 1 Beljaars, 2 Kantha-Clayson, "
"3 TKE, 4 measured variances, 5 Hanna."
),
)
kdef: int = Field(
0,
description="Horizontal turbulence from vertical mixing (0) or wind deformation (1).",
)
khinp: int = Field(
0,
description="Age, in hours, given to particles read from ``pinpf``. 0 keeps their own age.",
)
khmax: int = Field(
9999,
description="Maximum particle or trajectory age, in hours.",
)
kmix0: int = Field(
150,
description="Minimum mixed-layer depth, in meters.",
)
kmixd: int = Field(
3,
description=(
"Mixed-layer depth source: 0 from the meteorology, 1 from the "
"temperature profile, 2 from the TKE profile, 3 from a modified "
"Richardson number."
),
)
kpuff: int = Field(
0,
description="Horizontal puff growth: linear (0) or empirical (1).",
)
krand: Literal[0, 1, 2, 3, 4, 10, 11, 12, 13] = Field(
4,
description=(
"How HYSPLIT draws the random numbers for turbulence. 0 picks 2 "
"when ``numpar`` is 5000 or less and 1 otherwise. 1 uses a "
"precomputed table. 2 draws them during the run and is the only "
"mode that uses ``seed``. 3 uses no random numbers (a diagnostic "
"mode). 4 draws them during the run from a clock-based seed, so "
"every run differs; there are about 5000 possible seeds. 10 to 13 "
"are modes 0 to 3 with a random initial seed. HYSPLIT does not "
"check this value and other values silently break the turbulence, "
"so PYSTILT rejects them."
),
)
seed: int | None = Field(
None,
description=(
"Random seed for a reproducible run. Different seeds give "
"different runs. Requires ``krand: 2``. The bundled ``hycs_std`` "
"ignores the seed under ``krand`` 4 and 10 to 13, and under 1 uses "
"it only for the initial turbulent velocity. PYSTILT writes it to "
"``SETUP.CFG`` as ``-(abs(seed) + 1)``, because HYSPLIT reseeds its "
"generator only from a negative value and gives every ``SEED`` of "
"0 or more the same stream. Realization ``k`` of a variant runs "
"with ``seed + k``, so realization 0 shares the unperturbed run's "
"seed, as STILT-R's error run does."
),
)
krnd: int = Field(6, description="Enhanced-merging interval, in hours.")
kspl: int = Field(
1,
description="Standard puff-splitting interval, in hours.",
)
kwet: int = Field(
1,
description="Precipitation from the meteorology (1) or from an external ARL file (2).",
)
kzmix: int = Field(
0,
description=(
"Vertical mixing adjustment: 0 none, 1 a single PBL-average value, "
"2 scale by ``tvmix``."
),
)
maxdim: int = Field(
1,
description="Maximum number of pollutant species carried on one particle.",
)
maxpar: int | None = Field(
None,
description=(
"Maximum number of particles in a simulation. Unset gives every "
"release location its full share: ``numpar`` for one location, "
"``n * (numpar // n + 1)`` for n. A set value that stops a "
"multipoint release before its last location fails the run."
),
)
mgmin: int = Field(
10,
description="Minimum meteorological subgrid size, in grid points.",
)
mhrs: int = Field(
9999,
description="Trajectory restart duration limit, in hours.",
)
nbptyp: int = Field(
1,
description="Number of particle-size bins per pollutant type.",
)
ncycl: int = Field(
0,
description="Cycle time of the particle dump file, in hours.",
)
ndump: int = Field(
0,
description="Interval between particle dumps, in hours. 0 writes none.",
)
ninit: int = Field(
1,
description=(
"Particle initialization from ``pinpf``: 0 none, 1 once at the "
"start, 2 add every hour, 3 replace every hour."
),
)
nstr: int = Field(0, description="Trajectory restart interval, in hours.")
nturb: int = Field(
0,
description="Turbulence on (0) or off (1).",
)
nver: int = Field(0, description="Trajectory vertical split number.")
outdt: int = Field(
0,
description=(
"Interval between particle outputs in ``PARTICLE_STILT.DAT``, in "
"minutes. 0 writes every time step and a negative value writes none. "
"Only at 0 does PYSTILT check where the particles of a run that "
"stopped early ended, to tell met that ran out from particles that "
"left the met's domain."
),
)
p10f: float = Field(
1.0,
description="Dust threshold velocity sensitivity factor.",
)
pinbc: str = Field(
"",
description="Particle input file for time-varying boundary conditions.",
)
pinpf: str = Field(
"",
description="Particle input file for initialization or boundary-condition runs.",
)
poutf: str = Field(
"",
description="Particle output file name.",
)
qcycle: float = Field(
0.0,
description="Emission cycling period, in hours. 0 turns cycling off.",
)
rhb: int = Field(
80,
description="Relative humidity that defines a cloud base, in percent.",
)
rht: int = Field(
60,
description="Relative humidity that defines a cloud top, in percent.",
)
splitf: float = Field(
1.0,
description=(
"Factor for the automatic horizontal splitting size. A negative "
"value turns the automatic sizing off."
),
)
tkerd: float = Field(
0.18,
description="Ratio w'²/(u'²+v'²) of TKE components when unstable.",
)
tkern: float = Field(
0.18,
description="Ratio w'²/(u'²+v'²) of TKE components when stable.",
)
tlfrac: float = Field(
0.1,
description=(
"Fraction of the vertical Lagrangian timescale used as the time "
"step of the STILT dispersion scheme."
),
)
tout: int = Field(
0,
description="Trajectory output interval, in minutes.",
)
tratio: float = Field(
0.75,
description="Advection stability ratio (fraction of a grid cell per time step).",
)
tvmix: float = Field(
1.0,
description="Vertical mixing scale factor, used by the ``kzmix`` scaling modes.",
)
vscale: float = Field(
200.0,
description="Vertical Lagrangian timescale, in seconds.",
)
vscaleu: float = Field(
200.0,
description="Vertical Lagrangian timescale in an unstable boundary layer, in seconds.",
)
vscales: float = Field(
-1.0,
description=(
"Vertical Lagrangian timescale in a stable boundary layer, in "
"seconds. -1 uses the Hanna timescale, which varies with the "
"turbulence, and then ``vscaleu`` is not used."
),
)
w_option: int = Field(
0,
description=(
"Vertical motion method: 0 the meteorology's vertical velocity, "
"1 isobaric, 2 isentropic, 3 constant density, 4 constant sigma."
),
)
wbbh: float = Field(
0.0,
description=(
"Height at which the fixed vertical velocity switches from rise to "
"fall, in meters. Used by vertical motion option 9."
),
)
wbwf: float = Field(
0.0,
description="Fixed fall velocity, in m/s. Used by vertical motion options 9 and 10.",
)
wbwr: float = Field(
0.0,
description="Fixed rise velocity, in m/s. Used by vertical motion option 9.",
)
wvert: bool = Field(
False,
description="Interpolate WRF fields vertically with the WRF scheme instead of HYSPLIT's.",
)
z_top: float = Field(
25000.0,
description="Top of HYSPLIT's computational domain, in meters above ground.",
)
ziscale: float | list[float] = Field(
1.0,
description=(
"Factor applied to the mixed-layer height, written to HYSPLIT's "
"``ZICONTROL`` file. 1.0 leaves it unscaled. A single value applies "
"to every hour of the run. A list gives one factor per hour from "
"the release (at most 150), and later hours are unscaled. HYSPLIT "
"applies ``kmix0`` after the factor, so the mixed layer never drops "
"below ``kmix0``. A negative value uses the meteorology's own PBL "
"height where the met files carry one."
),
)
siguverr: float | None = Field(
None,
description="Standard deviation of the horizontal wind error, in m/s.",
)
tluverr: float | None = Field(
None,
description="Correlation timescale of the horizontal wind error, in minutes.",
)
zcoruverr: float | None = Field(
None,
description="Vertical correlation length of the horizontal wind error, in meters.",
)
horcoruverr: float | None = Field(
None,
description="Horizontal correlation length of the horizontal wind error, in km.",
)
sigzierr: float | None = Field(
None,
description="Standard deviation of the mixed-layer height error, in percent.",
)
tlzierr: float | None = Field(
None,
description="Correlation timescale of the mixed-layer height error, in minutes.",
)
horcorzierr: float | None = Field(
None,
description="Horizontal correlation length of the mixed-layer height error, in km.",
)
@field_validator("varsiwant", mode="before")
@classmethod
def _refuse_zfx1(cls, value: Any) -> Any:
"""Refuse ``zfx1``, which ``hycs_std`` stops on, and name ``zcfx``, which it reads."""
if isinstance(value, list) and any(
isinstance(name, str) and name.lower() == "zfx1" for name in value
):
raise ValueError(
"'zfx1' is not a varsiwant name: hycs_std stops with 'Weird "
"variable'. HYSPLIT reads 'zcfx' for the vertical displacement "
"by convective flux."
)
return value
@field_validator("ziscale", mode="before")
@classmethod
def _flatten_ziscale(cls, value: Any) -> Any:
"""
Read STILT-R's one-element nested list, ``[[0.8, 0.9]]``, as ``[0.8, 0.9]``.
Both mean the same factors, so they must hash the same.
"""
if isinstance(value, list) and value and isinstance(value[0], list):
if len(value) != 1:
raise ValueError(
"Per-simulation ziscale lists are not supported. Pass one shared "
"list of hourly factors for all simulations."
)
return value[0]
return value
@model_validator(mode="after")
def _validate_error_params(self) -> Self:
"""Require each error group to be set in full or not at all."""
for name, fields in (("XY", WIND_ERROR_SETTINGS), ("ZI", ZI_ERROR_SETTINGS)):
unset = [getattr(self, f) is None for f in fields]
if any(unset) and not all(unset):
raise ValueError(
f"Inconsistent {name} error parameters: all must be set or all None"
)
return self
@model_validator(mode="after")
def _validate_ziscale(self) -> Self:
"""Reject ziscale values that are empty, zero, or longer than 150 hours."""
hourly = self.hourly_ziscale()
if not hourly:
raise ValueError("ziscale cannot be empty; use 1.0 for no scaling.")
if any(v == 0.0 for v in hourly):
raise ValueError(
"ziscale of 0 would collapse the mixed layer to kmix0. STILT-R "
"uses 0 to mean unset; use 1.0 for no scaling."
)
if any(v != 1.0 for v in hourly) and len(hourly) > MAX_ZISCALE_HOURS:
raise ValueError(
f"ziscale gives {len(hourly)} hourly factors, but HYSPLIT holds at "
f"most {MAX_ZISCALE_HOURS}. A scalar ziscale is repeated for "
"every hour, so it needs abs(n_hours) <= "
f"{MAX_ZISCALE_HOURS}; for longer runs give a list of up to "
f"{MAX_ZISCALE_HOURS} factors, after which the mixed layer "
"is unscaled."
)
return self
@model_validator(mode="after")
def _validate_seed(self) -> Self:
"""Require ``krand=2`` when a seed is set."""
if self.seed is not None and self.krand != 2:
raise ValueError(
f"seed={self.seed} requires krand=2 (got krand={self.krand}): the "
"bundled hycs_std discards the namelist seed under krand=4 and "
"10-13 and uses it only for the initial turbulent velocity under "
"krand=1. Set krand=2, or drop the seed."
)
return self
@model_validator(mode="after")
def _validate_hnf_plume(self) -> Self:
"""Require the variables the near-field plume model reads when ``hnf_plume`` is on."""
if self.hnf_plume:
missing = set(HNF_PLUME_COLUMNS) - set(self.varsiwant)
if missing:
raise ValueError(
f"hnf_plume=True requires varsiwant to include: {sorted(missing)}"
)
return self
[docs]
def maxpar_for(self, locations: int) -> int:
"""
Return the ``MAXPAR`` HYSPLIT gets for a release from *locations* horizontal locations.
HYSPLIT gives each location ``numpar // locations`` particles, which
its rounding often raises by one, and stops releasing when it reaches
``MAXPAR``. A set ``maxpar`` is used as given. Unset, one location
gets ``numpar`` (HYSPLIT may ask for one more and keeps ``numpar``),
and several get ``locations * (numpar // locations + 1)``, so the
last locations are not left short.
Examples
--------
>>> HysplitConfig(numpar=1000).maxpar_for(40)
1040
>>> HysplitConfig(numpar=1000).maxpar_for(1)
1000
"""
if self.maxpar is not None:
return self.maxpar
if locations <= 1:
return self.numpar
return locations * (self.numpar // locations + 1)
[docs]
def hourly_ziscale(self) -> list[float]:
"""Return ``ziscale`` as one factor per hour: a single value repeated over the run, or the list as given."""
if isinstance(self.ziscale, int | float):
return [float(self.ziscale)] * max(abs(self.n_hours), 1)
return [float(v) for v in self.ziscale]
[docs]
def realizations(self, n: int) -> list[Self]:
"""
Return *n* realizations of this config, realization ``k`` with ``seed + k``.
Realization 0 uses the configured seed, as STILT-R's single error run
does.
Raises
------
ValueError
If *n* is more than one and the runs would repeat one another:
only ``krand=4`` (HYSPLIT seeds each run from the clock) and
``krand=2`` with a seed give each realization its own draws.
"""
if n > 1 and not (
self.krand == 4 or (self.krand == 2 and self.seed is not None)
):
raise ValueError(
f"realizations={n} requires krand=4 or krand=2 with a seed (got "
f"krand={self.krand}, seed={self.seed}): under krand=4 HYSPLIT "
"seeds each run from the clock; under krand=2 PYSTILT gives each "
"realization its own seed. Any other mode would repeat the same "
"perturbation."
)
return super().realizations(n)
__all__ = [
"MAX_ZISCALE_HOURS",
"WIND_ERROR_SETTINGS",
"ZI_ERROR_SETTINGS",
"HysplitConfig",
]