"""era5_to_lagranto.py — convert ERA5/MERRA-2 pressure-level data to the netCDF layout LAGRANTO.MPAS `caltra` demands, and build analytic test fields. Contract enforced here (derived from code/lib/ioinp.f, verified 2026-07-28): * filename latlon.current.YYYY-MM-DD_00:00:00.nc (one file per day) * dims AND coord vars named exactly lon, lat, lev, time * variables declared in CDL order (time, lev, lat, lon) * uzonal, umeridional, w in m/s ; height in m (geometric) * exactly 4 times per file at 00/06/12/18 UTC * levels ordered SURFACE-FIRST (pressure descending 1000->1 hPa) so that ioinp's ps = height[raw index 1] really is the topography. Ordering top-first makes ps the model top and every parcel starts "underground". """ import numpy as np from netCDF4 import Dataset RD, G0, RE = 287.05, 9.80665, 6371229.0 FILL = -999.0 def omega_to_w(omega_pa_s, t_kelvin, p_pa): """ERA5 vertical velocity omega [Pa/s] -> w [m/s] (w = -omega/(rho g)).""" rho = p_pa / (RD * t_kelvin) return -omega_pa_s / (rho * G0) def geopotential_to_geometric(z_m2_s2): """ERA5 geopotential [m^2/s^2] -> geometric height [m] above the geoid.""" gph = z_m2_s2 / G0 # geopotential height return RE * gph / (RE - gph) # geometric height def _order_surface_first(levels_hpa, *fields): """Return levels sorted DESCENDING (surface first) plus fields reordered.""" idx = np.argsort(levels_hpa)[::-1] return levels_hpa[idx], [f[:, idx, :, :] for f in fields] def write_caltra_input(path, date_str, lon, lat, levels_hpa, uzonal, umeridional, w, height, extra=None, history=""): """Write one daily caltra input file. date_str : 'YYYY-MM-DD' (goes in the filename; caltra always uses _00:00:00) arrays : shape (4, nlev, nlat, nlon), times 00/06/12/18 UTC extra : {varname: array} traced fields, same shape """ for name, a in [("uzonal", uzonal), ("umeridional", umeridional), ("w", w), ("height", height)]: if a.shape != (4, len(levels_hpa), len(lat), len(lon)): raise ValueError("%s has shape %s, expected %s" % ( name, a.shape, (4, len(levels_hpa), len(lat), len(lon)))) if levels_hpa[0] < levels_hpa[-1]: raise ValueError("levels must be DESCENDING (surface first); see module docstring") ds = Dataset(path, "w", format="NETCDF3_CLASSIC") ds.createDimension("lon", len(lon)) ds.createDimension("lat", len(lat)) ds.createDimension("lev", len(levels_hpa)) ds.createDimension("time", 4) for nm, vals, units in [("lon", lon, "degrees_east"), ("lat", lat, "degrees_north"), ("lev", levels_hpa, "hPa"), ("time", np.array([0, 6, 12, 18], "f4"), "hours since 00 UTC")]: v = ds.createVariable(nm, "f4", (nm,)); v.units = units; v[:] = vals payload = {"uzonal": (uzonal, "m s-1"), "umeridional": (umeridional, "m s-1"), "w": (w, "m s-1"), "height": (height, "m")} payload.update({k: (a, "unknown") for k, a in (extra or {}).items()}) for nm, (arr, units) in payload.items(): v = ds.createVariable(nm, "f4", ("time", "lev", "lat", "lon")) v.units = units v.missing_value = FILL v[:] = arr.astype("f4") ds.date = date_str ds.history = history ds.close() return path def analytic_constant_flow(lon, lat, levels_hpa, u0=20.0, v0=0.0, w0=0.0, t_iso=220.0): """Spatially uniform wind on a US-standard-ish height column. Chosen because the resulting trajectory is EXACT in closed form, so it tests caltra's integrator rather than any meteorology: dlon/dt = u0 / (1.112e5 * cos(lat)) [deg/s] (lat fixed when v0=0) dz/dt = w0 NOT a weather field and not usable for transport science. """ nlev, nlat, nlon = len(levels_hpa), len(lat), len(lon) p = np.asarray(levels_hpa, float) * 100.0 # isothermal hypsometric column -> monotonic height, surface-first ordering z_col = (RD * t_iso / G0) * np.log(p[0] / p) shape = (4, nlev, nlat, nlon) height = np.broadcast_to(z_col[None, :, None, None], shape).copy() return (np.full(shape, u0), np.full(shape, v0), np.full(shape, w0), height, np.full(shape, t_iso)) def verify_caltra_input(path, expect_release_m=None): """Post-write audit of a caltra input file. Returns {ok, errors, warnings, info}. Catches the failure mode verified on 2026-07-28: with levels ordered top-first, ioinp.f takes height[raw index 1] as the topography, so parcels launched at 20 km were silently moved to 44.5 km while caltra reported 0 parcels leaving the domain and exit code 0. Nothing in the output flags it, so the file must be audited before the run. """ from netCDF4 import Dataset errors, warnings, info = [], [], {} ds = Dataset(path) try: for d in ("lon", "lat", "lev", "time"): if d not in ds.dimensions: errors.append("missing dimension '%s'" % d) if d not in ds.variables: errors.append("missing coordinate variable '%s'" % d) for v in ("uzonal", "umeridional", "w", "height"): if v not in ds.variables: errors.append("missing required variable '%s'" % v) elif ds.variables[v].dimensions != ("time", "lev", "lat", "lon"): errors.append("%s dims are %s, ioinp requires ('time','lev','lat','lon')" % (v, ds.variables[v].dimensions)) if len(ds.dimensions.get("time", range(0))) != 4: errors.append("time length is %d, caltra maps exactly 4 records (00/06/12/18)" % len(ds.dimensions["time"])) if errors: return {"ok": False, "errors": errors, "warnings": warnings, "info": info} lev = ds.variables["lev"][:] h = ds.variables["height"][:] info["lev_first"], info["lev_last"] = float(lev[0]), float(lev[-1]) # topography as ioinp computes it: raw level index 1, BEFORE any swap ps = h[0, 0, :, :] info["implied_topography_m"] = [float(ps.min()), float(ps.max())] info["height_range_m"] = [float(h.min()), float(h.max())] if lev[0] < lev[-1]: errors.append( "levels are TOP-FIRST (lev[0]=%g < lev[-1]=%g hPa). ioinp will treat " "height[level 1] = %.0f m as the topography; parcels get silently " "lifted to there and the run still exits 0. Order pressure DESCENDING." % (lev[0], lev[-1], ps.max())) if ps.max() > 2000.0: errors.append("implied topography reaches %.0f m — implausible; level " "ordering is almost certainly inverted" % ps.max()) if expect_release_m is not None: if expect_release_m <= ps.max(): errors.append("release height %.0f m is at/below the implied topography " "%.0f m — parcels will be pinned to the surface" % (expect_release_m, ps.max())) if expect_release_m > h.max(): errors.append("release height %.0f m is above the top of the height " "field (%.0f m)" % (expect_release_m, h.max())) for v in ("uzonal", "umeridional", "w"): a = ds.variables[v][:] lim = 200.0 if v != "w" else 50.0 if abs(a).max() > lim: warnings.append("%s reaches %.1f m/s — check units (w must be m/s, " "not Pa/s)" % (v, abs(a).max())) lon = ds.variables["lon"][:] if len(lon) > 1: dlon = (lon[-1] - lon[0]) / (len(lon) - 1) info["periodic_closear"] = bool(abs((lon[-1] - lon[0] - 360.0) + dlon) < 1e-4) finally: ds.close() return {"ok": not errors, "errors": errors, "warnings": warnings, "info": info}