Dynamic tropopause map¶
Ertel potential vorticity is computed on 23 pressure levels from the GFS temperature and wind, using the isobaric form with spherical derivatives and the curvature term in the relative vorticity. Each column is searched from the top down for the first level below 2 PVU, and potential temperature and wind are interpolated to that crossing. Low values mark troughs, where stratospheric air extends down. Black contours show cyclonic vorticity in the 925-850 hPa layer.
Data: GFS 0.5°, 00 UTC 27 October 2026, six-hour forecast, read with fastgrib from NOAA's archive on AWS by byte range. The loading and finite-difference code is shared with the other diagnostics in docs/examples/gfs3d.py.
docs/examples/dynamic_tropopause.py
"""The dynamic tropopause: potential temperature and wind on the 2 PVU surface,
with cyclonic vorticity in the 925-850 hPa layer.
Ertel potential vorticity is computed on each of 23 pressure levels from the
GFS temperature and wind (isobaric form, hydrostatic):
PV = -g [ (f + zeta) dtheta/dp - (dv/dp dtheta/dx - du/dp dtheta/dy) ]
Each column is then searched from the top down for the first level where PV
falls below 2 PVU, and potential temperature and wind are interpolated
linearly in PV to that crossing. Searching from the top skips the low-level
PV made by latent heating. Low potential temperature on the tropopause marks
troughs, where stratospheric air reaches down; high values mark ridges.
"""
import os
import numpy as np
import fastplot.crs as ccrs
import fastplot.pyplot as plt
import gfs3d
grib = gfs3d.open_gfs()
lat, lon = gfs3d.axes(grib)
p = np.array(gfs3d.LEVELS, np.float64)
t, u, v = (gfs3d.levels(grib, n) for n in ('TMP', 'UGRD', 'VGRD'))
th = gfs3d.theta(t, p)
zeta = gfs3d.vorticity(u, v, lat, lon)
pv = -gfs3d.G * ((gfs3d.coriolis(lat) + zeta) * gfs3d.ddp(th, p)
- (gfs3d.ddp(v, p) * gfs3d.ddx(th, lat, lon)
- gfs3d.ddp(u, p) * gfs3d.ddy(th, lat))) * 1e6 # PVU
# Top down: the first level below 2 PVU, and the one above it.
pv_top, th_top, u_top, v_top = (a[::-1] for a in (gfs3d.smooth(pv, 2), th, u, v))
below = pv_top < 2.0
k = np.argmax(below, axis=0)
found = below.any(axis=0) & (k > 0)
k = np.clip(k, 1, None)
take = lambda a, kk: np.take_along_axis(a, kk[None], axis=0)[0]
w = (2.0 - take(pv_top, k)) / (take(pv_top, k - 1) - take(pv_top, k))
dt = {name: np.where(found, take(a, k) + w * (take(a, k - 1) - take(a, k)), np.nan)
for name, a in (('theta', th_top), ('u', u_top), ('v', v_top))}
low = slice(gfs3d.LEVELS.index(925), gfs3d.LEVELS.index(850) + 1)
zeta_low = gfs3d.smooth(zeta[low].mean(axis=0), 6) * 1e5 # 10^-5 s^-1
data = ccrs.PlateCarree()
fig = plt.figure(figsize=(11, 7.6), dpi=150)
ax = fig.add_subplot(projection=ccrs.LambertConformal(central_longitude=-115,
central_latitude=45,
standard_parallels=(35, 60)))
ax.set_extent([-170, -62, 22, 68], crs=data)
filled = ax.contourf(lon, lat, dt['theta'], levels=np.arange(282, 375, 4), cmap='turbo',
extend='both', transform=data)
ax.contour(lon, lat, zeta_low, levels=np.arange(8, 41, 4), colors='black', linewidths=0.8,
transform=data)
step = 14
kt = 1.94384
ax.barbs(lon[::step], lat[::step], dt['u'][::step, ::step] * kt, dt['v'][::step, ::step] * kt,
length=5, linewidth=0.45, color='#222222', transform=data)
ax.coastlines(linewidth=0.8, color='#111111')
ax.gridlines(linewidth=0.25, color='#555555', linestyle='--')
fig.colorbar(filled, ax=ax, shrink=0.75, pad=0.02, label='Potential temperature on 2 PVU (K)')
ax.set_title('Dynamic tropopause, 925-850 hPa cyclonic vorticity', loc='left', fontsize=11)
ax.set_title('GFS 06 UTC 27 Oct 2026', loc='right', fontsize=9, color='#555555')
ax.text(0, -0.02, 'Black contours: 925-850 hPa relative vorticity from 8 x 10$^{-5}$ s$^{-1}$, '
'every 4. Barbs: wind on the tropopause (kt).', fontsize=8, color='#444444',
va='top', transform=ax.transAxes)
fig.savefig(os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'images',
'dynamic_tropopause.png'), bbox_inches='tight')
docs/examples/gfs3d.py, shared by the diagnostics
"""Shared by the diagnostic examples: GFS fields on pressure levels, read with
fastgrib, and the finite differences on a sphere the diagnostics are built
from.
The fields come from the GFS 0.5° run of 00 UTC 27 October 2026, six-hour
forecast, in NOAA's archive on AWS. fastgrib reads the file's index and
fetches only the records asked for. Set GFS_FILE_0P50 to read a local copy.
"""
import os
import numpy as np
import fastgrib
URL = 'https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20261027/00/atmos/gfs.t00z.pgrb2.0p50.f006'
LEVELS = (1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350,
300, 250, 200, 150, 100, 70, 50) # hPa, bottom to top
R_EARTH = 6.371e6 # m
G = 9.80665 # m s-2
OMEGA = 7.2921e-5 # s-1
KAPPA = 0.2857 # R_d / c_p
def open_gfs():
return fastgrib.open(os.environ.get('GFS_FILE_0P50', URL))
def levels(grib, name, levels=LEVELS):
"""One variable on the given pressure levels, as (level, lat, lon), bottom
to top. All records are fetched in one call."""
recs = grib.select(name=name, level_kind='isobaricInhPa',
level_value=lambda v: v in set(levels))
recs = sorted(recs, key=lambda r: -r.level_value)
if [r.level_value for r in recs] != list(levels):
raise ValueError('%s: levels %s' % (name, [r.level_value for r in recs]))
data = grib.read(recs).astype(np.float64)
return data[None] if data.ndim == 2 else data # one level: still (1, lat, lon)
def axes(grib):
"""The grid's 1-D latitudes (south to north) and longitudes (east from 0)."""
lat, lon = grib.latlon()
return lat[:, 0].astype(np.float64), lon[0].astype(np.float64)
def ddx(f, lat, lon):
"""d/dx on the sphere along the last axis, periodic in longitude."""
dlam = np.radians(lon[1] - lon[0])
dfdl = (np.roll(f, -1, axis=-1) - np.roll(f, 1, axis=-1)) / (2 * dlam)
coslat = np.maximum(np.cos(np.radians(lat)), 1e-6)[:, None]
return dfdl / (R_EARTH * coslat)
def ddy(f, lat):
"""d/dy on the sphere along the second-last axis."""
return np.gradient(f, np.radians(lat), axis=-2) / R_EARTH
def ddp(f, p_hpa):
"""d/dp along the first axis, p in Pa (levels need not be evenly spaced)."""
return np.gradient(f, np.asarray(p_hpa, np.float64) * 100.0, axis=0)
def coriolis(lat):
return 2 * OMEGA * np.sin(np.radians(lat))[:, None]
def vorticity(u, v, lat, lon):
"""Relative vorticity, with the curvature term of a spherical grid."""
tan = np.tan(np.radians(lat))[:, None]
return ddx(v, lat, lon) - ddy(u, lat) + u * tan / R_EARTH
def theta(t, p_hpa):
"""Potential temperature of temperature t (K) at pressure p (hPa, along axis 0
or a scalar)."""
p = np.asarray(p_hpa, np.float64)
if p.ndim:
p = p[:, None, None]
return t * (1000.0 / p) ** KAPPA
def smooth(f, passes=1):
"""A 9-point smoother on the last two axes, periodic in longitude."""
for _ in range(passes):
g = np.concatenate([f[..., :1, :], f, f[..., -1:, :]], axis=-2)
g = np.concatenate([g[..., -1:], g, g[..., :1]], axis=-1)
sides = g[..., 2:, 1:-1] + g[..., :-2, 1:-1] + g[..., 1:-1, 2:] + g[..., 1:-1, :-2]
corners = g[..., 2:, 2:] + g[..., 2:, :-2] + g[..., :-2, 2:] + g[..., :-2, :-2]
f = (4 * g[..., 1:-1, 1:-1] + 2 * sides + corners) / 16
return f
def bilinear(f, lat, lon, plat, plon):
"""f (..., lat, lon) at points (plat, plon) by bilinear interpolation."""
j = np.interp(plat, lat, np.arange(lat.size))
i = ((np.asarray(plon) - lon[0]) % 360.0) / (lon[1] - lon[0])
j0 = np.clip(np.floor(j).astype(int), 0, lat.size - 2)
i0 = np.floor(i).astype(int) % lon.size
i1 = (i0 + 1) % lon.size
wj, wi = j - j0, i - np.floor(i)
return ((1 - wj) * ((1 - wi) * f[..., j0, i0] + wi * f[..., j0, i1])
+ wj * ((1 - wi) * f[..., j0 + 1, i0] + wi * f[..., j0 + 1, i1]))