From 740fc3ab6e790c820c16a0f3ec89002097d1c5b1 Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Thu, 11 Jun 2026 21:49:30 +0000
Subject: [PATCH 1/7] gitignore: ignore WPS_GEOG and geog downloads

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 .gitignore | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/.gitignore b/.gitignore
index 0d1c3f67..fdfe91bd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -72,3 +72,9 @@ met_data
 
 # Python compiled bytecode folders
 __pycache__
+
+# WPS geographical static data (large downloads — never commit)
+WPS_GEOG/
+._WPS_GEOG
+geog_*.tar.gz
+geog_*.tar.bz2

From c24534c1de5f0cf379121ccb1202a4d095f970ad Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Mon, 15 Jun 2026 19:26:14 +0000
Subject: [PATCH 2/7] gitignore: ignore stray HOME/per-user dirs created in the
 repo

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 .gitignore | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/.gitignore b/.gitignore
index fdfe91bd..e86b337b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -78,3 +78,16 @@ WPS_GEOG/
 ._WPS_GEOG
 geog_*.tar.gz
 geog_*.tar.bz2
+
+# Stray per-user / HOME files created in the repo when a tool ran with HOME
+# pointing here (conda cache+envs, shell history, etc.). Never commit; the real
+# ones live in $HOME and these copies are safe to delete.
+.conda/
+.cache/
+.local/
+.config/
+.claude/
+.condarc
+.viminfo
+.python_history
+.bash_history

From e4a586678f100ea2f1bd492212b6327edd267e7a Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Mon, 15 Jun 2026 18:51:07 +0000
Subject: [PATCH 3/7] usp-utils: robust conda env activation + add install
 script shebang

- setup_environment.sh: stop inferring env existence from 'conda activate's exit
  code (which gives a misleading "Couldn't find the environment" when the env
  exists but a transient activate returns non-zero). Now: ensure the 'conda'
  shell function is loaded (source conda.sh if only the binary is on PATH),
  check the env in 'conda env list', then activate. Clearer not-found message.
- install_conda_environment.sh: add '#!/usr/bin/env bash' shebang. It uses
  bash-only features (BASH_SOURCE, '[ == ]'); running it with sh/dash failed
  with "Bad substitution" / "unexpected operator" and left SCRIPT_DIR empty.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 usp-utils/install_conda_environment.sh |  7 ++++++-
 usp-utils/setup_environment.sh         | 26 ++++++++++++++++----------
 2 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/usp-utils/install_conda_environment.sh b/usp-utils/install_conda_environment.sh
index 3d3e4004..eeb9bca3 100755
--- a/usp-utils/install_conda_environment.sh
+++ b/usp-utils/install_conda_environment.sh
@@ -1,4 +1,9 @@
-#get directory where the sourced script is located
+#!/usr/bin/env bash
+# Run with bash (./install_conda_environment.sh or bash install_conda_environment.sh).
+# It uses bash-only features (BASH_SOURCE, [ == ]); running it with `sh`/dash
+# fails with "Bad substitution" / "unexpected operator".
+
+# get directory where this script is located
 SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
 
 CYAN='\033[0;36m'
diff --git a/usp-utils/setup_environment.sh b/usp-utils/setup_environment.sh
index d1b6ba20..5cf0aee7 100644
--- a/usp-utils/setup_environment.sh
+++ b/usp-utils/setup_environment.sh
@@ -29,18 +29,24 @@ fi
 export PYTHONPATH="$SCRIPT_DIR/libs/py:$PYTHONPATH"
 echo -e "${INFO} New enviroment variable set: PYTHONPATH=$PYTHONPATH"
 
-conda &> /dev/null
-status=$?
-if [ ! $status == "0" ]; then
-    echo -e "${WARNING} Couldn't find 'conda' binary, 'cgfd-usp-mpas' conda environment will not be activated. Python scripts might not work."
+if ! command -v conda &> /dev/null; then
+    echo -e "${WARNING} Couldn't find 'conda', the 'cgfd-usp-mpas' conda environment will not be activated. Python scripts might not work."
     return 1
 fi
 
-conda activate cgfd-usp-mpas &> /dev/null
-status=$?
-if [ ! $status == "0" ]; then
-    echo -e "${WARNING} Couldn't find the 'cgfd-usp-mpas' conda environment. It will not be activated and Python scripts might not work. If you wish to install the conda environment please run the script in $SCRIPT_DIR/install_conda_environment.sh"
-    return 1
+# 'conda activate' needs the 'conda' shell function; load it if only the binary
+# is on PATH (e.g. in a shell where 'conda init' hasn't run).
+if ! type conda 2> /dev/null | grep -q 'function'; then
+    __conda_sh="$(conda info --base 2> /dev/null)/etc/profile.d/conda.sh"
+    [ -f "$__conda_sh" ] && source "$__conda_sh"
 fi
 
-echo -e "${INFO} Conda enviroment 'cgfd-usp-mpas' activated"
+# Decide from the actual list of environments, not from 'activate's exit code
+# (the env may exist even if a transient activate call returns non-zero).
+if conda env list 2> /dev/null | awk '{print $1}' | grep -qx 'cgfd-usp-mpas'; then
+    conda activate cgfd-usp-mpas
+    echo -e "${INFO} Conda enviroment 'cgfd-usp-mpas' activated"
+else
+    echo -e "${WARNING} The 'cgfd-usp-mpas' conda environment was not found. It will not be activated and Python scripts might not work. To install it, run: bash $SCRIPT_DIR/install_conda_environment.sh"
+    return 1
+fi

From 62be7ce74e74880d35508b63431c9fbddc9f6864 Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Sat, 20 Jun 2026 17:00:30 +0000
Subject: [PATCH 4/7] pre_proc: real-data + static-field preparation pipelines
 and docs
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Add the pre-processing utilities for MPAS init_atmosphere under usp-utils/pre_proc/.

real_data/ — download + intermediate-file generation for real-data init:
- download_{era5,gfs,oisst}.py and {era5,gfs,gfs_sst,oisst,oisst_clim}_to_intermediate.py
- prepare_{era5,gfs,oisst}.sh end-to-end wrappers
- README plus per-scenario recipes (operational/hindcast x global/regional)

static_fields/ — static-field generation (init_atmosphere config_init_case = 7):
- download_static_data.sh (mesh + MPAS geog bundle, integrity-checked)
- STATIC_FIELDS_GUIDE.md walkthrough

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 usp-utils/pre_proc/real_data/README.md        | 296 +++++++++++++
 usp-utils/pre_proc/real_data/download_era5.py | 166 ++++++++
 usp-utils/pre_proc/real_data/download_gfs.py  | 391 ++++++++++++++++++
 .../pre_proc/real_data/download_oisst.py      | 132 ++++++
 .../real_data/era5_to_intermediate.py         | 176 ++++++++
 .../real_data/gfs_sst_to_intermediate.py      |  98 +++++
 .../pre_proc/real_data/gfs_to_intermediate.py | 128 ++++++
 .../real_data/oisst_clim_to_intermediate.py   | 126 ++++++
 .../real_data/oisst_to_intermediate.py        | 115 ++++++
 usp-utils/pre_proc/real_data/prepare_era5.sh  |  55 +++
 usp-utils/pre_proc/real_data/prepare_gfs.sh   |  81 ++++
 usp-utils/pre_proc/real_data/prepare_oisst.sh |  55 +++
 .../pre_proc/real_data/recipes/README.md      |  25 ++
 .../real_data/recipes/hindcast_global.md      |  36 ++
 .../real_data/recipes/hindcast_regional.md    |  45 ++
 .../real_data/recipes/operational_global.md   |  41 ++
 .../real_data/recipes/operational_regional.md |  41 ++
 .../static_fields/STATIC_FIELDS_GUIDE.md      | 224 ++++++++++
 .../static_fields/download_static_data.sh     | 254 ++++++++++++
 19 files changed, 2485 insertions(+)
 create mode 100644 usp-utils/pre_proc/real_data/README.md
 create mode 100755 usp-utils/pre_proc/real_data/download_era5.py
 create mode 100755 usp-utils/pre_proc/real_data/download_gfs.py
 create mode 100755 usp-utils/pre_proc/real_data/download_oisst.py
 create mode 100755 usp-utils/pre_proc/real_data/era5_to_intermediate.py
 create mode 100755 usp-utils/pre_proc/real_data/gfs_sst_to_intermediate.py
 create mode 100755 usp-utils/pre_proc/real_data/gfs_to_intermediate.py
 create mode 100755 usp-utils/pre_proc/real_data/oisst_clim_to_intermediate.py
 create mode 100755 usp-utils/pre_proc/real_data/oisst_to_intermediate.py
 create mode 100755 usp-utils/pre_proc/real_data/prepare_era5.sh
 create mode 100755 usp-utils/pre_proc/real_data/prepare_gfs.sh
 create mode 100755 usp-utils/pre_proc/real_data/prepare_oisst.sh
 create mode 100644 usp-utils/pre_proc/real_data/recipes/README.md
 create mode 100644 usp-utils/pre_proc/real_data/recipes/hindcast_global.md
 create mode 100644 usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
 create mode 100644 usp-utils/pre_proc/real_data/recipes/operational_global.md
 create mode 100644 usp-utils/pre_proc/real_data/recipes/operational_regional.md
 create mode 100644 usp-utils/pre_proc/static_fields/STATIC_FIELDS_GUIDE.md
 create mode 100755 usp-utils/pre_proc/static_fields/download_static_data.sh

diff --git a/usp-utils/pre_proc/real_data/README.md b/usp-utils/pre_proc/real_data/README.md
new file mode 100644
index 00000000..7b8781b5
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/README.md
@@ -0,0 +1,296 @@
+# `real_data/` — GFS / ERA5 → WPS intermediate met files for MPAS
+
+Scripts to fetch meteorological data and turn it into the **WPS intermediate-format**
+file (`GFS:YYYY-MM-DD_HH`) that `init_atmosphere` reads to build the real-data
+`init.nc`. MPAS has no built-in ungrib, so this is done in Python with **pywinter**
+— no need to build WPS/ungrib.
+
+```
+download_gfs.py   →  GFS GRIB2    →  gfs_to_intermediate.py   ┐
+                                                              ├→  GFS:YYYY-MM-DD_HH
+download_era5.py  →  ERA5 GRIB    →  era5_to_intermediate.py  ┘
+```
+
+The intermediate file keeps the `GFS:` prefix (= `config_met_prefix='GFS'`)
+regardless of the source, so the `init_atmosphere` setup is identical either way.
+This directory only produces that **met intermediate file** — running
+`init_atmosphere` with it (namelist, streams, linking `*.static.nc`, etc.) is covered
+by the model usage tutorial, not repeated here. Run these after the static fields
+are ready (see `../static_fields/`).
+
+## Which source? (data availability)
+
+Two products feed init_atmosphere: the **atmosphere** (3D state → `GFS:` intermediate,
+init case 7) and the **SST/sea-ice** lower boundary (→ `SST:` intermediate, init case 8,
+see "SST / sea-ice update").
+
+**Atmosphere:**
+
+| Source | Coverage | Resolution | Notes |
+|---|---|---|---|
+| GFS `--source aws` (default) | **2021-03-23 → present** | 0.25° | S3 byte-range subset, no rate limit |
+| GFS `--source nomads` | **last ~10 days** | 0.25/0.5/1° | grib_filter; can hit `Over Rate Limit` |
+| GFS `--source rda` | **2015-01-15 → present** | 0.25° | NCAR RDA ds084.1, full files; may need a free RDA account |
+| ERA5 (`download_era5.py`) | **1940 → present** | 0.25° | Copernicus CDS; needs a CDS account |
+
+**SST / sea-ice:**
+
+| Source | Coverage | Notes |
+|---|---|---|
+| GFS (`gfs_sst_to_intermediate.py`, or `--fields sst`) | same as GFS above | from the GFS GRIB; analysis **or** forecast |
+| OISST observed (`download_oisst.py`) | **1981-09-01 → present** | NOAA OI SST v2.1 daily 0.25° (NCEI) |
+| OISST climatology (`oisst_clim_to_intermediate.py`) | **any date** (LTM 1991-2020) | day-of-year climatology via OPeNDAP; for forecasts / typical-month runs |
+
+### GFS is a *forecast* model — analysis vs forecast
+GFS runs four times a day, at **00, 06, 12 and 18 UTC**. Each of those runs is a
+**cycle**: `--cycle 00` selects the 00 UTC run of `--date`. Within a cycle, `--fhour`
+is the **forecast hour** (lead time, in hours after the cycle time).
+
+A cycle produces an **analysis** at forecast hour 0 (`--fhour 0`, the default) and
+**forecasts** at later hours (`--fhour 3, 6, … up to 384`, i.e. 16 days). So:
+- **Analysis** (`--fhour 0`) = best estimate of the atmosphere *at the cycle time*. Use it
+  for the initial conditions. Only exists for the past (and the current cycle).
+- **Forecast** (`--fhour > 0`) = a prediction for times *after* the cycle, i.e. it can give
+  data for the **future** (times that have not happened yet) — that is what makes an
+  operational MPAS forecast possible.
+
+By contrast **ERA5 (reanalysis) and OISST observed only exist for the past.** For the
+future (or any date with no observed SST) use the **OISST climatology**.
+
+There is **no GFS 0.25° before 2015-01-15** — for anything older, use ERA5 (atmosphere)
+and OISST (SST). `download_gfs.py` fails fast with the exact limit for out-of-range dates.
+
+## Scripts
+
+| Script | What it does |
+|--------|--------------|
+| `download_gfs.py` | Downloads the GFS fields/levels MPAS needs, with a `*.provenance.json`. Sources `aws` / `nomads` / `rda` (see table above). `--fields sst` for a small surface-only (SST update) download. |
+| `gfs_to_intermediate.py` | Converts one GFS GRIB2 file to a WPS intermediate file via pywinter. Flips latitude to S→N, writes pressure in **Pa**, prints the level count for `config_nfglevels`. |
+| `download_era5.py` | Downloads ERA5 pressure- and single-level GRIB via the CDS API (cdsapi), with provenance. For dates GFS can't reach. |
+| `era5_to_intermediate.py` | ERA5 counterpart of the GFS converter (same intermediate field names; handles geopotential→height and snow units). |
+| `download_oisst.py` | Downloads NOAA OISST v2.1 daily 0.25° SST/sea-ice NetCDF (NCEI) over a date range. For SST update files (scientific pipeline). |
+| `oisst_to_intermediate.py` | Converts OISST daily files to `SST:` intermediate files (SST→K, SEAICE, LANDSEA) — one per day. |
+| `oisst_clim_to_intermediate.py` | Builds `SST:` files from the OISST daily **climatology** (LTM 1991-2020, via OPeNDAP) for target dates with no observed SST (forecasts / typical-month runs). |
+| `gfs_sst_to_intermediate.py` | Extracts SST/SEAICE/LANDSEA from a GFS GRIB into an `SST:` intermediate file (operational SST update). |
+
+One-shot wrappers (download + convert in one call, inside the conda env):
+
+| Wrapper | Covers |
+|---------|--------|
+| `prepare_gfs.sh` | GFS atmosphere and/or SST — `--product atm` (default), `sst`, or `both` (one download, both products). |
+| `prepare_era5.sh` | ERA5 atmosphere (download pl+sl → `GFS:` intermediate). |
+| `prepare_oisst.sh` | OISST SST over a date range — observed (default) or `--climatology` (→ daily `SST:` intermediates). |
+
+## Dependencies
+
+The `cgfd-usp-mpas` conda env: `cfgrib` + `eccodes` (conda-forge) and `pywinter`
+(pip); plus `cdsapi` (pip) for the ERA5 path. Run with
+`conda run -n cgfd-usp-mpas python <script>.py ...`, or use the `prepare_*.sh` wrappers.
+ERA5 also needs a CDS account and `~/.cdsapirc`
+(https://cds.climate.copernicus.eu/how-to-api); RDA may need a free account at
+https://rda.ucar.edu (set `RDA_API_TOKEN`). Never hard-code credentials.
+
+## Recipes (by use case)
+
+Step-by-step data-prep recipes live in [`recipes/`](recipes/). Pick by **what you
+simulate** (future vs past) and **mesh type** (global vs regional — regional needs
+lateral boundary data):
+
+| | Global mesh | Regional mesh (LBCs) |
+|---|---|---|
+| **Operational forecast** (future) | [operational_global](recipes/operational_global.md) | [operational_regional](recipes/operational_regional.md) |
+| **Hindcast / case study** (past) | [hindcast_global](recipes/hindcast_global.md) | [hindcast_regional](recipes/hindcast_regional.md) |
+
+In one line each:
+- **Operational** → atmosphere from GFS **analysis**; future SST from GFS **forecast**
+  hours or OISST **climatology** (the future has no observed SST).
+- **Hindcast** → atmosphere from **ERA5** (or recent GFS analysis); SST from **OISST
+  observed** daily.
+- **Regional** adds lateral-boundary data at a cadence (GFS forecast hours, or ERA5 6-hourly).
+
+Minimal examples:
+```sh
+./prepare_gfs.sh   --date 2026-06-16 --cycle 00 --product both   # operational, global
+./prepare_era5.sh  --date 2014-09-10 --time 00                   # hindcast atmosphere
+./prepare_oisst.sh --start 2014-09-10 --end 2014-09-15           # hindcast SST (daily)
+```
+
+> **SST for an N-day forecast:** `--product both` gives SST only at the *initial* time.
+> Short forecasts usually hold SST fixed (`config_sst_update=.false.`). To evolve it,
+> loop GFS forecast hours (`--product sst --fhour 024 048 …`, each tagged with its valid
+> date) or use the OISST climatology — see [operational_global](recipes/operational_global.md).
+
+Calling the Python scripts directly (instead of the wrappers) lets you inspect each stage;
+every `download_*.py` prints the exact next command. Output goes under `<repo>/met_data/`
+(`gfs/`, `era5/`, `oisst/`, `oisst_clim/`) by default; `--outdir` changes it.
+
+Key options:
+- `download_gfs.py`: `--date YYYY-MM-DD`, `--cycle {00,06,12,18}` (UTC run hour),
+  `--fhour` (forecast lead time in hours; 0 = analysis), `--res {0p25,0p50,1p00}`,
+  `--source {aws,nomads,rda}`,
+  `--fields {full,sst}` (`sst` = small surface-only download for SST update files),
+  `--outdir`.
+- `download_era5.py`: `--date YYYY-MM-DD`, `--time HH`, `--area N W S E` (default
+  global), `--outdir`.
+- `download_oisst.py` / `oisst_to_intermediate.py`: `--start`/`--end` date range,
+  `--hour HH`, `--outdir`.
+- `oisst_clim_to_intermediate.py` (climatology): `--start`/`--end` target dates,
+  `--hour HH`, `--outdir` (no download — reads the LTM over OPeNDAP).
+- converters: `--outdir`, `--prefix` (default `GFS` for atmosphere, `SST` for SST);
+  GFS uses `--grib`, ERA5 uses `--pl`/`--sl`.
+
+The converter sets the date tag from the data's valid time, so the output name
+(`GFS:YYYY-MM-DD_HH`) is what `config_start_time` must match.
+
+## Verify the init.nc (recommended after every run)
+
+A clean `init_atmosphere` run does **not** guarantee a usable initial state — it can
+finish with "0 errors" and still write a physically broken `init.nc` (e.g. from a
+units/interpolation problem upstream), which then blows up the model (`w → NaN`) on
+the first step. After building the `init.nc`, sanity-check the **prognostic** state:
+
+```python
+import numpy as np, xarray as xr
+ds = xr.open_dataset("x1.10242.init.nc")
+Rd, cp = 287.0, 1004.5; cv = cp - Rd; gamma = cp/cv; p0 = 1e5
+th, rho = ds["theta"].values[0,:,0], ds["rho"].values[0,:,0]   # surface level
+print("theta surf:", th.min(), th.max())   # expect ~270–320 K  (NOT ~1100)
+print("rho   surf:", rho.min(), rho.max())  # expect ~1.0–1.3    (NOT ~0.01)
+p = p0*(Rd*rho*th/p0)**gamma / 100.0        # EOS pressure from prognostic state
+print("EOS(prognostic) hPa:", p.mean())     # expect ~950–1010 hPa (NOT ~10)
+```
+`rho_base`/`theta_base` and `surface_pressure` can look fine even when the
+prognostic `theta`/`rho` are broken, so check `theta`/`rho` specifically.
+
+## Field mapping (GRIB → WPS intermediate)
+
+Both converters emit the same intermediate field names (WPS GFS Vtable):
+
+| Field | Intermediate | Kind | GFS source | ERA5 source |
+|-------|--------------|------|------------|-------------|
+| geopotential height | `GHT`     | 3D / pressure | HGT       | z / g0 |
+| temperature         | `TT`      | 3D / pressure | TMP       | t |
+| wind u / v          | `UU`/`VV` | 3D / pressure | UGRD/VGRD | u / v |
+| relative humidity   | `RH`      | 3D / pressure | RH        | r |
+| mean sea level pres.| `PMSL`    | 2D            | PRMSL     | msl |
+| surface pressure    | `PSFC`    | 2D            | PRES sfc  | sp |
+| terrain height      | `SOILHGT` | 2D            | HGT sfc   | z(sfc) / g0 |
+| skin temperature    | `SKINTEMP`| 2D            | TMP sfc   | skt |
+| land-sea mask       | `LANDSEA` | 2D            | LAND      | lsm |
+| sea-ice fraction    | `SEAICE`  | 2D            | ICEC      | siconc |
+| snow water equiv.   | `SNOW`    | 2D            | WEASD     | sd × 1000 |
+| soil temperature    | `ST<tag>` | 2D (4 layers) | TSOIL     | stl1–4 |
+| soil moisture       | `SM<tag>` | 2D (4 layers) | SOILW     | swvl1–4 |
+
+Soil layer tags are the source's actual layer depths (GFS: `000010/010040/040100/
+100200`; ERA5: `000007/007028/028100/100289`, cm). All fields share the global
+regular lat/lon grid; pywinter takes the geometry once (`Geo0`).
+
+## SST / sea-ice update files
+
+For runs longer than a day or two you usually want the lower boundary (SST and
+sea ice) to evolve in time instead of staying fixed at the initial value. MPAS
+supports this with a **surface update file** (`x1.<mesh>.sfc_update.nc`), built
+independently of the atmospheric initial conditions. Two sources here:
+
+- **Operational** (GFS atmosphere): `gfs_sst_to_intermediate.py` — reuses the GFS
+  GRIB you already downloaded (SST = surface skin temperature, SEAICE = ICEC,
+  LANDSEA = LAND).
+- **Scientific / downscaling** (ERA5 atmosphere): `download_oisst.py` +
+  `oisst_to_intermediate.py` — NOAA OISST v2.1 daily 0.25° (an independent,
+  reference daily SST that matches the 0.25° atmosphere resolution).
+
+The whole flow has three steps:
+
+**1. Build `SST:YYYY-MM-DD_HH` intermediate files** (fields `SST`, `SEAICE`,
+`LANDSEA`), one per update time:
+```sh
+# operational (GFS), SST only — lightweight: downloads just the surface
+# SST/sea-ice/land fields (~1 MB, not the full GFS), then converts:
+./prepare_gfs.sh --date 2026-06-11 --product sst
+#   by hand:
+python download_gfs.py --date 2026-06-11 --fields sst        # -> *.f000.sst (~1 MB)
+python gfs_sst_to_intermediate.py --grib <the .f000.sst file>
+# (if you already downloaded the full GFS for the atmosphere, just convert it:
+#  python gfs_sst_to_intermediate.py --grib <the .f000 file>)
+
+# scientific (OISST): daily files over the run period
+./prepare_oisst.sh --start 2014-09-01 --end 2014-09-30
+#   by hand:
+python download_oisst.py --start 2014-09-01 --end 2014-09-30
+python oisst_to_intermediate.py --indir <repo>/met_data/oisst --start 2014-09-01 --end 2014-09-30
+```
+> `download_gfs.py --fields sst` writes a distinct `*.f000.sst` file (only TMP/ICEC/
+> LAND at the surface), so it never collides with a full atmosphere download in the
+> same directory. `--product both` downloads the full GFS once and makes both
+> `GFS:` and `SST:`.
+
+**2. Run `init_atmosphere` in case 8** to interpolate them onto the mesh and write
+`x1.<mesh>.sfc_update.nc`. Key `namelist.init_atmosphere` settings (static fields
+already done; reuse the same mesh `*.static.nc`):
+```
+&nhyd_model       config_init_case = 8
+                  config_start_time = 'YYYY-MM-DD_HH:00:00'
+                  config_stop_time  = 'YYYY-MM-DD_HH:00:00'   ! end of the run
+&dimensions       config_nsoillevels = 4
+&data_sources     config_sfc_prefix = 'SST'
+                  config_fg_interval = 86400                  ! update cadence (s); 86400 = daily
+&preproc_stages   config_static_interp = .false.
+                  config_vertical_grid = .false.
+                  config_met_interp    = .false.
+                  config_input_sst     = .true.               ! REQUIRED for case 8
+                  config_frac_seaice   = .true.
+```
+The `surface` output stream writes `x1.<mesh>.sfc_update.nc` at `output_interval`
+equal to the cadence; its input stream reads the `*.static.nc`.
+
+**3. Enable the update in the model run** (`namelist.atmosphere`):
+```
+&physics          config_sst_update = .true.
+```
+and add an input stream in `streams.atmosphere` reading the update file at the same
+cadence:
+```xml
+<stream name="surface" type="input" filename_template="x1.<mesh>.sfc_update.nc"
+        input_interval="86400" >
+    <var name="sst"/>
+    <var name="xice"/>
+</stream>
+```
+
+> Notes. `SST` is written in **Kelvin** (OISST is degC → +273.15; GFS skin temp is
+> already K). `LANDSEA` is 1 = land / 0 = water; SST over land is a filler and is
+> masked out via `LANDSEA`. OISST starts 1981-09-01; very recent days come as
+> `_preliminary` files (handled automatically). OISST 0.25° matches the GFS/ERA5
+> atmosphere resolution — finer GHRSST products (e.g. MUR 0.01°) are overkill for a
+> 120–240 km mesh.
+
+## Known issues / gotchas
+
+- **GFS availability — you can't go arbitrarily far back.** GFS 0.25° only exists
+  from **2015-01-15**. `--source aws` covers 2021-03-23+ (when the bucket adopted the
+  `gfs.*/atmos/` layout), `--source rda` covers 2015-01-15+ (NCAR RDA ds084.1, full
+  files), `--source nomads` only the last ~10 days. `download_gfs.py` checks the date
+  against the chosen source and **errors early stating the limit** — a rejected old
+  date is expected, not a bug. For pre-2015 use ERA5.
+- **ERA5 path is newer and not yet validated end-to-end here.** `download_era5.py`
+  and `era5_to_intermediate.py` mirror the validated GFS path, but the ERA5 branch
+  hasn't been run through `init_atmosphere` + model in this repo yet. Always run the
+  *Verify the init.nc* check. Watch for: geopotential is converted to height (÷ g0),
+  snow depth (m w.e.) to kg m⁻² (× 1000), and pressure written in Pa.
+- **Credentials.** ERA5 needs `~/.cdsapirc` (CDS account); RDA may need a free account
+  and `RDA_API_TOKEN`. If RDA returns a login/HTML page instead of GRIB, the script
+  says so — set the token.
+- **NOMADS rate limit.** `--source nomads` can return `Over Rate Limit` (temporary IP
+  block). Prefer `--source aws`; if hit, wait ~1 h.
+- **Pressure units must be Pa.** The WPS intermediate format (and pywinter's `V3dp`)
+  expect pressure in **Pa**. Levels arrive from cfgrib as `isobaricInhPa` (hPa), so the
+  converters multiply by 100. In hPa, every 3D level would be 100× too low:
+  `init_atmosphere` still runs, but `init.nc` gets `theta`/`rho` ~100×/~4× off and the
+  model goes `w → NaN`; `surface_pressure` still looks fine (it comes from PSFC/PMSL,
+  already in Pa), which makes it sneaky — hence the *Verify the init.nc* check.
+- **`config_ztop` vs level coverage.** `config_nfglevels` must equal the value the
+  converter prints (isobaric levels + 1). The downloaders ship levels up to 1 hPa
+  (~48 km), so `config_ztop=30000` (30 km) works. A dataset that stops lower combined
+  with a higher top triggers
+  `extrap_type == 2 not implemented for target_z >= zf(1,nz)` — lower `config_ztop`.
diff --git a/usp-utils/pre_proc/real_data/download_era5.py b/usp-utils/pre_proc/real_data/download_era5.py
new file mode 100755
index 00000000..4f17fa35
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/download_era5.py
@@ -0,0 +1,166 @@
+#!/usr/bin/env python3
+"""
+download_era5.py — fetch ERA5 reanalysis fields needed to initialize MPAS, for
+dates the GFS archive does not cover (GFS 0.25 deg only exists from 2015-01-15;
+ERA5 goes back to 1940).
+
+Downloads two GRIB files via the CDS API (cdsapi):
+  - pressure levels (reanalysis-era5-pressure-levels): z, t, u, v, r
+  - single levels   (reanalysis-era5-single-levels):   sp, msl, skt, z (orog),
+                     lsm, siconc, sd (snow), stl1-4 (soil T), swvl1-4 (soil moist)
+They are then turned into the WPS intermediate format by era5_to_intermediate.py
+(consumed by init_atmosphere with config_met_prefix='GFS', same as the GFS path).
+
+Requirements (in the cgfd-usp-mpas env):
+  - cdsapi (pip)
+  - a CDS account + ~/.cdsapirc with your API key
+    (https://cds.climate.copernicus.eu/how-to-api). Never hard-code the key.
+
+Usage:
+    python download_era5.py --date 2014-09-10 --time 00 [--area N W S E] [--outdir <dir>]
+
+Note: ERA5 is a reanalysis (analysis only, hourly) — there is no forecast hour.
+"""
+
+import argparse
+import json
+import sys
+from datetime import datetime, timezone
+from pathlib import Path
+
+# --- Fields MPAS needs (CDS variable names) -----------------------------------
+# Pressure-level fields (-> 3D intermediate). NOTE: ERA5 'geopotential' is in
+# m^2/s^2; the converter divides by g0 to get geopotential height (GHT, m).
+PL_VARIABLES = [
+    "geopotential",            # z   -> GHT (after /g0)
+    "temperature",             # t   -> TT
+    "u_component_of_wind",     # u   -> UU
+    "v_component_of_wind",     # v   -> VV
+    "relative_humidity",       # r   -> RH
+]
+
+# Single-level fields (-> 2D intermediate).
+SL_VARIABLES = [
+    "surface_pressure",            # sp     -> PSFC
+    "mean_sea_level_pressure",     # msl    -> PMSL
+    "skin_temperature",            # skt    -> SKINTEMP
+    "geopotential",                # z(sfc) -> SOILHGT (after /g0)
+    "land_sea_mask",               # lsm    -> LANDSEA
+    "sea_ice_cover",               # siconc -> SEAICE
+    "snow_depth",                  # sd     -> SNOW (m w.e.; converter -> kg m-2)
+    "soil_temperature_level_1",    # stl1   -> ST000007
+    "soil_temperature_level_2",    # stl2   -> ST007028
+    "soil_temperature_level_3",    # stl3   -> ST028100
+    "soil_temperature_level_4",    # stl4   -> ST100289
+    "volumetric_soil_water_layer_1",  # swvl1 -> SM000007
+    "volumetric_soil_water_layer_2",  # swvl2 -> SM007028
+    "volumetric_soil_water_layer_3",  # swvl3 -> SM028100
+    "volumetric_soil_water_layer_4",  # swvl4 -> SM100289
+]
+
+# Full ERA5 isobaric set (hPa); like the GFS path, includes the upper levels so a
+# 30 km model top has data to interpolate from.
+PRESSURE_LEVELS = [
+    1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550,
+    500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20,
+    10, 7, 5, 3, 2, 1,
+]
+
+
+def _retrieve(client, dataset, request, target):
+    print(f"[INFO] CDS retrieve {dataset} -> {target}")
+    client.retrieve(dataset, request, str(target))
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="Download ERA5 fields for MPAS init.")
+    ap.add_argument("--date", required=True, help="Analysis date, YYYY-MM-DD (UTC)")
+    ap.add_argument("--time", default="00", help="Analysis hour HH (UTC), e.g. 00")
+    ap.add_argument("--area", nargs=4, type=float, metavar=("N", "W", "S", "E"),
+                    default=None, help="Optional bounding box (default: global)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: $MPAS_ROOT/met_data/era5/<date><HH>)")
+    args = ap.parse_args()
+
+    try:
+        d = datetime.strptime(args.date, "%Y-%m-%d").date()
+    except ValueError:
+        print("[ERROR] --date must be YYYY-MM-DD", file=sys.stderr)
+        return 2
+    hh = args.time.zfill(2)
+
+    try:
+        import cdsapi
+    except ImportError:
+        print("[ERROR] cdsapi not installed. In the cgfd-usp-mpas env: pip install cdsapi, "
+              "and set up ~/.cdsapirc (https://cds.climate.copernicus.eu/how-to-api).",
+              file=sys.stderr)
+        return 1
+
+    ymd = args.date.replace("-", "")
+    repo_root = Path(__file__).resolve().parents[3]
+    outdir = (Path(args.outdir) if args.outdir
+              else repo_root / "met_data" / "era5" / f"{ymd}{hh}")
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    common = {
+        "product_type": "reanalysis",
+        "year": [d.strftime("%Y")],
+        "month": [d.strftime("%m")],
+        "day": [d.strftime("%d")],
+        "time": [f"{hh}:00"],
+        "data_format": "grib",          # older cdsapi: use "format" instead
+        "download_format": "unarchived",
+    }
+    if args.area:
+        common["area"] = [args.area[0], args.area[1], args.area[2], args.area[3]]
+
+    pl_file = outdir / f"era5_pl_{ymd}{hh}.grib"
+    sl_file = outdir / f"era5_sl_{ymd}{hh}.grib"
+
+    try:
+        client = cdsapi.Client()
+    except Exception as e:  # noqa: BLE001
+        print(f"[ERROR] cdsapi client init failed (missing ~/.cdsapirc?): {e}",
+              file=sys.stderr)
+        return 1
+
+    try:
+        if not (pl_file.exists() and pl_file.stat().st_size > 0):
+            _retrieve(client, "reanalysis-era5-pressure-levels",
+                      {**common, "variable": PL_VARIABLES,
+                       "pressure_level": [str(p) for p in PRESSURE_LEVELS]}, pl_file)
+        else:
+            print(f"[OK] {pl_file.name} already present — skipping")
+        if not (sl_file.exists() and sl_file.stat().st_size > 0):
+            _retrieve(client, "reanalysis-era5-single-levels",
+                      {**common, "variable": SL_VARIABLES}, sl_file)
+        else:
+            print(f"[OK] {sl_file.name} already present — skipping")
+    except Exception as e:  # noqa: BLE001
+        print(f"[ERROR] CDS retrieve failed: {e}", file=sys.stderr)
+        return 1
+
+    # --- Provenance ----------------------------------------------------------
+    meta = {
+        "dataset": "ERA5",
+        "products": ["reanalysis-era5-pressure-levels",
+                     "reanalysis-era5-single-levels"],
+        "source": "Copernicus CDS (cdsapi)",
+        "valid_time": f"{args.date} {hh}:00 UTC",
+        "pressure_variables": PL_VARIABLES,
+        "single_variables": SL_VARIABLES,
+        "pressure_levels_hpa": PRESSURE_LEVELS,
+        "area": common.get("area", "global"),
+        "files": [str(pl_file), str(sl_file)],
+        "download_date": datetime.now(timezone.utc).isoformat(),
+    }
+    (outdir / f"era5_{ymd}{hh}.provenance.json").write_text(json.dumps(meta, indent=2))
+    print(f"[OK] ERA5 files in {outdir}")
+    print(f"[INFO] Next: python era5_to_intermediate.py "
+          f"--pl {pl_file} --sl {sl_file}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/download_gfs.py b/usp-utils/pre_proc/real_data/download_gfs.py
new file mode 100755
index 00000000..28ec13f5
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/download_gfs.py
@@ -0,0 +1,391 @@
+#!/usr/bin/env python3
+"""
+download_gfs.py — fetch NCEP GFS (0.25 deg) fields needed to initialize MPAS.
+
+Downloads only the variables/levels MPAS init_atmosphere needs (when the source
+supports subsetting) and logs provenance. The downloaded GRIB2 is later turned
+into the WPS intermediate format (consumed by init_atmosphere with
+config_met_prefix='GFS') by gfs_to_intermediate.py.
+
+Sources (--source) and their temporal coverage — GFS 0.25 deg:
+  aws    : AWS Open Data (noaa-gfs-bdp-pds), byte-range subset via the .idx.
+           No rate limit. Covers 2021-03-23 -> present (this script's
+           gfs.YYYYMMDD/HH/atmos/ layout). DEFAULT.
+  nomads : NOMADS grib_filter subset. Keeps only the most recent ~10 days; can
+           hit 'Over Rate Limit'.
+  rda    : NCAR RDA ds084.1 historical archive. Covers 2015-01-15 -> present
+           (the earliest GFS 0.25 deg). Full files (no subset); may need a free
+           RDA account (set RDA_API_TOKEN).
+
+There is no GFS 0.25 deg before 2015-01-15. For older cases use a reanalysis
+(see download_era5.py — ERA5 goes back to 1940). Out-of-range dates fail fast
+with a message stating the limit.
+
+Usage:
+    python download_gfs.py --date 2024-09-10 --cycle 00 [--fhour 0]
+                           [--res 0p25] [--source aws|nomads|rda] [--outdir <dir>]
+
+Notes:
+- GFS analysis = forecast hour 0 (default). For a forecast lead time, set --fhour.
+"""
+
+import argparse
+import json
+import os
+import shutil
+import sys
+import urllib.parse
+import urllib.request
+from datetime import date, datetime, timezone
+from pathlib import Path
+
+# --- Temporal coverage of each source (GFS 0.25 deg) --------------------------
+# AWS has this script's gfs.*/atmos/ layout from 2021-03-23 on; NCAR RDA ds084.1
+# is the GFS 0.25 deg historical archive from 2015-01-15 (nothing earlier exists
+# at 0.25 deg); NOMADS keeps only the most recent ~10 days.
+AWS_START = date(2021, 3, 23)
+RDA_START = date(2015, 1, 15)
+NOMADS_DAYS = 10
+
+# --- Fields MPAS needs from GFS (NOMADS grib_filter variable names) -----------
+# 3D (on isobaric levels) + 2D surface / soil fields.
+VARIABLES = [
+    "HGT",    # geopotential height (-> GHT) ; also surface terrain
+    "TMP",    # temperature (-> TT) ; also skin/2 m
+    "UGRD",   # u-wind (-> UU) ; also 10 m
+    "VGRD",   # v-wind (-> VV) ; also 10 m
+    "RH",     # relative humidity (-> RH)
+    "PRMSL",  # mean sea level pressure (-> PMSL)
+    "PRES",   # surface pressure (-> PSFC)
+    "LAND",   # land-sea mask (-> LANDSEA)
+    "ICEC",   # sea-ice cover (-> SEAICE)
+    "WEASD",  # water-equivalent snow depth (-> SNOW)
+    "TSOIL",  # soil temperature (-> ST<layer>)
+    "SOILW",  # volumetric soil moisture (-> SM<layer>)
+]
+
+# Full GFS isobaric set (hPa) for the 3D fields. Includes the upper levels
+# (down to 1 hPa, ~48 km) so a 30 km model top has data to interpolate from.
+# Levels absent from a given product are simply not returned; the converter
+# reports the actual count to set config_nfglevels.
+PRESSURE_LEVELS = [
+    1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500,
+    450, 400, 350, 300, 250, 200, 150, 100, 70, 50, 40, 30, 20, 15, 10,
+    7, 5, 3, 2, 1,
+]
+
+# Non-isobaric levels (grib_filter checkbox names).
+SURFACE_LEVELS = [
+    "surface",
+    "mean_sea_level",
+    "2_m_above_ground",
+    "10_m_above_ground",
+    "0-0.1_m_below_ground",
+    "0.1-0.4_m_below_ground",
+    "0.4-1_m_below_ground",
+    "1-2_m_below_ground",
+]
+
+NOMADS = "https://nomads.ncep.noaa.gov/cgi-bin"
+_UA = {"User-Agent": "MPAS-Research/usp-utils (research)"}
+
+
+# Minimal field set for an SST/sea-ice update file (--fields sst): just the
+# surface SST proxy + sea-ice + land mask.
+SST_VARIABLES = ["TMP", "ICEC", "LAND"]
+SST_SURFACE_LEVELS = ["surface"]
+
+
+def build_filter_url(date: str, cycle: str, fhour: int, res: str,
+                     fields: str = "full") -> tuple[str, str]:
+    """Return (url, grib_filename) for the NOMADS grib_filter request."""
+    ymd = date.replace("-", "")
+    res_tag = {"0p25": "0p25", "0p50": "0p50", "1p00": "1p00"}[res]
+    grib_name = f"gfs.t{cycle}z.pgrb2.{res_tag}.f{fhour:03d}"
+    sst_only = fields == "sst"
+    varis = SST_VARIABLES if sst_only else VARIABLES
+    plevs = [] if sst_only else PRESSURE_LEVELS
+    slevs = SST_SURFACE_LEVELS if sst_only else SURFACE_LEVELS
+    params = [
+        ("dir", f"/gfs.{ymd}/{cycle}/atmos"),
+        ("file", grib_name),
+    ]
+    for v in varis:
+        params.append((f"var_{v}", "on"))
+    for p in plevs:
+        params.append((f"lev_{p}_mb", "on"))
+    for lv in slevs:
+        params.append((f"lev_{lv}", "on"))
+    url = f"{NOMADS}/filter_gfs_{res_tag}.pl?" + urllib.parse.urlencode(params)
+    return url, grib_name
+
+
+# --- AWS S3 source (no rate limit; byte-range subset via the .idx file) -------
+AWS_BUCKET = "https://noaa-gfs-bdp-pds.s3.amazonaws.com"
+_IDX_3D = {"HGT", "TMP", "UGRD", "VGRD", "RH"}                  # all 'X mb' levels
+_IDX_2D = {("PRMSL", "mean sea level"), ("PRES", "surface"), ("HGT", "surface"),
+           ("TMP", "surface"), ("LAND", "surface"), ("ICEC", "surface"),
+           ("WEASD", "surface")}
+_IDX_SOIL = {"TSOIL", "SOILW"}                                  # any 'below ground'
+# Surface messages for an SST/sea-ice update file (--fields sst)
+_IDX_SST = {("TMP", "surface"), ("ICEC", "surface"), ("LAND", "surface")}
+
+
+def _aws_urls(date, cycle, fhour, res):
+    ymd = date.replace("-", "")
+    base = f"{AWS_BUCKET}/gfs.{ymd}/{cycle}/atmos"
+    name = f"gfs.t{cycle}z.pgrb2.{res}.f{fhour:03d}"
+    return f"{base}/{name}", f"{base}/{name}.idx", name
+
+
+def _get(url, rng=None, timeout=300):
+    headers = dict(_UA)
+    if rng:
+        headers["Range"] = rng
+    with urllib.request.urlopen(urllib.request.Request(url, headers=headers),
+                                timeout=timeout) as r:
+        return r.read()
+
+
+def fetch_aws(date, cycle, fhour, res, fields="full"):
+    """Download only the MPAS-needed GFS messages from AWS S3 via the .idx
+    (HTTP byte-range). Returns (grib_bytes, name). No rate limit.
+    fields='sst' fetches only the surface SST/sea-ice/land messages."""
+    grib_url, idx_url, name = _aws_urls(date, cycle, fhour, res)
+    lines = [ln for ln in _get(idx_url, timeout=120).decode("utf-8", "replace").splitlines()
+             if ln.strip()]
+    starts = [int(ln.split(":")[1]) for ln in lines]
+    sst_only = fields == "sst"
+    sel = []
+    for i, ln in enumerate(lines):
+        p = ln.split(":")
+        var, lev = p[3], p[4]
+        if sst_only:
+            keep = (var, lev) in _IDX_SST
+        else:
+            keep = ((var in _IDX_3D and lev.endswith("mb")) or ((var, lev) in _IDX_2D)
+                    or (var in _IDX_SOIL and "below ground" in lev))
+        if keep:
+            sel.append((starts[i], starts[i + 1] - 1 if i + 1 < len(starts) else None))
+    if not sel:
+        raise RuntimeError("no matching messages in .idx")
+    sel.sort()
+    merged = [list(sel[0])]
+    for s, e in sel[1:]:
+        if merged[-1][1] is not None and s == merged[-1][1] + 1:
+            merged[-1][1] = e
+        else:
+            merged.append([s, e])
+    chunks = [_get(grib_url, rng=f"bytes={s}-{'' if e is None else e}") for s, e in merged]
+    return b"".join(chunks), name
+
+
+# --- NCAR RDA ds084.1 source (historical GFS 0.25 deg, 2015-01-15+) -----------
+# Full GRIB2 files (no byte-range subset); the converter reads only the fields it
+# needs. Direct download from data.rda.ucar.edu; some access may require a free
+# RDA account (https://rda.ucar.edu) — set RDA_API_TOKEN to send a bearer token.
+RDA_BASE = "https://data.rda.ucar.edu/ds084.1"
+
+
+def _rda_url(date, cycle, fhour):
+    ymd = date.replace("-", "")
+    name = f"gfs.0p25.{ymd}{cycle}.f{fhour:03d}.grib2"
+    return f"{RDA_BASE}/{ymd[:4]}/{ymd}/{name}", name
+
+
+def fetch_rda_to(dest, date, cycle, fhour):
+    """Stream a full GFS 0.25 deg file from NCAR RDA ds084.1 to dest.
+    Returns (url, name). Tries anonymous access; set RDA_API_TOKEN if RDA
+    requires authentication."""
+    url, name = _rda_url(date, cycle, fhour)
+    headers = dict(_UA)
+    tok = os.environ.get("RDA_API_TOKEN")
+    if tok:
+        headers["Authorization"] = f"Bearer {tok}"
+    req = urllib.request.Request(url, headers=headers)
+    with urllib.request.urlopen(req, timeout=900) as r, open(dest, "wb") as f:
+        shutil.copyfileobj(r, f)
+    return url, name
+
+
+def check_coverage(source: str, d: date) -> str | None:
+    """Return an error message if date `d` is outside `source`'s coverage, else None."""
+    today = datetime.now(timezone.utc).date()
+    if source == "aws" and d < AWS_START:
+        return (f"AWS GFS archive (gfs.*/atmos/ layout) starts {AWS_START.isoformat()}. "
+                f"For {RDA_START.isoformat()}..{AWS_START.isoformat()} use --source rda; "
+                f"for older dates use a reanalysis (download_era5.py).")
+    if source == "rda" and d < RDA_START:
+        return (f"GFS 0.25 deg only exists from {RDA_START.isoformat()} (NCAR RDA ds084.1); "
+                f"nothing earlier exists at 0.25 deg. For older dates use a reanalysis "
+                f"(download_era5.py — ERA5 goes back to 1940).")
+    if source == "nomads" and (today - d).days > NOMADS_DAYS:
+        return (f"NOMADS keeps only ~{NOMADS_DAYS} days (requested {(today - d).days} days ago). "
+                f"Use --source aws ({AWS_START.isoformat()}+) or rda ({RDA_START.isoformat()}+), "
+                f"or a reanalysis (download_era5.py) for older dates.")
+    return None
+
+
+SRC_LABEL = {
+    "aws": "AWS S3 (noaa-gfs-bdp-pds, idx byte-range)",
+    "nomads": "NOMADS grib_filter",
+    "rda": "NCAR RDA ds084.1 (full file)",
+}
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="Download GFS fields for MPAS init.")
+    ap.add_argument("--date", required=True, help="Cycle date, YYYY-MM-DD (UTC)")
+    ap.add_argument("--cycle", default="00", choices=["00", "06", "12", "18"])
+    ap.add_argument("--fhour", type=int, default=0, help="Forecast hour (0 = analysis)")
+    ap.add_argument("--res", default="0p25", choices=["0p25", "0p50", "1p00"])
+    ap.add_argument("--source", default="aws", choices=["aws", "nomads", "rda"],
+                    help="aws = S3 byte-range subset (2021-03-23+, no rate limit, default); "
+                         "nomads = grib_filter (~last 10 days); "
+                         "rda = NCAR RDA ds084.1 historical (2015-01-15+, full files)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: $MPAS_ROOT/met_data/gfs/<date><cycle>)")
+    ap.add_argument("--fields", default="full", choices=["full", "sst"],
+                    help="full = all fields MPAS needs (default); "
+                         "sst = only surface SST/sea-ice/land (small download, for "
+                         "SST update files via gfs_sst_to_intermediate.py)")
+    args = ap.parse_args()
+
+    # Validate the date and that the chosen source actually covers it.
+    try:
+        d = datetime.strptime(args.date, "%Y-%m-%d").date()
+    except ValueError:
+        print("[ERROR] --date must be YYYY-MM-DD", file=sys.stderr)
+        return 2
+    cov_err = check_coverage(args.source, d)
+    if cov_err:
+        print(f"[ERROR] {cov_err}", file=sys.stderr)
+        return 1
+    if args.source == "rda" and args.res != "0p25":
+        print("[ERROR] NCAR RDA ds084.1 is GFS 0.25 deg only; use --res 0p25.",
+              file=sys.stderr)
+        return 1
+
+    ymd = args.date.replace("-", "")
+    repo_root = Path(__file__).resolve().parents[3]
+    outdir = Path(args.outdir) if args.outdir else repo_root / "met_data" / "gfs" / f"{ymd}{args.cycle}"
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    if args.source == "aws":
+        _, _, grib_name = _aws_urls(args.date, args.cycle, args.fhour, args.res)
+        url = "AWS S3 byte-range (.idx)"
+    elif args.source == "rda":
+        url, grib_name = _rda_url(args.date, args.cycle, args.fhour)
+        if args.fields == "sst":
+            print("[INFO] RDA serves full files (no subset); downloading the whole "
+                  "file and using only its SST/sea-ice/land fields.")
+    else:
+        url, grib_name = build_filter_url(args.date, args.cycle, args.fhour, args.res,
+                                          fields=args.fields)
+    # SST-only subsets get a distinct name so they never collide with a full
+    # download in the same directory (subsettable sources only).
+    if args.fields == "sst" and args.source in ("aws", "nomads"):
+        grib_name += ".sst"
+    dest = outdir / grib_name
+
+    if dest.exists() and dest.stat().st_size > 0:
+        print(f"[OK] {grib_name} already present — skipping ({dest})")
+        return 0
+
+    print(f"[INFO] Downloading {grib_name} from {args.source} -> {dest}")
+
+    if args.source == "rda":
+        # Full file, streamed to disk (no subset).
+        try:
+            url, grib_name = fetch_rda_to(dest, args.date, args.cycle, args.fhour)
+        except urllib.error.HTTPError as e:
+            dest.unlink(missing_ok=True)
+            if e.code in (401, 403):
+                print("[ERROR] RDA requires authentication. Create a free account at "
+                      "https://rda.ucar.edu and set RDA_API_TOKEN.", file=sys.stderr)
+            else:
+                print(f"[ERROR] RDA download failed: HTTP {e.code} "
+                      f"(date/cycle not in ds084.1?)", file=sys.stderr)
+            return 1
+        except Exception as e:  # noqa: BLE001
+            dest.unlink(missing_ok=True)
+            print(f"[ERROR] RDA download failed: {e}", file=sys.stderr)
+            return 1
+        with open(dest, "rb") as f:
+            head = f.read(4)
+        if head != b"GRIB":
+            dest.unlink(missing_ok=True)
+            print("[ERROR] RDA response is not GRIB (login/HTML page?). "
+                  "Set RDA_API_TOKEN if the dataset needs authentication.", file=sys.stderr)
+            return 1
+        size_bytes = dest.stat().st_size
+        print(f"[OK] {grib_name} ({size_bytes / 1e6:.1f} MB)")
+    else:
+        if args.source == "aws":
+            try:
+                data, _ = fetch_aws(args.date, args.cycle, args.fhour, args.res,
+                                    fields=args.fields)
+            except urllib.error.HTTPError as e:
+                print(f"[ERROR] AWS download failed: HTTP {e.code} "
+                      f"(date/cycle not available?)", file=sys.stderr)
+                return 1
+            except Exception as e:  # noqa: BLE001
+                print(f"[ERROR] AWS download failed: {e}", file=sys.stderr)
+                return 1
+        else:
+            # NOMADS grib_filter (descriptive User-Agent is good etiquette).
+            req = urllib.request.Request(url, headers={
+                "User-Agent": "MPAS-Research/usp-utils (research; NOMADS grib_filter)"})
+            try:
+                with urllib.request.urlopen(req, timeout=600) as r:
+                    data = r.read()
+            except urllib.error.HTTPError as e:
+                data = e.read() or b""
+            except Exception as e:  # noqa: BLE001
+                print(f"[ERROR] Download failed: {e}", file=sys.stderr)
+                return 1
+
+        # NOMADS returns an HTML body (often with a 302) instead of GRIB on errors.
+        if data[:4] != b"GRIB":
+            if b"Over Rate Limit" in data or b"abusive-user-block" in data:
+                print("[ERROR] NOMADS rate limit reached ('Over Rate Limit'). This is a "
+                      "temporary IP block from too many requests.", file=sys.stderr)
+                print("        Wait ~1 hour and retry, and avoid rapid repeated downloads. "
+                      "See https://www.weather.gov/abusive-user-block", file=sys.stderr)
+            else:
+                print("[ERROR] Response is not GRIB (bad date/cycle, or data not yet "
+                      "available / purged from NOMADS).", file=sys.stderr)
+                print("        First bytes:", data[:120], file=sys.stderr)
+            return 1
+
+        dest.write_bytes(data)
+        size_bytes = len(data)
+        print(f"[OK] {grib_name} ({size_bytes / 1e6:.1f} MB)")
+
+    # --- Provenance ----------------------------------------------------------
+    sst_only = args.fields == "sst"
+    meta = {
+        "dataset": "GFS",
+        "product": f"pgrb2.{args.res}",
+        "source": SRC_LABEL[args.source],
+        "fields": args.fields,
+        "cycle": f"{args.date} {args.cycle}:00 UTC",
+        "forecast_hour": args.fhour,
+        "variables": SST_VARIABLES if sst_only else VARIABLES,
+        "pressure_levels_hpa": [] if sst_only else PRESSURE_LEVELS,
+        "surface_levels": SST_SURFACE_LEVELS if sst_only else SURFACE_LEVELS,
+        "subset": args.source != "rda",  # rda = full file
+        "url": url,
+        "file": str(dest),
+        "size_bytes": size_bytes,
+        "download_date": datetime.now(timezone.utc).isoformat(),
+    }
+    (outdir / f"{grib_name}.provenance.json").write_text(json.dumps(meta, indent=2))
+    print(f"[INFO] Provenance written to {grib_name}.provenance.json")
+    next_script = "gfs_sst_to_intermediate.py" if sst_only else "gfs_to_intermediate.py"
+    print(f"[INFO] Next: python {next_script} --grib {dest}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/download_oisst.py b/usp-utils/pre_proc/real_data/download_oisst.py
new file mode 100755
index 00000000..14f3567a
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/download_oisst.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+"""
+download_oisst.py — fetch NOAA OISST v2.1 (AVHRR) daily 0.25 deg SST + sea-ice
+NetCDF from NCEI, for building MPAS SST/sea-ice update files.
+
+OISST is the daily 0.25 deg optimum-interpolation SST analysis (Sep 1981 ->
+present); it is the recommended lower-boundary SST for the scientific /
+downscaling pipeline (ERA5 atmosphere + OISST SST). The downloaded files are
+turned into WPS intermediate `SST:` files by oisst_to_intermediate.py, then fed
+to init_atmosphere config_init_case=8 (see the README "SST / sea-ice update").
+
+Each daily file (oisst-avhrr-v02r01.YYYYMMDD.nc) holds: sst (degC), ice
+(fraction), anom, err on a 1440x720 global 0.25 deg grid. Files in the last
+~2 weeks are served with a `_preliminary` suffix; this script falls back to it.
+
+Usage:
+    python download_oisst.py --start 2014-09-01 [--end 2014-09-15] [--outdir <dir>]
+"""
+
+import argparse
+import json
+import sys
+import urllib.error
+import urllib.request
+from datetime import date, datetime, timedelta, timezone
+from pathlib import Path
+
+OISST_BASE = ("https://www.ncei.noaa.gov/data/"
+              "sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr")
+OISST_START = date(1981, 9, 1)
+_UA = {"User-Agent": "MPAS-Research/usp-utils (research)"}
+
+
+def _urls(d: date):
+    ym = d.strftime("%Y%m")
+    ymd = d.strftime("%Y%m%d")
+    final = f"{OISST_BASE}/{ym}/oisst-avhrr-v02r01.{ymd}.nc"
+    prelim = f"{OISST_BASE}/{ym}/oisst-avhrr-v02r01.{ymd}_preliminary.nc"
+    return final, prelim
+
+
+def _fetch_to(url, dest) -> bool:
+    """Stream url to dest. Return True on success, False on 404."""
+    import shutil
+    req = urllib.request.Request(url, headers=_UA)
+    try:
+        with urllib.request.urlopen(req, timeout=300) as r, open(dest, "wb") as f:
+            shutil.copyfileobj(r, f)
+        return True
+    except urllib.error.HTTPError as e:
+        if e.code == 404:
+            return False
+        raise
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="Download OISST v2.1 daily SST for MPAS SST update.")
+    ap.add_argument("--start", required=True, help="First date, YYYY-MM-DD (UTC)")
+    ap.add_argument("--end", default=None, help="Last date, YYYY-MM-DD (default: = start)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: $MPAS_ROOT/met_data/oisst)")
+    args = ap.parse_args()
+
+    try:
+        d0 = datetime.strptime(args.start, "%Y-%m-%d").date()
+        d1 = datetime.strptime(args.end, "%Y-%m-%d").date() if args.end else d0
+    except ValueError:
+        print("[ERROR] dates must be YYYY-MM-DD", file=sys.stderr)
+        return 2
+    if d1 < d0:
+        print("[ERROR] --end is before --start", file=sys.stderr)
+        return 2
+    if d0 < OISST_START:
+        print(f"[ERROR] OISST v2.1 starts {OISST_START.isoformat()}.", file=sys.stderr)
+        return 1
+
+    repo_root = Path(__file__).resolve().parents[3]
+    outdir = Path(args.outdir) if args.outdir else repo_root / "met_data" / "oisst"
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    files, d = [], d0
+    while d <= d1:
+        final, prelim = _urls(d)
+        name = Path(final).name
+        dest = outdir / name
+        if dest.exists() and dest.stat().st_size > 0:
+            print(f"[OK] {name} already present — skipping")
+            files.append(str(dest))
+            d += timedelta(days=1)
+            continue
+        print(f"[INFO] {d.isoformat()} -> {dest}")
+        try:
+            if _fetch_to(final, dest):
+                files.append(str(dest))
+            else:
+                # fall back to the preliminary file (recent dates)
+                pdest = outdir / Path(prelim).name
+                if _fetch_to(prelim, pdest):
+                    print(f"[INFO]   using preliminary file for {d.isoformat()}")
+                    files.append(str(pdest))
+                else:
+                    print(f"[WARN]   no OISST file for {d.isoformat()} (final or preliminary)",
+                          file=sys.stderr)
+        except Exception as e:  # noqa: BLE001
+            print(f"[ERROR] download failed for {d.isoformat()}: {e}", file=sys.stderr)
+            return 1
+        d += timedelta(days=1)
+
+    if not files:
+        print("[ERROR] no files downloaded", file=sys.stderr)
+        return 1
+
+    meta = {
+        "dataset": "OISST",
+        "product": "NOAA OI SST v2.1 (AVHRR), daily 0.25 deg",
+        "source": "NCEI",
+        "base_url": OISST_BASE,
+        "period": {"start": d0.isoformat(), "end": d1.isoformat()},
+        "variables_used": ["sst (degC)", "ice (fraction)"],
+        "files": files,
+        "download_date": datetime.now(timezone.utc).isoformat(),
+    }
+    (outdir / f"oisst_{d0.strftime('%Y%m%d')}_{d1.strftime('%Y%m%d')}.provenance.json"
+     ).write_text(json.dumps(meta, indent=2))
+    print(f"[OK] {len(files)} OISST file(s) in {outdir}")
+    print(f"[INFO] Next: python oisst_to_intermediate.py --indir {outdir} "
+          f"--start {d0.isoformat()} --end {d1.isoformat()}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/era5_to_intermediate.py b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
new file mode 100755
index 00000000..76d73273
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
@@ -0,0 +1,176 @@
+#!/usr/bin/env python3
+"""
+era5_to_intermediate.py — convert ERA5 GRIB files (from download_era5.py) into a
+WPS intermediate-format file (e.g. GFS:2014-09-10_00) that MPAS init_atmosphere
+reads for real-data initialization (config_met_prefix='GFS').
+
+This is the ERA5 counterpart of gfs_to_intermediate.py (same pywinter output,
+same intermediate field names), for dates the GFS archive does not cover.
+
+Run inside the conda env that has cfgrib + pywinter:
+    conda run -n cgfd-usp-mpas python era5_to_intermediate.py \
+        --pl era5_pl_YYYYMMDDHH.grib --sl era5_sl_YYYYMMDDHH.grib
+
+ERA5 specifics handled here:
+- Geopotential (z, m^2/s^2) is divided by g0=9.80665 to get height (GHT/SOILHGT, m).
+- Snow depth (sd, m of water equivalent) is scaled to kg m-2 (SNOW) by *1000.
+- Pressure levels are written in Pa (WPS requirement; see the units note in README).
+- Latitude is flipped to ascending (S->N), as WPS expects (SW-corner start).
+- Soil layers use the ERA5 depths (0-7, 7-28, 28-100, 100-289 cm).
+
+NOTE: this converter mirrors the validated GFS path but the ERA5 branch has not
+yet been validated end-to-end here — verify the resulting init.nc with the
+sanity check in the README before running the model.
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+import pywinter.winter as pw
+
+G0 = 9.80665  # standard gravity (m s-2), for geopotential -> height
+
+# 3D isobaric: cfgrib shortName -> (intermediate name, transform or None)
+ISOBARIC = {
+    "z": ("GHT", lambda a: a / G0),
+    "t": ("TT", None),
+    "u": ("UU", None),
+    "v": ("VV", None),
+    "r": ("RH", None),
+}
+# 2D surface (typeOfLevel='surface'): shortName -> (intermediate name, transform)
+SURFACE = {
+    "sp": ("PSFC", None),
+    "skt": ("SKINTEMP", None),
+    "z": ("SOILHGT", lambda a: a / G0),
+    "lsm": ("LANDSEA", None),
+    "siconc": ("SEAICE", None),
+    "sd": ("SNOW", lambda a: a * 1000.0),  # m w.e. -> kg m-2
+}
+# Soil layers: ERA5 shortNames and the matching WPS depth tags (top-bottom, cm)
+SOIL_T = ["stl1", "stl2", "stl3", "stl4"]
+SOIL_M = ["swvl1", "swvl2", "swvl3", "swvl4"]
+SOIL_TAGS = ["000007", "007028", "028100", "100289"]
+
+
+def _open(path, **keys):
+    return xr.open_dataset(
+        path, engine="cfgrib",
+        backend_kwargs={"filter_by_keys": keys, "indexpath": ""},
+    )
+
+
+def _open_var(path, sn, **keys):
+    """Open a single variable by shortName; return the DataArray or None."""
+    try:
+        ds = _open(path, shortName=sn, **keys)
+    except Exception as e:  # noqa: BLE001
+        print(f"[WARN] could not read {sn}: {e}")
+        return None
+    return ds[sn] if sn in ds else None
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="ERA5 GRIB -> WPS intermediate (pywinter).")
+    ap.add_argument("--pl", required=True, help="ERA5 pressure-level GRIB file")
+    ap.add_argument("--sl", required=True, help="ERA5 single-level GRIB file")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: the --pl file's directory)")
+    ap.add_argument("--prefix", default="GFS",
+                    help="Intermediate prefix (config_met_prefix; keep 'GFS')")
+    args = ap.parse_args()
+
+    pl, sl = Path(args.pl), Path(args.sl)
+    for f in (pl, sl):
+        if not f.exists():
+            print(f"[ERROR] not found: {f}", file=sys.stderr)
+            return 1
+    outdir = Path(args.outdir) if args.outdir else pl.parent
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    print(f"[INFO] Reading {pl} (pressure levels)")
+    iso = _open(pl, typeOfLevel="isobaricInhPa")
+
+    # Grid geometry (flip latitude to ascending = S->N, SW-corner start)
+    lat = iso.latitude.values
+    lon = iso.longitude.values
+    flip = lat[0] > lat[-1]
+    dlat = abs(float(lat[1] - lat[0]))
+    dlon = abs(float(lon[1] - lon[0]))
+    geo = pw.Geo0(float(min(lat)), float(min(lon)), dlat, dlon)
+
+    def arr2d(a):
+        a = np.asarray(a, dtype="float32")
+        return a[::-1, :] if flip else a
+
+    def arr3d(a):
+        a = np.asarray(a, dtype="float32")
+        return a[:, ::-1, :] if flip else a
+
+    variables = []
+
+    # --- 3D isobaric ---
+    plevs_hpa = np.asarray(iso.isobaricInhPa.values, dtype="float32")  # hPa
+    plevs = plevs_hpa * 100.0  # WPS intermediate / pywinter V3dp expect Pa
+    print(f"[INFO] {plevs.size} isobaric levels ({int(plevs_hpa.max())}..{int(plevs_hpa.min())} hPa)"
+          f" -> set config_nfglevels = {plevs.size + 1}")
+    for sn, (name, fn) in ISOBARIC.items():
+        if sn in iso:
+            a = arr3d(iso[sn].values)
+            variables.append(pw.V3dp(name, fn(a) if fn else a, plevs))
+        else:
+            print(f"[WARN] isobaric field missing: {sn}")
+
+    # --- 2D surface (open each var individually to avoid cfgrib hypercube clashes) ---
+    print(f"[INFO] Reading {sl} (single levels)")
+    for sn, (name, fn) in SURFACE.items():
+        da = _open_var(sl, sn, typeOfLevel="surface")
+        if da is not None:
+            a = arr2d(da.values)
+            variables.append(pw.V2d(name, fn(a) if fn else a))
+        else:
+            print(f"[WARN] surface field missing: {sn}")
+    msl = _open_var(sl, "msl", typeOfLevel="meanSea")
+    if msl is not None:
+        variables.append(pw.V2d("PMSL", arr2d(msl.values)))
+    else:
+        print("[WARN] surface field missing: msl")
+
+    # --- Soil layers (stack the 4 ERA5 layers; pywinter Vsl -> ST/SM<tag>) ---
+    st = [_open_var(sl, s, typeOfLevel="depthBelowLandLayer") for s in SOIL_T]
+    sm = [_open_var(sl, s, typeOfLevel="depthBelowLandLayer") for s in SOIL_M]
+    if all(x is not None for x in st):
+        variables.append(pw.Vsl("ST", np.stack([arr2d(x.values) for x in st]), SOIL_TAGS))
+    else:
+        print("[WARN] soil temperature layers incomplete — skipping ST")
+    if all(x is not None for x in sm):
+        variables.append(pw.Vsl("SM", np.stack([arr2d(x.values) for x in sm]), SOIL_TAGS))
+    else:
+        print("[WARN] soil moisture layers incomplete — skipping SM")
+
+    # Date tag (YYYY-MM-DD_HH) from the analysis valid time
+    t = iso.valid_time.values if "valid_time" in iso else iso.time.values
+    vt = pd.to_datetime(np.atleast_1d(t)[0])
+    date_tag = vt.strftime("%Y-%m-%d_%H")
+
+    print(f"[INFO] Writing {len(variables)} fields -> {outdir}/{args.prefix}:{date_tag}")
+    pw.cinter(args.prefix, date_tag, geo, variables, str(outdir) + "/")
+
+    out = outdir / f"{args.prefix}:{date_tag}"
+    if out.exists():
+        print(f"[OK] {out} ({out.stat().st_size / 1e6:.1f} MB)")
+        print(f"[INFO] Link it in the run dir and set config_met_prefix='{args.prefix}', "
+              f"config_start_time='{date_tag}:00:00'")
+        print("[INFO] Verify the init.nc with the sanity check in the README before "
+              "running the model.")
+        return 0
+    print(f"[ERROR] expected output not found: {out}", file=sys.stderr)
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/gfs_sst_to_intermediate.py b/usp-utils/pre_proc/real_data/gfs_sst_to_intermediate.py
new file mode 100755
index 00000000..1e9e5e8a
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/gfs_sst_to_intermediate.py
@@ -0,0 +1,98 @@
+#!/usr/bin/env python3
+"""
+gfs_sst_to_intermediate.py — extract SST / sea-ice / land-sea mask from a GFS
+GRIB2 file (the same file download_gfs.py already fetches for the atmosphere)
+and write a WPS intermediate `SST:YYYY-MM-DD_HH` file for MPAS SST/sea-ice
+updates. This is the operational counterpart of oisst_to_intermediate.py.
+
+The three fields MPAS expects for an SST update come from the GFS surface
+records already in the GRIB:
+    TMP:surface  -> SST     (skin temperature; over water this is SST, and
+                             init_atmosphere masks land via LANDSEA)
+    ICEC:surface -> SEAICE  (sea-ice fraction)
+    LAND:surface -> LANDSEA (land-sea mask, 1 = land, 0 = water)
+
+Feed the resulting `SST:` files to init_atmosphere with config_init_case=8
+(see the README "SST / sea-ice update").
+
+Run inside the conda env that has cfgrib + pywinter:
+    conda run -n cgfd-usp-mpas python gfs_sst_to_intermediate.py --grib <file.f000>
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+import pywinter.winter as pw
+
+
+def _open(grib, **keys):
+    return xr.open_dataset(
+        grib, engine="cfgrib",
+        backend_kwargs={"filter_by_keys": keys, "indexpath": ""},
+    )
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="GFS GRIB -> WPS intermediate SST: file.")
+    ap.add_argument("--grib", required=True, help="GFS GRIB2 file (e.g. *.f000)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: the GRIB's directory)")
+    ap.add_argument("--prefix", default="SST", help="Intermediate prefix (config_sfc_prefix)")
+    args = ap.parse_args()
+
+    grib = Path(args.grib)
+    if not grib.exists():
+        print(f"[ERROR] not found: {grib}", file=sys.stderr)
+        return 1
+    outdir = Path(args.outdir) if args.outdir else grib.parent
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    print(f"[INFO] Reading surface fields from {grib}")
+    sfc = _open(grib, typeOfLevel="surface")
+
+    lat = sfc.latitude.values
+    lon = sfc.longitude.values
+    flip = lat[0] > lat[-1]
+    dlat = abs(float(lat[1] - lat[0]))
+    dlon = abs(float(lon[1] - lon[0]))
+    geo = pw.Geo0(float(min(lat)), float(min(lon)), dlat, dlon)
+
+    def arr2d(da):
+        a = np.asarray(da.values, dtype="float32")
+        return a[::-1, :] if flip else a
+
+    # GFS surface shortNames: t (skin temp), siconc (sea ice), lsm (land mask)
+    mapping = {"t": "SST", "siconc": "SEAICE", "lsm": "LANDSEA"}
+    variables, added = [], set()
+    for sn, name in mapping.items():
+        if sn in sfc:
+            variables.append(pw.V2d(name, arr2d(sfc[sn])))
+            added.add(name)
+        else:
+            print(f"[WARN] surface field missing: {sn} (-> {name})", file=sys.stderr)
+    if "SST" not in added:
+        print("[ERROR] no surface temperature (SST) in the GRIB — cannot build SST file.",
+              file=sys.stderr)
+        return 1
+
+    vt = pd.to_datetime(np.atleast_1d(sfc.valid_time.values)[0])
+    date_tag = vt.strftime("%Y-%m-%d_%H")
+
+    print(f"[INFO] Writing {len(variables)} fields -> {outdir}/{args.prefix}:{date_tag}")
+    pw.cinter(args.prefix, date_tag, geo, variables, str(outdir) + "/")
+    out = outdir / f"{args.prefix}:{date_tag}"
+    if out.exists():
+        print(f"[OK] {out.name} ({out.stat().st_size / 1e6:.1f} MB)")
+        print(f"[INFO] Next: run init_atmosphere config_init_case=8 with "
+              f"config_sfc_prefix='{args.prefix}' (see README 'SST / sea-ice update').")
+        return 0
+    print(f"[ERROR] expected output not found: {out}", file=sys.stderr)
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/gfs_to_intermediate.py b/usp-utils/pre_proc/real_data/gfs_to_intermediate.py
new file mode 100755
index 00000000..e81163a2
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/gfs_to_intermediate.py
@@ -0,0 +1,128 @@
+#!/usr/bin/env python3
+"""
+gfs_to_intermediate.py — convert a GFS GRIB2 file (from download_gfs.py) into a
+WPS intermediate-format file (e.g. GFS:2024-09-10_00) that MPAS init_atmosphere
+reads for real-data initialization (config_met_prefix='GFS').
+
+Uses cfgrib (to read the GRIB) and pywinter (to write the WPS intermediate
+format) — no need to build WPS/ungrib.
+
+Run inside the conda env that has cfgrib + pywinter:
+    conda run -n cgfd-usp-mpas python gfs_to_intermediate.py --grib <file.f000>
+                                       [--outdir <dir>]
+
+The GFS grid is global regular lat/lon, stored N->S; the latitude axis is
+flipped to S->N so the intermediate file starts at the SW corner, as WPS expects.
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+import pywinter.winter as pw
+
+# 3D isobaric fields: GRIB shortName -> intermediate name
+ISOBARIC = {"gh": "GHT", "t": "TT", "u": "UU", "v": "VV", "r": "RH"}
+# 2D surface fields (typeOfLevel='surface'): shortName -> intermediate name
+SURFACE = {"sp": "PSFC", "orog": "SOILHGT", "t": "SKINTEMP",
+           "sdwe": "SNOW", "lsm": "LANDSEA", "siconc": "SEAICE"}
+# Soil layers (depthBelowLandLayer tops, m) -> WPS depth tag (top-bottom, cm)
+SOIL_TAGS = ["000010", "010040", "040100", "100200"]
+
+
+def _open(grib, **keys):
+    return xr.open_dataset(
+        grib, engine="cfgrib",
+        backend_kwargs={"filter_by_keys": keys, "indexpath": ""},
+    )
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="GFS GRIB -> WPS intermediate (pywinter).")
+    ap.add_argument("--grib", required=True, help="GFS GRIB2 file (e.g. *.f000)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir for the intermediate file "
+                         "(default: the GRIB's directory)")
+    ap.add_argument("--prefix", default="GFS", help="Intermediate prefix (config_met_prefix)")
+    args = ap.parse_args()
+
+    grib = Path(args.grib)
+    if not grib.exists():
+        print(f"[ERROR] not found: {grib}", file=sys.stderr)
+        return 1
+    outdir = Path(args.outdir) if args.outdir else grib.parent
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    print(f"[INFO] Reading {grib}")
+    iso = _open(grib, typeOfLevel="isobaricInhPa")
+    sfc = _open(grib, typeOfLevel="surface")
+    msl = _open(grib, typeOfLevel="meanSea")
+    soil = _open(grib, typeOfLevel="depthBelowLandLayer")
+
+    # Grid geometry (flip latitude to ascending = S->N, SW-corner start)
+    lat = iso.latitude.values
+    lon = iso.longitude.values
+    flip = lat[0] > lat[-1]
+    dlat = abs(float(lat[1] - lat[0]))
+    dlon = abs(float(lon[1] - lon[0]))
+    geo = pw.Geo0(float(min(lat)), float(min(lon)), dlat, dlon)
+
+    def arr2d(da):
+        a = np.asarray(da.values, dtype="float32")
+        return a[::-1, :] if flip else a
+
+    def arr3d(da):
+        a = np.asarray(da.values, dtype="float32")
+        return a[:, ::-1, :] if flip else a
+
+    variables = []
+
+    # --- 3D isobaric ---
+    plevs_hpa = np.asarray(iso.isobaricInhPa.values, dtype="float32")  # hPa
+    plevs = plevs_hpa * 100.0  # WPS intermediate / pywinter V3dp expect pressure in Pa
+    print(f"[INFO] {plevs.size} isobaric levels ({int(plevs_hpa.max())}..{int(plevs_hpa.min())} hPa)"
+          f" -> set config_nfglevels = {plevs.size + 1}")
+    for sn, name in ISOBARIC.items():
+        if sn in iso:
+            variables.append(pw.V3dp(name, arr3d(iso[sn]), plevs))
+        else:
+            print(f"[WARN] isobaric field missing: {sn}")
+
+    # --- 2D surface + mean sea level ---
+    for sn, name in SURFACE.items():
+        if sn in sfc:
+            variables.append(pw.V2d(name, arr2d(sfc[sn])))
+        else:
+            print(f"[WARN] surface field missing: {sn}")
+    if "prmsl" in msl:
+        variables.append(pw.V2d("PMSL", arr2d(msl["prmsl"])))
+
+    # --- Soil layers (pywinter Vsl: name 'ST'/'SM' + numeric-string level tags
+    #     -> intermediate fields ST000010, SM000010, ...) ---
+    if "st" in soil:
+        variables.append(pw.Vsl("ST", arr3d(soil["st"]), SOIL_TAGS))
+    if "soilw" in soil:
+        variables.append(pw.Vsl("SM", arr3d(soil["soilw"]), SOIL_TAGS))
+
+    # Date tag (YYYY-MM-DD_HH) from the analysis valid time
+    vt = pd.to_datetime(np.atleast_1d(iso.valid_time.values)[0])
+    date_tag = vt.strftime("%Y-%m-%d_%H")
+
+    print(f"[INFO] Writing {len(variables)} fields -> {outdir}/{args.prefix}:{date_tag}")
+    pw.cinter(args.prefix, date_tag, geo, variables, str(outdir) + "/")
+
+    out = outdir / f"{args.prefix}:{date_tag}"
+    if out.exists():
+        print(f"[OK] {out} ({out.stat().st_size/1e6:.1f} MB)")
+        print(f"[INFO] Link it in the run dir and set config_met_prefix='{args.prefix}', "
+              f"config_start_time='{date_tag}:00:00'")
+        return 0
+    print(f"[ERROR] expected output not found: {out}", file=sys.stderr)
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/oisst_clim_to_intermediate.py b/usp-utils/pre_proc/real_data/oisst_clim_to_intermediate.py
new file mode 100755
index 00000000..8b768bd2
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/oisst_clim_to_intermediate.py
@@ -0,0 +1,126 @@
+#!/usr/bin/env python3
+"""
+oisst_clim_to_intermediate.py — build MPAS SST/sea-ice update files from the
+NOAA OISST daily *climatology* (long-term mean, 1991-2020), for use as a
+climatological lower boundary in forecasts / idealized runs where there is no
+observed SST for the target dates (e.g. a forecast into the future, or a run
+representing a "typical" month).
+
+For each target date it picks the matching day-of-year from the climatology and
+writes an `SST:YYYY-MM-DD_HH` intermediate (fields SST, SEAICE, LANDSEA) tagged
+with that target date — exactly like the observed pipeline, so the rest of the
+SST-update workflow (init_atmosphere config_init_case=8) is identical.
+
+The 0.25° daily climatology is read on demand over OPeNDAP (only the needed
+day-of-year slices are fetched, a few MB each), so there is no multi-GB download.
+Source: NOAA PSL (https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres.html).
+
+Run inside the conda env that has xarray (with OPeNDAP) + pywinter:
+    conda run -n cgfd-usp-mpas python oisst_clim_to_intermediate.py \
+        --start 2026-07-01 --end 2026-07-31 [--hour 00] [--outdir <dir>]
+"""
+
+import argparse
+import calendar
+import sys
+from datetime import datetime, timedelta
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
+import pywinter.winter as pw
+
+PSL = "https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres"
+SST_LTM = f"{PSL}/sst.day.mean.ltm.1991-2020.nc"
+ICE_LTM = f"{PSL}/icec.day.mean.ltm.1991-2020.nc"
+LAND_SST_FILL_K = 273.15
+
+
+def clim_index(d) -> int:
+    """0-based index into a 365-day (non-leap) day-of-year climatology.
+    Feb 29 maps to Feb 28; later leap-year days are shifted back by one."""
+    doy = d.timetuple().tm_yday
+    if calendar.isleap(d.year) and doy >= 60:
+        doy = 59 if doy == 60 else doy - 1
+    return doy - 1
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="OISST daily climatology -> WPS SST: files.")
+    ap.add_argument("--start", required=True, help="First target date, YYYY-MM-DD")
+    ap.add_argument("--end", default=None, help="Last target date (default: = start)")
+    ap.add_argument("--hour", default="00", help="Hour tag HH for the daily field (default 00)")
+    ap.add_argument("--outdir", default=None,
+                    help="Output dir (default: $MPAS_ROOT/met_data/oisst_clim)")
+    ap.add_argument("--prefix", default="SST", help="Intermediate prefix (config_sfc_prefix)")
+    ap.add_argument("--sst-url", default=SST_LTM, help="Override SST climatology URL/path")
+    ap.add_argument("--ice-url", default=ICE_LTM, help="Override sea-ice climatology URL/path")
+    args = ap.parse_args()
+
+    try:
+        d0 = datetime.strptime(args.start, "%Y-%m-%d").date()
+        d1 = datetime.strptime(args.end, "%Y-%m-%d").date() if args.end else d0
+    except ValueError:
+        print("[ERROR] dates must be YYYY-MM-DD", file=sys.stderr)
+        return 2
+    if d1 < d0:
+        print("[ERROR] --end is before --start", file=sys.stderr)
+        return 2
+
+    repo_root = Path(__file__).resolve().parents[3]
+    outdir = Path(args.outdir) if args.outdir else repo_root / "met_data" / "oisst_clim"
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    print(f"[INFO] Opening OISST climatology (OPeNDAP): {args.sst_url}")
+    try:
+        sst_ds = xr.open_dataset(args.sst_url, decode_times=False)
+        ice_ds = xr.open_dataset(args.ice_url, decode_times=False)
+    except Exception as e:  # noqa: BLE001
+        print(f"[ERROR] could not open climatology: {e}", file=sys.stderr)
+        return 1
+    sst_var = "sst" if "sst" in sst_ds else list(sst_ds.data_vars)[0]
+    ice_var = "icec" if "icec" in ice_ds else list(ice_ds.data_vars)[0]
+
+    lat = sst_ds["lat"].values
+    lon = sst_ds["lon"].values
+    flip = lat[0] > lat[-1]
+    dlat = abs(float(lat[1] - lat[0]))
+    dlon = abs(float(lon[1] - lon[0]))
+    geo = pw.Geo0(float(min(lat)), float(min(lon)), dlat, dlon)
+
+    def arr2d(a):
+        a = np.asarray(a, dtype="float32")
+        return a[::-1, :] if flip else a
+
+    n, d = 0, d0
+    while d <= d1:
+        idx = clim_index(d)
+        sst_c = np.asarray(sst_ds[sst_var].isel(time=idx).values, dtype="float64")
+        ice_c = np.asarray(ice_ds[ice_var].isel(time=idx).values, dtype="float64")
+        land = ~np.isfinite(sst_c)
+        sst_k = np.where(land, LAND_SST_FILL_K, sst_c + 273.15)
+        seaice = np.nan_to_num(ice_c, nan=0.0)
+        landsea = land.astype("float32")
+        variables = [
+            pw.V2d("SST", arr2d(sst_k)),
+            pw.V2d("SEAICE", arr2d(seaice)),
+            pw.V2d("LANDSEA", arr2d(landsea)),
+        ]
+        date_tag = f"{d.isoformat()}_{args.hour.zfill(2)}"
+        pw.cinter(args.prefix, date_tag, geo, variables, str(outdir) + "/")
+        out = outdir / f"{args.prefix}:{date_tag}"
+        if out.exists():
+            print(f"[OK] {out.name} (clim day-of-year idx {idx + 1})")
+            n += 1
+        else:
+            print(f"[ERROR] expected output not found: {out}", file=sys.stderr)
+        d += timedelta(days=1)
+
+    print(f"[OK] wrote {n} climatological {args.prefix}: file(s) -> {outdir}")
+    print(f"[INFO] Next: run init_atmosphere config_init_case=8 with "
+          f"config_sfc_prefix='{args.prefix}' (see README 'SST / sea-ice update').")
+    return 0 if n else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/oisst_to_intermediate.py b/usp-utils/pre_proc/real_data/oisst_to_intermediate.py
new file mode 100755
index 00000000..95d78946
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/oisst_to_intermediate.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+"""
+oisst_to_intermediate.py — convert NOAA OISST v2.1 daily NetCDF (from
+download_oisst.py) into WPS intermediate `SST:YYYY-MM-DD_HH` files for MPAS
+SST/sea-ice updates.
+
+Each daily OISST file becomes one intermediate file with the three 2D fields
+MPAS expects for an SST update (SST, SEAICE, LANDSEA). Feed the resulting
+`SST:` files to init_atmosphere with config_init_case=8 (see the README
+"SST / sea-ice update").
+
+Run inside the conda env that has xarray + pywinter:
+    conda run -n cgfd-usp-mpas python oisst_to_intermediate.py --indir <dir> \
+        [--start YYYY-MM-DD --end YYYY-MM-DD] [--hour 00]
+
+OISST specifics handled here:
+- sst is in degC -> converted to Kelvin (+273.15); land (missing sst) is filled
+  with a constant and masked by LANDSEA.
+- ice is a 0-1 fraction -> SEAICE.
+- LANDSEA is derived from the sst mask (1 = land, 0 = water), as OISST has no
+  separate mask variable in the daily file.
+- The grid is global 0.25 deg; latitude is flipped to ascending (S->N) if needed.
+"""
+
+import argparse
+import sys
+from datetime import datetime
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
+import pywinter.winter as pw
+
+LAND_SST_FILL_K = 273.15  # finite filler over land; ignored via LANDSEA
+
+
+def _date_from_name(p: Path):
+    """Parse YYYYMMDD out of oisst-avhrr-v02r01.YYYYMMDD[...].nc."""
+    digits = "".join(c for c in p.stem.split(".")[-1] if c.isdigit())[:8]
+    return datetime.strptime(digits, "%Y%m%d").date()
+
+
+def convert_one(path: Path, outdir: Path, prefix: str, hour: str) -> bool:
+    ds = xr.open_dataset(path)
+    sst = ds["sst"].isel(time=0, zlev=0)        # degC, NaN over land
+    ice = ds["ice"].isel(time=0, zlev=0)        # fraction, NaN over land
+    lat = sst["lat"].values
+    lon = sst["lon"].values
+    flip = lat[0] > lat[-1]
+    dlat = abs(float(lat[1] - lat[0]))
+    dlon = abs(float(lon[1] - lon[0]))
+    geo = pw.Geo0(float(min(lat)), float(min(lon)), dlat, dlon)
+
+    def arr2d(a):
+        a = np.asarray(a, dtype="float32")
+        return a[::-1, :] if flip else a
+
+    sst_v = sst.values
+    land = ~np.isfinite(sst_v)                  # land where sst is missing
+    sst_k = np.where(land, LAND_SST_FILL_K, sst_v + 273.15)
+    seaice = np.nan_to_num(ice.values, nan=0.0)
+    landsea = land.astype("float32")            # 1 = land, 0 = water
+
+    variables = [
+        pw.V2d("SST", arr2d(sst_k)),
+        pw.V2d("SEAICE", arr2d(seaice)),
+        pw.V2d("LANDSEA", arr2d(landsea)),
+    ]
+    date_tag = f"{_date_from_name(path).isoformat()}_{hour}"
+    pw.cinter(prefix, date_tag, geo, variables, str(outdir) + "/")
+    out = outdir / f"{prefix}:{date_tag}"
+    if out.exists():
+        print(f"[OK] {out.name} ({out.stat().st_size / 1e6:.1f} MB)")
+        return True
+    print(f"[ERROR] expected output not found: {out}", file=sys.stderr)
+    return False
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description="OISST NetCDF -> WPS intermediate SST: files.")
+    ap.add_argument("--indir", required=True, help="Directory with oisst-*.nc files")
+    ap.add_argument("--start", default=None, help="First date YYYY-MM-DD (default: all files)")
+    ap.add_argument("--end", default=None, help="Last date YYYY-MM-DD")
+    ap.add_argument("--outdir", default=None, help="Output dir (default: --indir)")
+    ap.add_argument("--hour", default="00", help="Hour tag HH for the daily field (default 00)")
+    ap.add_argument("--prefix", default="SST", help="Intermediate prefix (config_sfc_prefix)")
+    args = ap.parse_args()
+
+    indir = Path(args.indir)
+    outdir = Path(args.outdir) if args.outdir else indir
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    files = sorted(indir.glob("oisst-avhrr-v02r01.*.nc"))
+    if args.start:
+        d0 = datetime.strptime(args.start, "%Y-%m-%d").date()
+        d1 = datetime.strptime(args.end, "%Y-%m-%d").date() if args.end else d0
+        files = [f for f in files if d0 <= _date_from_name(f) <= d1]
+    if not files:
+        print(f"[ERROR] no OISST files found in {indir} for the requested range",
+              file=sys.stderr)
+        return 1
+
+    n = 0
+    for f in files:
+        print(f"[INFO] {f.name}")
+        if convert_one(f, outdir, args.prefix, args.hour.zfill(2)):
+            n += 1
+    print(f"[OK] wrote {n} {args.prefix}: intermediate file(s) -> {outdir}")
+    print(f"[INFO] Next: run init_atmosphere config_init_case=8 with "
+          f"config_sfc_prefix='{args.prefix}' (see README 'SST / sea-ice update').")
+    return 0 if n else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/usp-utils/pre_proc/real_data/prepare_era5.sh b/usp-utils/pre_proc/real_data/prepare_era5.sh
new file mode 100755
index 00000000..09a4af2a
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/prepare_era5.sh
@@ -0,0 +1,55 @@
+#!/usr/bin/env bash
+#
+# prepare_era5.sh — one-shot ERA5 atmospheric data for MPAS real-data init:
+# download pressure- and single-level GRIB (cdsapi) + convert to the WPS
+# intermediate file (GFS:YYYY-MM-DD_HH) that init_atmosphere reads.
+#
+# Use for the scientific/downscaling pipeline (ERA5 atmosphere; pair with OISST
+# for SST via prepare_oisst.sh). ERA5 covers 1940 -> present.
+#
+# Usage:
+#   ./prepare_era5.sh --date YYYY-MM-DD [--time 00] [--area "N W S E"]
+#                     [--outdir DIR] [--env ENV]
+#
+# Requires a CDS account + ~/.cdsapirc (https://cds.climate.copernicus.eu/how-to-api).
+# Runs the Python steps in the conda env with cdsapi + cfgrib + pywinter
+# (default: cgfd-usp-mpas). Prints config_nfglevels for namelist.init_atmosphere.
+
+set -euo pipefail
+SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+
+DATE="" TIME="00" AREA="" OUTDIR="" ENV="cgfd-usp-mpas"
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --date)   DATE="$2"; shift 2 ;;
+        --time)   TIME="$2"; shift 2 ;;
+        --area)   AREA="$2"; shift 2 ;;
+        --outdir) OUTDIR="$2"; shift 2 ;;
+        --env)    ENV="$2"; shift 2 ;;
+        -h|--help) sed -n '2,16p' "$0" | sed 's/^# \?//'; exit 0 ;;
+        *) echo "Unknown argument: $1" >&2; exit 2 ;;
+    esac
+done
+[ -n "$DATE" ] || { echo "ERROR: --date YYYY-MM-DD is required" >&2; exit 2; }
+
+run() { conda run -n "$ENV" python "$SCRIPT_DIR/$@"; }
+
+echo "[prepare_era5] download ${DATE} ${TIME}z"
+DL_ARGS=(--date "$DATE" --time "$TIME")
+[ -n "$OUTDIR" ] && DL_ARGS+=(--outdir "$OUTDIR")
+[ -n "$AREA" ] && DL_ARGS+=(--area $AREA)   # 4 tokens: N W S E
+run download_era5.py "${DL_ARGS[@]}"
+
+# locate the freshly downloaded GRIB files
+ymd=${DATE//-/}; hh=$(printf '%02d' "${TIME#0}")
+era5dir="${OUTDIR:-$( cd "$SCRIPT_DIR/../../.." && pwd )/met_data/era5/${ymd}${hh}}"
+pl="$era5dir/era5_pl_${ymd}${hh}.grib"
+sl="$era5dir/era5_sl_${ymd}${hh}.grib"
+[ -f "$pl" ] && [ -f "$sl" ] || { echo "ERROR: ERA5 GRIB not found in $era5dir" >&2; exit 1; }
+
+echo "[prepare_era5] convert -> WPS intermediate (GFS:)"
+run era5_to_intermediate.py --pl "$pl" --sl "$sl"
+
+echo "[prepare_era5] done. Set config_met_prefix='GFS' and the printed"
+echo "               config_nfglevels in namelist.init_atmosphere; link the"
+echo "               GFS:* intermediate (and the *.static.nc) into the run dir."
diff --git a/usp-utils/pre_proc/real_data/prepare_gfs.sh b/usp-utils/pre_proc/real_data/prepare_gfs.sh
new file mode 100755
index 00000000..9efb77af
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/prepare_gfs.sh
@@ -0,0 +1,81 @@
+#!/usr/bin/env bash
+#
+# prepare_gfs.sh — one-shot GFS met data for MPAS real-data init:
+# download (AWS S3 by default) + convert to WPS intermediate files.
+#
+# GFS feeds both products from a SINGLE download:
+#   atm  -> GFS:YYYY-MM-DD_HH  (atmosphere, for init_atmosphere case 7)
+#   sst  -> SST:YYYY-MM-DD_HH  (SST/sea-ice, for init_atmosphere case 8)
+#
+# Usage:
+#   ./prepare_gfs.sh --date YYYY-MM-DD [--cycle 00] [--fhour 0] [--res 0p25]
+#                    [--product atm|sst|both] [--source aws|nomads|rda]
+#                    [--outdir DIR] [--env ENV]
+#
+# --fhour is the forecast lead time in hours (0 = analysis). Use forecast hours
+# (e.g. 024 048 ...) of one cycle to get evolving SST over a forecast window;
+# the SST: file is auto-tagged with the valid date (cycle + fhour).
+#
+# Source coverage (GFS 0.25 deg): aws = 2021-03-23+, rda = 2015-01-15+,
+# nomads = last ~10 days. For older dates use ERA5 (prepare_era5.sh) for the
+# atmosphere and OISST (prepare_oisst.sh) for SST.
+#
+# Runs the Python steps in the conda env that has cfgrib + pywinter
+# (default: cgfd-usp-mpas). For atm it prints the config_nfglevels value to put
+# in namelist.init_atmosphere.
+
+set -euo pipefail
+SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+
+DATE="" CYCLE="00" FHOUR="0" RES="0p25" PRODUCT="atm" SOURCE="aws" OUTDIR="" ENV="cgfd-usp-mpas"
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --date)    DATE="$2"; shift 2 ;;
+        --cycle)   CYCLE="$2"; shift 2 ;;
+        --fhour)   FHOUR="$2"; shift 2 ;;
+        --res)     RES="$2"; shift 2 ;;
+        --product) PRODUCT="$2"; shift 2 ;;
+        --source)  SOURCE="$2"; shift 2 ;;
+        --outdir)  OUTDIR="$2"; shift 2 ;;
+        --env)     ENV="$2"; shift 2 ;;
+        -h|--help) sed -n '2,25p' "$0" | sed 's/^# \?//'; exit 0 ;;
+        *) echo "Unknown argument: $1" >&2; exit 2 ;;
+    esac
+done
+[ -n "$DATE" ] || { echo "ERROR: --date YYYY-MM-DD is required" >&2; exit 2; }
+fhh=$(printf 'f%03d' "$((10#$FHOUR))")   # base-10 so e.g. 024 is not read as octal
+case "$PRODUCT" in atm|sst|both) ;; *) echo "ERROR: --product must be atm|sst|both" >&2; exit 2 ;; esac
+
+run() { conda run -n "$ENV" python "$SCRIPT_DIR/$@"; }
+
+echo "[prepare_gfs] download ($SOURCE) ${DATE} ${CYCLE}z ${fhh} ${RES} [product=$PRODUCT]"
+DL_ARGS=(--date "$DATE" --cycle "$CYCLE" --fhour "$FHOUR" --res "$RES" --source "$SOURCE")
+# A pure SST product needs only the small surface subset (skip the 3D atmosphere).
+[ "$PRODUCT" = "sst" ] && DL_ARGS+=(--fields sst)
+[ -n "$OUTDIR" ] && DL_ARGS+=(--outdir "$OUTDIR")
+run download_gfs.py "${DL_ARGS[@]}"
+
+# locate the freshly downloaded GRIB (name depends on source, fhour and product)
+ymd=${DATE//-/}
+griddir="${OUTDIR:-$( cd "$SCRIPT_DIR/../../.." && pwd )/met_data/gfs/${ymd}${CYCLE}}"
+if [ "$SOURCE" = "rda" ]; then
+    grib="$griddir/gfs.0p25.${ymd}${CYCLE}.${fhh}.grib2"
+else
+    grib="$griddir/gfs.t${CYCLE}z.pgrb2.${RES}.${fhh}"
+    [ "$PRODUCT" = "sst" ] && grib="${grib}.sst"   # download_gfs.py --fields sst suffix
+fi
+[ -f "$grib" ] || { echo "ERROR: downloaded GRIB not found: $grib" >&2; exit 1; }
+
+if [ "$PRODUCT" = "atm" ] || [ "$PRODUCT" = "both" ]; then
+    echo "[prepare_gfs] convert -> atmosphere intermediate (GFS:)"
+    run gfs_to_intermediate.py --grib "$grib"
+fi
+if [ "$PRODUCT" = "sst" ] || [ "$PRODUCT" = "both" ]; then
+    echo "[prepare_gfs] convert -> SST/sea-ice intermediate (SST:)"
+    run gfs_sst_to_intermediate.py --grib "$grib"
+fi
+
+echo "[prepare_gfs] done."
+[ "$PRODUCT" != "sst" ] && echo "  atm: set config_met_prefix='GFS' + the printed config_nfglevels; link GFS:* + *.static.nc in the run dir."
+[ "$PRODUCT" != "atm" ] && echo "  sst: use config_init_case=8 with config_sfc_prefix='SST' (see README 'SST / sea-ice update')."
+exit 0
diff --git a/usp-utils/pre_proc/real_data/prepare_oisst.sh b/usp-utils/pre_proc/real_data/prepare_oisst.sh
new file mode 100755
index 00000000..5052fedf
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/prepare_oisst.sh
@@ -0,0 +1,55 @@
+#!/usr/bin/env bash
+#
+# prepare_oisst.sh — one-shot OISST SST/sea-ice data for MPAS SST update files.
+#
+# Two modes:
+#   (default)      observed daily OISST v2.1 (NCEI) over a date range  -> SST: files
+#   --climatology  OISST daily climatology (LTM 1991-2020, via OPeNDAP) -> SST: files
+#                  tagged with the target dates (for forecasts / typical-month runs,
+#                  where there is no observed SST for the dates).
+#
+# Usage:
+#   ./prepare_oisst.sh --start YYYY-MM-DD [--end YYYY-MM-DD] [--hour 00]
+#                      [--climatology] [--outdir DIR] [--env ENV]
+#
+# OISST observed covers 1981-09-01 -> present. Runs the Python steps in the conda
+# env with xarray + pywinter (default: cgfd-usp-mpas). See README "SST / sea-ice
+# update" for the init_atmosphere config_init_case=8 run.
+
+set -euo pipefail
+SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+
+START="" END="" HOUR="00" CLIM=0 OUTDIR="" ENV="cgfd-usp-mpas"
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --start)       START="$2"; shift 2 ;;
+        --end)         END="$2"; shift 2 ;;
+        --hour)        HOUR="$2"; shift 2 ;;
+        --climatology) CLIM=1; shift ;;
+        --outdir)      OUTDIR="$2"; shift 2 ;;
+        --env)         ENV="$2"; shift 2 ;;
+        -h|--help) sed -n '2,17p' "$0" | sed 's/^# \?//'; exit 0 ;;
+        *) echo "Unknown argument: $1" >&2; exit 2 ;;
+    esac
+done
+[ -n "$START" ] || { echo "ERROR: --start YYYY-MM-DD is required" >&2; exit 2; }
+END="${END:-$START}"
+
+run() { conda run -n "$ENV" python "$SCRIPT_DIR/$@"; }
+
+if [ "$CLIM" = "1" ]; then
+    echo "[prepare_oisst] climatology (LTM 1991-2020) ${START}..${END}"
+    CL_ARGS=(--start "$START" --end "$END" --hour "$HOUR")
+    [ -n "$OUTDIR" ] && CL_ARGS+=(--outdir "$OUTDIR")
+    run oisst_clim_to_intermediate.py "${CL_ARGS[@]}"
+    echo "[prepare_oisst] done (climatology)."
+else
+    oisstdir="${OUTDIR:-$( cd "$SCRIPT_DIR/../../.." && pwd )/met_data/oisst}"
+    echo "[prepare_oisst] download observed ${START}..${END}"
+    run download_oisst.py --start "$START" --end "$END" --outdir "$oisstdir"
+    echo "[prepare_oisst] convert -> WPS intermediate (SST:)"
+    run oisst_to_intermediate.py --indir "$oisstdir" --start "$START" --end "$END" --hour "$HOUR"
+    echo "[prepare_oisst] done. SST:* files in $oisstdir."
+fi
+echo "                Run init_atmosphere config_init_case=8 with"
+echo "                config_sfc_prefix='SST' (see README 'SST / sea-ice update')."
diff --git a/usp-utils/pre_proc/real_data/recipes/README.md b/usp-utils/pre_proc/real_data/recipes/README.md
new file mode 100644
index 00000000..08d18132
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/recipes/README.md
@@ -0,0 +1,25 @@
+# Data-preparation recipes (by use case)
+
+Each recipe lists exactly which data to download and convert (with the scripts in the
+parent directory) for one kind of MPAS run. They cover **data preparation only** — the
+mesh, the `init_atmosphere` runs (cases 7/8/9) and the model run are covered by the
+model tutorial (separate branch/directory).
+
+Pick by **what you simulate** (future vs past) and **mesh type** (global vs regional):
+
+| | Global mesh (no LBCs) | Regional mesh (needs LBCs) |
+|---|---|---|
+| **Operational forecast** (future) | [operational_global.md](operational_global.md) | [operational_regional.md](operational_regional.md) |
+| **Hindcast** (past) | [hindcast_global.md](hindcast_global.md) | [hindcast_regional.md](hindcast_regional.md) |
+
+Quick guide to the choices:
+- **Future (operational):** atmosphere from **GFS analysis**; future SST from **GFS
+  forecast** hours or **OISST climatology** (the future has no observed SST).
+- **Past (hindcast):** atmosphere from **ERA5** (or recent GFS analysis); SST from
+  **OISST observed** daily.
+- **Regional** additionally needs **lateral boundary** data at a regular cadence
+  (GFS forecast hours, or ERA5 at e.g. 6-hourly), which the LBC step (`init_atmosphere`
+  case 9) consumes.
+
+See the parent [README](../README.md) for the source/coverage tables, per-script options,
+the *SST / sea-ice update* mechanics, and the *Verify the init.nc* check.
diff --git a/usp-utils/pre_proc/real_data/recipes/hindcast_global.md b/usp-utils/pre_proc/real_data/recipes/hindcast_global.md
new file mode 100644
index 00000000..38fe8322
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/recipes/hindcast_global.md
@@ -0,0 +1,36 @@
+# Recipe — Hindcast, global mesh
+
+**Goal:** prepare the input *data* to simulate a **past** period on a **global** mesh
+(no lateral boundaries) — e.g. a global case study or a multi-month/year run. Scope here
+is data preparation only.
+
+| What | Source | When |
+|------|--------|------|
+| Atmosphere initial conditions | **ERA5** (1940→present) or **GFS analysis** (2015+) | once, at the run start |
+| SST / sea-ice | **OISST observed** daily | at the update cadence over the run |
+
+The past has observed data, so use a reanalysis (ERA5) for the atmosphere and the OISST
+daily analysis for SST. (For a recent date, 2021-03-23+, GFS analysis works too.)
+
+## 1. Atmosphere initial conditions
+```sh
+./prepare_era5.sh --date 2010-01-01 --time 00          # -> GFS:2010-01-01_00
+#   recent date alternative:
+# ./prepare_gfs.sh --date 2023-09-10 --cycle 00 --product atm
+```
+
+## 2. SST / sea-ice over the run (observed, daily)
+```sh
+./prepare_oisst.sh --start 2010-01-01 --end 2012-12-31   # one SST: per day
+```
+For a "typical" period rather than a specific year, use `--climatology` instead of the
+observed range. `config_fg_interval` (case-8 namelist) sets how often MPAS reads SST.
+
+## Hand-off to the model
+- `GFS:` → `init_atmosphere` **case 7** (global `*.static.nc`) → `x1.*.init.nc`.
+- `SST:` → `init_atmosphere` **case 8** → `x1.*.sfc_update.nc`; model `config_sst_update`
+  + the `sfc_update` stream. See README *SST / sea-ice update*.
+- **Verify the `init.nc`** (README *Verify the init.nc*) — this catches unit/interpolation
+  problems (e.g. the ERA5 path is newer; confirm `theta`/`rho` are physical).
+
+Global mesh ⇒ **no lateral boundary conditions**. Mesh and runs: separate tutorial.
diff --git a/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
new file mode 100644
index 00000000..4a8e52d2
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
@@ -0,0 +1,45 @@
+# Recipe — Hindcast / downscaling, regional (limited-area) mesh
+
+**Goal:** prepare the input *data* for a **downscaling / case study** on a **regional**
+mesh over a **past** period (the most common scientific use). Like the global hindcast
+but with **lateral boundary conditions (LBCs)** from ERA5 through the run. Scope here is
+data preparation only.
+
+| What | Source | When |
+|------|--------|------|
+| Atmosphere initial conditions | **ERA5** | once, at the run start |
+| Lateral boundary data | **ERA5** at the LBC cadence | every `N` h across the run |
+| SST / sea-ice | **OISST observed** daily | at the update cadence |
+
+## 1. Atmosphere initial conditions
+```sh
+./prepare_era5.sh --date 2014-09-10 --time 00          # -> GFS:2014-09-10_00
+```
+
+## 2. Boundary data through the run (the LBC step consumes these)
+One `GFS:` intermediate per boundary time = ERA5 at the LBC cadence (here 6-hourly):
+```sh
+for d in 2014-09-10 2014-09-11 2014-09-12; do
+    for t in 00 06 12 18; do
+        ./prepare_era5.sh --date $d --time $t          # GFS:<d>_<t>
+    done
+done
+```
+(Optionally pass `--area "N W S E"` to ERA5 to download only a box around the regional
+domain and keep the files small.)
+
+## 3. SST / sea-ice over the run (observed, daily)
+```sh
+./prepare_oisst.sh --start 2014-09-10 --end 2014-09-15
+```
+
+## Hand-off to the model
+- Initial conditions: `GFS:<init>` → `init_atmosphere` **case 7** (regional `*.static.nc`).
+- LBCs: the `GFS:` series at the boundary times → `init_atmosphere` **case 9** →
+  `lbc.*.nc`; model `config_apply_lbcs=.true.` + the `lbc` stream.
+- SST: `SST:` series → case 8 + `config_sst_update`.
+- **Verify the `init.nc`** (README *Verify the init.nc*).
+
+> **Out of scope here (separate tutorial / branch):** creating the limited-area mesh
+> (`create_region` from MPAS-Limited-Area) and running `init_atmosphere` case 7/9 and the
+> model. This directory only prepares the `GFS:`/`SST:` intermediates those steps read.
diff --git a/usp-utils/pre_proc/real_data/recipes/operational_global.md b/usp-utils/pre_proc/real_data/recipes/operational_global.md
new file mode 100644
index 00000000..37b78792
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/recipes/operational_global.md
@@ -0,0 +1,41 @@
+# Recipe — Operational forecast, global mesh
+
+**Goal:** prepare the input *data* for a real-time/near-real-time MPAS forecast on a
+**global** mesh (no lateral boundaries). Scope here is data preparation only; the mesh,
+the `init_atmosphere` runs and the model run are covered by the model tutorial.
+
+| What | Source | When |
+|------|--------|------|
+| Atmosphere initial conditions | latest **GFS analysis** (`--fhour 0`) | once, at the forecast start |
+| SST / sea-ice | GFS (init time) + GFS forecast or OISST climatology over the window | once, or per update step |
+
+GFS is a forecast model: the analysis (`--fhour 0`) is the best estimate at the cycle
+time; forecast hours (`--fhour > 0`) give the **future**. Use the latest available cycle.
+
+## 1. Atmosphere + initial SST (one GFS download)
+```sh
+./prepare_gfs.sh --date 2026-06-16 --cycle 00 --product both
+```
+Produces `GFS:2026-06-16_00` (atmosphere) and `SST:2026-06-16_00` (SST/sea-ice) under
+`met_data/gfs/2026061600/`.
+
+## 2. SST over the forecast window (only if you update SST)
+Short forecasts (≲1 week) usually keep SST **fixed** (`config_sst_update=.false.`) — then
+step 1 is enough. To let SST evolve (the future has no observed SST), use either:
+```sh
+# GFS forecast SST from the same cycle (each fhour -> SST: tagged with its valid date)
+for fh in 024 048 072 096 120; do
+    ./prepare_gfs.sh --date 2026-06-16 --cycle 00 --product sst --fhour $fh
+done
+# or the OISST day-of-year climatology for the forecast dates
+./prepare_oisst.sh --start 2026-06-17 --end 2026-06-21 --climatology
+```
+
+## Hand-off to the model
+- `GFS:` → `init_atmosphere` **case 7** (with the global `*.static.nc`) → `x1.*.init.nc`.
+- `SST:` → `init_atmosphere` **case 8** → `x1.*.sfc_update.nc`; model `config_sst_update`
+  + the `sfc_update` stream. See the README section *SST / sea-ice update*.
+- **Verify the `init.nc`** before running the model (README *Verify the init.nc*).
+
+Global mesh ⇒ **no lateral boundary conditions**. Mesh creation and the
+`init_atmosphere`/model runs are documented separately (other branch/tutorial).
diff --git a/usp-utils/pre_proc/real_data/recipes/operational_regional.md b/usp-utils/pre_proc/real_data/recipes/operational_regional.md
new file mode 100644
index 00000000..5540ff9a
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/recipes/operational_regional.md
@@ -0,0 +1,41 @@
+# Recipe — Operational forecast, regional (limited-area) mesh
+
+**Goal:** prepare the input *data* for an operational MPAS forecast on a **regional**
+(limited-area) mesh. A regional run needs, in addition to the global case, **lateral
+boundary conditions (LBCs)** through the forecast — and those come from GFS **forecast**
+hours of the cycle. Scope here is data preparation only.
+
+| What | Source | When |
+|------|--------|------|
+| Atmosphere initial conditions | **GFS analysis** (`--fhour 0`) | once, at the start |
+| Lateral boundary data | **GFS forecasts** of the same cycle, at the LBC cadence | every `N` h across the forecast |
+| SST / sea-ice | GFS or OISST climatology | once, or per update step |
+
+## 1. Atmosphere initial conditions
+```sh
+./prepare_gfs.sh --date 2026-06-16 --cycle 00 --product atm      # GFS:2026-06-16_00
+```
+
+## 2. Boundary data through the forecast (the LBC step consumes these)
+One `GFS:` intermediate per boundary-update time = forecast hours of the cycle at the LBC
+cadence (here every 6 h out to 120 h):
+```sh
+for fh in 006 012 018 024 030 036 042 048 ... 120; do
+    ./prepare_gfs.sh --date 2026-06-16 --cycle 00 --product atm --fhour $fh
+done   # each writes GFS:<valid date_hour> (cycle + fhour)
+```
+
+## 3. SST / sea-ice
+As in the global recipe (fixed at init, or GFS-forecast / OISST-climatology over the
+window): see [operational_global.md](operational_global.md) step 2.
+
+## Hand-off to the model
+- Initial conditions: `GFS:<init>` → `init_atmosphere` **case 7** (regional `*.static.nc`).
+- LBCs: the `GFS:` series at the boundary times → `init_atmosphere` **case 9** →
+  `lbc.*.nc`; model `config_apply_lbcs=.true.` + the `lbc` stream.
+- SST: as in the global recipe (case 8 + `config_sst_update`).
+- **Verify the `init.nc`** (README *Verify the init.nc*).
+
+> **Out of scope here (separate tutorial / branch):** creating the limited-area mesh
+> (`create_region` from MPAS-Limited-Area) and running `init_atmosphere` case 7/9 and the
+> model. This directory only prepares the `GFS:`/`SST:` intermediates those steps read.
diff --git a/usp-utils/pre_proc/static_fields/STATIC_FIELDS_GUIDE.md b/usp-utils/pre_proc/static_fields/STATIC_FIELDS_GUIDE.md
new file mode 100644
index 00000000..255d30a6
--- /dev/null
+++ b/usp-utils/pre_proc/static_fields/STATIC_FIELDS_GUIDE.md
@@ -0,0 +1,224 @@
+# Downloading meshes & terrain data for `init_atmosphere` static fields
+
+This guide covers the inputs and steps for the **static-field generation** stage
+of MPAS pre-processing — `init_atmosphere` with `config_init_case = 7`. That
+stage interpolates geographical/terrain data (topography, land use, soil,
+greenness, albedo, ...) onto an MPAS mesh, producing a *static* file that later
+real-data initialization and model runs build on.
+
+```
+   global mesh  (grid.nc + graph.info)            ┐
+                                                   ├─► init_atmosphere (case 7) ─► static file (e.g. x1.40962.init.nc)
+   MPAS geog static dataset (config_geog_data_path)┘
+```
+
+## Prerequisites
+
+- `init_atmosphere_model` compiled — see
+  [`../../install/README.md`](../../install/README.md).
+- The build environment sourced (so the MPICH `mpirun` is on `PATH`), e.g.:
+  `source ../../install/mpas_build_env.local.sh`. This path is just an
+  example — source whatever per-user build-env copy you created from the
+  `mpas_build_env.sh` template (you may have renamed or relocated it).
+
+---
+
+## Step 1 — Download the mesh and the MPAS geog data
+
+Use the helper (downloads, integrity-checks, **extracts into `WPS_GEOG/`**, and
+records `SHA256SUMS`):
+
+```sh
+./download_static_data.sh --mesh x1.40962 --geog mpas
+```
+
+- `--mesh NAME` — the global mesh (default `x1.40962`, ≈120 km). Nominal
+  resolutions of the quasi-uniform meshes (cell count → spacing):
+  `x1.2562`≈480 km, `x1.10242`≈240 km, `x1.40962`≈120 km, `x1.163842`≈60 km,
+  `x1.655362`≈30 km, `x1.2621442`≈15 km. The full list (incl.
+  variable-resolution meshes) is on the
+  [MPAS mesh page](https://mpas-dev.github.io/atmosphere/atmosphere_meshes.html).
+  The tarball contains `NAME.grid.nc`, `NAME.graph.info`, and pre-made
+  `NAME.graph.info.part.N` partitions (N = 2,4,6,8,...,128).
+- `--geog mpas` (default) — the **MPAS-curated** static bundle
+  ([`mpas_static.tar.bz2`](https://www2.mmm.ucar.edu/projects/mpas/site/downloads/static.html),
+  ~2.2 GB). **This bundle is complete for the standard real-data config**
+  (GMTED2010 topo + MODIS 30s land use + STATSGO soil + Noah-MP): it contains
+  every subdirectory `mpas_init_atm_static.F` reads for that config, including
+  `topo_gmted2010_30s`, `modis_landuse_20class_30s`, `soiltype_top_30s`, and
+  `soilgrids/{soilcomp,texture_layer1-4}`. The WRF `geog_high_res_mandatory`
+  bundle **lacks** `modis_landuse_20class_30s` and `soilgrids`, which is why
+  `mpas` is the default. Alternatives: `--geog high` (~2.6 GB, incomplete for
+  MPAS), `--geog low` (~150 MB), `--geog none` (skip).
+
+The mesh lands in `<dest>/grids/`; the geog tarball is cached in
+`<dest>/met_data/` and **extracted into `<dest>/WPS_GEOG/`** (override the base
+with `--dest DIR`). Point `config_geog_data_path` at that `WPS_GEOG/` directory.
+
+> **Extraction pitfall (the usual cause of failures).** The bundle carries a
+> single leading directory (`mpas_static/`). The helper strips it so datasets
+> land directly under `WPS_GEOG/<dataset>/`. If you extract the tarball by hand
+> *without* stripping it (`tar xjf mpas_static.tar.bz2 -C WPS_GEOG/`), you get
+> `WPS_GEOG/mpas_static/<dataset>/` and `init_atmosphere` fails with
+> `ERROR: Could not find an 'index' file in geotile directory ...`. Extract with
+> `--strip-components=1` (the helper does this for you):
+> ```sh
+> tar xjf met_data/mpas_static.tar.bz2 --strip-components=1 -C WPS_GEOG/
+> ```
+
+### Optional add-ons (only for non-default configs)
+
+The same MPAS downloads page offers add-ons that extract into the *same*
+`WPS_GEOG/`. Fetch them with `--optional` (comma list, or `all`):
+
+| `--optional` | File(s) | Needed when |
+|---|---|---|
+| `ugwp` | `topo_ugwp.tar.gz` + `ugwp_limb_tau.nc` | UGWP/GSL gravity-wave drag (`config_native_gwd_gsl_static = true` / UGWP suite) |
+| `15s`  | `modis_landuse_20class_15s.tar.bz2` | 15-arc-second (higher-res) land use |
+| `bnu`  | `bnu_soiltype_top.tar.bz2` | BNU soil category (`config_soilcat_data = 'BNU'`) |
+
+```sh
+./download_static_data.sh --mesh none --geog none --optional ugwp   # just the UGWP add-on
+```
+
+> `ugwp_limb_tau.nc` is **not** a geog dataset — it is a model **run-dir** input
+> (the `ugwp_ngw` stream). The helper downloads it into `met_data/`; copy it into
+> the `atmosphere_model` run directory if you enable the UGWP suite.
+
+**Download just one component** by setting the other to `none`:
+
+```sh
+./download_static_data.sh --mesh x1.10242 --geog none   # only a new mesh
+./download_static_data.sh --mesh none --geog mpas        # only the geog data
+```
+
+> **Provenance.** The script writes `SHA256SUMS` next to the downloads. The geog
+> bundle is large — download it once per machine and share it across projects
+> rather than re-downloading.
+
+---
+
+## Step 2 — Mesh partitioning (only if needed)
+
+`init_atmosphere` (and the model) decompose the mesh across MPI tasks using a
+`graph.info.part.N` file, where `N` is the number of tasks. The NCAR tarball
+already ships the common counts (2–128), selected via the namelist prefix:
+
+```
+config_block_decomp_file_prefix = 'x1.40962.graph.info.part.'
+```
+
+If you need a task count that is **not** provided, generate it with `gpmetis`
+(from METIS) in the `grids/` directory:
+
+```sh
+gpmetis x1.40962.graph.info 20      # -> x1.40962.graph.info.part.20
+```
+
+Running static generation on a single task needs no partition file at all.
+
+---
+
+## Step 3 — Configure `namelist.init_atmosphere` and streams
+
+Create a run directory (links the executable + copies default namelists/streams):
+
+```sh
+"$MPAS_ROOT/testing_and_setup/atmosphere/setup_run_dir.py" "$MPAS_ROOT/runs/static_x1.40962"
+cd "$MPAS_ROOT/runs/static_x1.40962"
+ln -s "$MPAS_ROOT/grids/x1.40962.grid.nc" .
+ln -s "$MPAS_ROOT/grids/x1.40962.graph.info.part."* .
+```
+
+Key `namelist.init_atmosphere` options for the static case:
+
+```
+&nhyd_model
+    config_init_case = 7                 ! static-field generation
+/
+&data_sources
+    config_geog_data_path = '/path/to/WPS_GEOG/'   ! from Step 1, trailing '/'
+    config_landuse_data   = 'MODIFIED_IGBP_MODIS_NOAH'
+    config_soilcat_data   = 'STATSGO'
+    config_topo_data      = 'GMTED2010'
+/
+&preproc_stages
+    config_static_interp         = true
+    config_native_gwd_static     = true   ! conventional GWDO (computed from GMTED2010)
+    config_native_gwd_gsl_static = false  ! UGWP/GSL drag — needs --optional ugwp; see note
+/
+&decomposition
+    config_block_decomp_file_prefix = 'x1.40962.graph.info.part.'
+/
+```
+
+> **Noah-MP soil composition is ON by default (hidden option).**
+> `config_noahmp_static` has `default_value="true"` but `in_defaults="false"` in
+> the Registry, so it is **not written** into the generated namelist yet is
+> active. It makes the static stage interpolate `soilgrids/soilcomp` and
+> `soilgrids/texture_layer1-4` — hence `soilgrids/` **must** be present under
+> `WPS_GEOG/`. The `mpas` bundle provides it; the WRF `high` bundle does not. If
+> you deliberately do not want these fields, set `config_noahmp_static = false`
+> explicitly.
+
+> **UGWP/GSL gravity-wave-drag static (optional).** Setting
+> `config_native_gwd_gsl_static = true` makes init_atmosphere also build the GSL
+> orographic-drag fields (`oro_data_ss`/`oro_data_ls`), which need the
+> `topo_ugwp_*` raw-topography datasets. Those are **MPAS-specific** and ship as
+> the **optional** `topo_ugwp.tar.gz` on the downloads page — fetch them with
+> `./download_static_data.sh --optional ugwp` (and copy `ugwp_limb_tau.nc` into
+> the model run dir). If a `topo_ugwp_*` tile is missing you get non-fatal
+> `ERROR: Error reading topography tile ...` and the run still finishes, but the
+> fields are empty. For standard runs keep this **false** (the conventional GWDO
+> above is sufficient).
+
+In `streams.init_atmosphere`, the `input` stream reads `x1.40962.grid.nc` and
+the `output` stream writes the static file (default `x1.40962.init.nc`; you may
+rename the template to `x1.40962.static.nc` to make its role explicit).
+
+---
+
+## Step 4 — Run
+
+```sh
+# single task (fine for static generation on moderate meshes):
+./init_atmosphere_model
+
+# or multiple tasks (needs a matching graph.info.part.N from Step 2):
+mpirun -np 4 ./init_atmosphere_model
+```
+
+Check `log.init_atmosphere.0000.out` for `Finished running the init_atmosphere
+core`. The output static file appears in the run directory.
+
+---
+
+## Output and next steps
+
+The static file (e.g. `x1.40962.init.nc`) carries the mesh **plus** the
+interpolated terrain/land-surface fields. It is the input for the next
+pre-processing stage — real-data initialization (a different `config_init_case`
+using meteorological data in `met_data/`) — which produces the actual initial
+conditions for `atmosphere_model`. Those stages are covered separately.
+
+---
+
+## Troubleshooting
+
+- **`ERROR: Could not find an 'index' file in geotile directory .../<dataset>/`**
+  — the dataset directory is missing or nested one level too deep. Almost always
+  the bundle was extracted **without** stripping its `mpas_static/` top
+  directory (see the extraction pitfall in Step 1), so the data sits in
+  `WPS_GEOG/mpas_static/<dataset>/`. Re-extract with `--strip-components=1`
+  (or just rerun `download_static_data.sh`). If it specifically names
+  `soilgrids/soilcomp/`, your geog set lacks `soilgrids` — use `--geog mpas`,
+  not the WRF `high`/`low` bundles.
+- **`init_atmosphere` can't find the geog data** — `config_geog_data_path` must
+  point to the extracted `WPS_GEOG/` directory and end with a trailing `/`.
+- **Decomposition / partition error** — the `graph.info.part.N` file for your
+  task count is missing; download a mesh that includes it, run on 1 task, or
+  build it with `gpmetis` (Step 2).
+- **`config_topo_data`/`config_landuse_data`/`config_soilcat_data` mismatch** —
+  those names must match subdirectories present under `WPS_GEOG/`; the `mpas`
+  bundle provides GMTED2010 topo, `MODIFIED_IGBP_MODIS_NOAH` land use, and
+  `STATSGO` soil.
diff --git a/usp-utils/pre_proc/static_fields/download_static_data.sh b/usp-utils/pre_proc/static_fields/download_static_data.sh
new file mode 100755
index 00000000..e01c5678
--- /dev/null
+++ b/usp-utils/pre_proc/static_fields/download_static_data.sh
@@ -0,0 +1,254 @@
+#!/usr/bin/env bash
+#
+# download_static_data.sh
+#
+# Download the inputs needed for MPAS init_atmosphere static-field generation
+# (config_init_case = 7):
+#   1. a global mesh from the NCAR MPAS-Atmosphere mesh repository
+#      (grid.nc + graph.info + pre-made graph.info.part.N partitions);
+#   2. the MPAS geographical ("static"/terrain) dataset that init_atmosphere
+#      interpolates onto the mesh (topography, land use, soil, greenness, ...).
+#
+# Sources:
+#   meshes        -> https://www2.mmm.ucar.edu/projects/mpas/atmosphere_meshes/
+#   static (mpas) -> https://www2.mmm.ucar.edu/projects/mpas/site/downloads/static.html
+#   WPS geog      -> https://www2.mmm.ucar.edu/wrf/src/wps_files/
+#
+# The default MPAS bundle (mpas_static.tar.bz2) is COMPLETE for the standard
+# real-data config (GMTED2010 topo + MODIS 30s land use + STATSGO soil +
+# Noah-MP). It contains every subdirectory mpas_init_atm_static.F reads for that
+# config, including soilgrids/{soilcomp,texture_layer1-4} (required because
+# config_noahmp_static defaults to TRUE even though it is not written into the
+# generated namelist) and modis_landuse_20class_30s. The WRF geog_high_res
+# bundle does NOT contain those, which is why 'mpas' is the default here.
+#
+# Optional add-ons (page above) — only for non-default configs:
+#   ugwp -> topo_ugwp.tar.gz + ugwp_limb_tau.nc  (UGWP/GSL gravity-wave drag,
+#           config_native_gwd_gsl_static = true / UGWP suite)
+#   15s  -> modis_landuse_20class_15s.tar.bz2     (15-arc-second land use)
+#   bnu  -> bnu_soiltype_top.tar.bz2              (BNU soil category, config_soilcat_data='BNU')
+#
+# IMPORTANT (extraction): the bundles carry a single leading directory
+# (e.g. mpas_static/ or WPS_GEOG/). This script strips it so datasets land
+# directly under WPS_GEOG/<dataset>, matching config_geog_data_path. Extracting
+# the tarball "as is" produces WPS_GEOG/mpas_static/<dataset> and init_atmosphere
+# then fails with "Could not find an 'index' file" — that is the usual pitfall.
+#
+# Usage:
+#   ./download_static_data.sh [--mesh NAME] [--geog mpas|high|low|none]
+#                             [--optional LIST|all] [--dest DIR] [--no-extract]
+#
+#   --mesh NAME     Mesh to download (default: x1.40962), or 'none' to skip.
+#                   See the mesh page for the full list (x1.10242, x1.40962,
+#                   x1.163842, x1.655362, x1.2621442, variable-res x4/x5..., etc.).
+#   --geog SET      Geographical data set:
+#                     mpas (~2.2 GB, DEFAULT) — MPAS-curated bundle. Complete.
+#                     high (~2.6 GB) — WRF 'geog_high_res_mandatory' (INCOMPLETE
+#                          for MPAS: lacks modis_landuse_20class_30s, soilgrids).
+#                     low  (~150 MB) — WRF low-res mandatory.
+#                     none — skip the geog.
+#   --optional LIST Comma-separated add-ons to also fetch into WPS_GEOG:
+#                   any of ugwp,15s,bnu  (or 'all'). Default: none.
+#   --dest DIR      Base directory. Default: repo root above usp-utils.
+#                   Tarballs cache into <dest>/met_data; datasets extract into
+#                   <dest>/WPS_GEOG (= config_geog_data_path).
+#   --no-extract    Download/verify the archives but do not extract them.
+#
+# Download just one component by setting the other to 'none', e.g.:
+#   ./download_static_data.sh --mesh x1.10242 --geog none    # only a new mesh
+#   ./download_static_data.sh --mesh none --geog mpas         # only the geog data
+#   ./download_static_data.sh --mesh none --geog none --optional ugwp  # only add-ons
+#
+# Idempotent: existing, valid downloads are skipped; each archive's extraction is
+# marked in WPS_GEOG/.extracted_<archive>; checksums recorded in SHA256SUMS.
+
+set -uo pipefail
+
+CYAN='\033[0;36m'; YELLOW='\033[1;33m'; RED='\033[0;31m'; GREEN='\033[0;32m'; NC='\033[0m'
+INFO="${CYAN}[INFO]${NC}"; WARN="${YELLOW}[WARN]${NC}"; ERR="${RED}[ERROR]${NC}"; OK="${GREEN}[OK]${NC}"
+
+SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+DEFAULT_ROOT=$( cd -- "$SCRIPT_DIR/../../.." &> /dev/null && pwd )   # repo root
+
+# --- Defaults / argument parsing ---------------------------------------------
+MESH="x1.40962"
+GEOG="mpas"
+OPTIONAL=""
+DEST="$DEFAULT_ROOT"
+EXTRACT=1
+
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --mesh)       MESH="$2"; shift 2 ;;
+        --geog)       GEOG="$2"; shift 2 ;;
+        --optional)   OPTIONAL="$2"; shift 2 ;;
+        --dest)       DEST="$2"; shift 2 ;;
+        --no-extract) EXTRACT=0; shift ;;
+        -h|--help)    grep '^#' "$0" | sed 's/^# \?//'; exit 0 ;;
+        *) echo -e "${ERR} Unknown argument: $1"; exit 2 ;;
+    esac
+done
+
+MESH_BASE="https://www2.mmm.ucar.edu/projects/mpas/atmosphere_meshes"
+MPAS_BASE="https://www2.mmm.ucar.edu/projects/mpas"
+GEOG_BASE="https://www2.mmm.ucar.edu/wrf/src/wps_files"
+DL_DIR="$DEST/met_data"       # tarball download cache
+GEOG_DIR="$DEST/WPS_GEOG"     # extraction target == config_geog_data_path
+
+case "$GEOG" in
+    mpas) GEOG_FILE="mpas_static.tar.bz2";            GEOG_URL="$MPAS_BASE/mpas_static.tar.bz2" ;;
+    high) GEOG_FILE="geog_high_res_mandatory.tar.gz"; GEOG_URL="$GEOG_BASE/geog_high_res_mandatory.tar.gz" ;;
+    low)  GEOG_FILE="geog_low_res_mandatory.tar.gz";  GEOG_URL="$GEOG_BASE/geog_low_res_mandatory.tar.gz" ;;
+    none) GEOG_FILE=""; GEOG_URL="" ;;
+    *) echo -e "${ERR} --geog must be mpas, high, low, or none"; exit 2 ;;
+esac
+
+# expand --optional all
+if [ "$OPTIONAL" = "all" ]; then OPTIONAL="ugwp,15s,bnu"; fi
+
+if [ "$MESH" = "none" ] && [ -z "$GEOG_FILE" ] && [ -z "$OPTIONAL" ]; then
+    echo -e "${ERR} Nothing to do: --mesh none, --geog none, no --optional."
+    exit 2
+fi
+
+# --- Helpers ------------------------------------------------------------------
+# integrity test by archive type (gzip or bzip2); plain files (.nc) pass through.
+verify_archive() {
+    case "$1" in
+        *.bz2) bzip2 -t "$1" 2>/dev/null ;;
+        *.gz)  gzip  -t "$1" 2>/dev/null ;;
+        *)     [ -s "$1" ] ;;
+    esac
+}
+
+# fetch URL -> dest file (idempotent + integrity-checked). Returns non-zero on failure.
+fetch() {
+    local url="$1" dest="$2" name; name="$(basename "$dest")"
+    if [ -f "$dest" ] && verify_archive "$dest"; then
+        echo -e "${OK} ${name} already present and valid — skipping"
+        return 0
+    fi
+    echo -e "${INFO} Downloading ${name}"
+    if ! curl -fL -C - --retry 3 --retry-delay 5 -o "$dest" "$url"; then
+        echo -e "${ERR} Download failed: ${url}"; rm -f "$dest"; return 1
+    fi
+    if ! verify_archive "$dest"; then
+        echo -e "${ERR} Integrity check failed (corrupt/truncated): ${name}"; rm -f "$dest"; return 1
+    fi
+    echo -e "${OK} ${name}"
+}
+
+# extract an archive into WPS_GEOG, stripping a single leading directory if the
+# archive has exactly one top-level component (mpas_static/, WPS_GEOG/, ...).
+# Idempotent via a per-archive marker file. Returns non-zero on failure.
+extract_into_geog() {
+    local tar="$1" name; name="$(basename "$tar")"
+    local marker="$GEOG_DIR/.extracted_${name}"
+    mkdir -p "$GEOG_DIR"
+    if [ -f "$marker" ]; then
+        echo -e "${OK} ${name} already extracted into WPS_GEOG — skipping"
+        return 0
+    fi
+    local tops n
+    tops=$(tar tf "$tar" 2>/dev/null | sed 's#^\./##; /^$/d' | cut -d/ -f1 | sort -u | grep -c .)
+    n="$tops"
+    echo -e "${INFO} Extracting ${name} into WPS_GEOG (this takes a while)"
+    if [ "$n" -eq 1 ]; then
+        tar xf "$tar" --strip-components=1 -C "$GEOG_DIR" || return 1
+    else
+        tar xf "$tar" -C "$GEOG_DIR" || return 1
+    fi
+    touch "$marker"
+    echo -e "${OK} extracted ${name} -> WPS_GEOG/"
+}
+
+record_sha() { ( cd "$1" && sha256sum "$2" 2>/dev/null >> SHA256SUMS && sort -u -o SHA256SUMS SHA256SUMS ); }
+
+failures=0
+
+# --- Mesh ---------------------------------------------------------------------
+if [ "$MESH" != "none" ]; then
+    MESH_DIR="$DEST/grids"
+    mkdir -p "$MESH_DIR"
+    echo -e "${INFO} Mesh '${MESH}' -> ${MESH_DIR}"
+    mesh_tar="$MESH_DIR/${MESH}.tar.gz"
+    if fetch "${MESH_BASE}/${MESH}.tar.gz" "$mesh_tar"; then
+        if [ "$EXTRACT" -eq 1 ]; then
+            if [ -f "$MESH_DIR/${MESH}.grid.nc" ]; then
+                echo -e "${OK} ${MESH}.grid.nc already extracted"
+            else
+                echo -e "${INFO} Extracting ${MESH}.tar.gz"
+                tar xzf "$mesh_tar" -C "$MESH_DIR" && echo -e "${OK} extracted ${MESH} (grid.nc + graph.info[.part.N])"
+            fi
+        fi
+        record_sha "$MESH_DIR" "${MESH}.tar.gz"
+    else
+        failures=$((failures+1))
+    fi
+else
+    echo -e "${INFO} Skipping mesh download (--mesh none)"
+fi
+
+# --- geog (default bundle) ----------------------------------------------------
+if [ -n "$GEOG_FILE" ]; then
+    mkdir -p "$DL_DIR"
+    echo -e "${INFO} geog (${GEOG}) -> extract into ${GEOG_DIR}  (large: ~2.2 GB mpas / ~2.6 GB high)"
+    geog_tar="$DL_DIR/${GEOG_FILE}"
+    if fetch "${GEOG_URL}" "$geog_tar"; then
+        [ "$EXTRACT" -eq 1 ] && { extract_into_geog "$geog_tar" || failures=$((failures+1)); }
+        record_sha "$DL_DIR" "$GEOG_FILE"
+    else
+        failures=$((failures+1))
+    fi
+fi
+
+# --- optional add-ons ---------------------------------------------------------
+if [ -n "$OPTIONAL" ]; then
+    mkdir -p "$DL_DIR"
+    IFS=',' read -ra OPTS <<< "$OPTIONAL"
+    for opt in "${OPTS[@]}"; do
+        case "$opt" in
+            ugwp)
+                # terrain stats (geog) + non-grid GW flux file (a run-dir input)
+                for f in topo_ugwp.tar.gz; do
+                    if fetch "$MPAS_BASE/$f" "$DL_DIR/$f"; then
+                        [ "$EXTRACT" -eq 1 ] && { extract_into_geog "$DL_DIR/$f" || failures=$((failures+1)); }
+                        record_sha "$DL_DIR" "$f"
+                    else failures=$((failures+1)); fi
+                done
+                if fetch "$MPAS_BASE/ugwp_limb_tau.nc" "$DL_DIR/ugwp_limb_tau.nc"; then
+                    record_sha "$DL_DIR" "ugwp_limb_tau.nc"
+                    echo -e "${INFO} ugwp_limb_tau.nc is a RUN-DIR input (ugwp_ngw stream), not a geog dataset: copy it into the atmosphere run dir. Kept in ${DL_DIR}."
+                else failures=$((failures+1)); fi
+                ;;
+            15s)
+                f="modis_landuse_20class_15s.tar.bz2"
+                if fetch "$MPAS_BASE/$f" "$DL_DIR/$f"; then
+                    [ "$EXTRACT" -eq 1 ] && { extract_into_geog "$DL_DIR/$f" || failures=$((failures+1)); }
+                    record_sha "$DL_DIR" "$f"
+                else failures=$((failures+1)); fi
+                ;;
+            bnu)
+                f="bnu_soiltype_top.tar.bz2"
+                if fetch "$MPAS_BASE/$f" "$DL_DIR/$f"; then
+                    [ "$EXTRACT" -eq 1 ] && { extract_into_geog "$DL_DIR/$f" || failures=$((failures+1)); }
+                    record_sha "$DL_DIR" "$f"
+                else failures=$((failures+1)); fi
+                ;;
+            *) echo -e "${WARN} unknown --optional item '${opt}' (expected ugwp,15s,bnu) — skipping" ;;
+        esac
+    done
+fi
+
+echo ""
+if [ "$failures" -eq 0 ]; then
+    [ "$MESH" != "none" ]                      && echo -e "${OK} Mesh in ${DEST}/grids"
+    { [ -n "$GEOG_FILE" ] || [ -n "$OPTIONAL" ]; } && echo -e "${OK} geog datasets in ${GEOG_DIR}"
+    [ -n "$GEOG_FILE" ] && [ "$EXTRACT" -eq 1 ] && \
+        echo -e "${INFO} Set config_geog_data_path = '${GEOG_DIR}/' in namelist.init_atmosphere"
+    echo -e "${INFO} Next: see usp-utils/pre_proc/static_fields/STATIC_FIELDS_GUIDE.md to run init_atmosphere (case 7)."
+    exit 0
+else
+    echo -e "${ERR} ${failures} step(s) failed — see messages above."
+    exit 1
+fi

From 834b90c16c03929cb4eded53e95be899a9c42252 Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Sun, 21 Jun 2026 15:14:27 +0000
Subject: [PATCH 5/7] pre_proc/real_data: name ERA5 intermediates ERA5: and
 expand regional downscaling docs

Name and docs:
- era5_to_intermediate.py defaults --prefix to 'ERA5' (was 'GFS'); the ERA5
  intermediate is now ERA5:YYYY-MM-DD_HH, set config_met_prefix='ERA5'. GFS keeps
  'GFS:'. prepare_era5.sh and the docs follow suit; GFS-source mentions unchanged.

Regional downscaling docs:
- hindcast_regional.md: explain --area as a download box (not the simulated
  region, which is the mesh) sized to mesh + margin; add a window-based recipe
  (START/END) covering a single month, a year, or multiple years, with the
  6-hourly LBC loop + closing boundary and daily OISST SST.
- README.md: clarify the --area bullet and link the recipe.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 usp-utils/pre_proc/real_data/README.md        | 38 ++++----
 .../real_data/era5_to_intermediate.py         | 12 +--
 usp-utils/pre_proc/real_data/prepare_era5.sh  |  8 +-
 .../real_data/recipes/hindcast_global.md      |  5 +-
 .../real_data/recipes/hindcast_regional.md    | 88 +++++++++++++++----
 5 files changed, 105 insertions(+), 46 deletions(-)

diff --git a/usp-utils/pre_proc/real_data/README.md b/usp-utils/pre_proc/real_data/README.md
index 7b8781b5..bca71f83 100644
--- a/usp-utils/pre_proc/real_data/README.md
+++ b/usp-utils/pre_proc/real_data/README.md
@@ -1,18 +1,18 @@
 # `real_data/` — GFS / ERA5 → WPS intermediate met files for MPAS
 
 Scripts to fetch meteorological data and turn it into the **WPS intermediate-format**
-file (`GFS:YYYY-MM-DD_HH`) that `init_atmosphere` reads to build the real-data
-`init.nc`. MPAS has no built-in ungrib, so this is done in Python with **pywinter**
-— no need to build WPS/ungrib.
+file (`GFS:` or `ERA5:` `YYYY-MM-DD_HH`) that `init_atmosphere` reads to build the
+real-data `init.nc`. MPAS has no built-in ungrib, so this is done in Python with
+**pywinter** — no need to build WPS/ungrib.
 
 ```
-download_gfs.py   →  GFS GRIB2    →  gfs_to_intermediate.py   ┐
-                                                              ├→  GFS:YYYY-MM-DD_HH
-download_era5.py  →  ERA5 GRIB    →  era5_to_intermediate.py  ┘
+download_gfs.py   →  GFS GRIB2    →  gfs_to_intermediate.py   →  GFS:YYYY-MM-DD_HH
+download_era5.py  →  ERA5 GRIB    →  era5_to_intermediate.py  →  ERA5:YYYY-MM-DD_HH
 ```
 
-The intermediate file keeps the `GFS:` prefix (= `config_met_prefix='GFS'`)
-regardless of the source, so the `init_atmosphere` setup is identical either way.
+Each source writes its own prefix — GFS → `GFS:`, ERA5 → `ERA5:` — which is the
+`config_met_prefix` you set in `namelist.init_atmosphere`. The field names and the
+rest of the setup are identical either way; only the prefix differs.
 This directory only produces that **met intermediate file** — running
 `init_atmosphere` with it (namelist, streams, linking `*.static.nc`, etc.) is covered
 by the model usage tutorial, not repeated here. Run these after the static fields
@@ -20,8 +20,9 @@ are ready (see `../static_fields/`).
 
 ## Which source? (data availability)
 
-Two products feed init_atmosphere: the **atmosphere** (3D state → `GFS:` intermediate,
-init case 7) and the **SST/sea-ice** lower boundary (→ `SST:` intermediate, init case 8,
+Two products feed init_atmosphere: the **atmosphere** (3D state → `GFS:`/`ERA5:`
+intermediate, init case 7) and the **SST/sea-ice** lower boundary (→ `SST:`
+intermediate, init case 8,
 see "SST / sea-ice update").
 
 **Atmosphere:**
@@ -78,7 +79,7 @@ One-shot wrappers (download + convert in one call, inside the conda env):
 | Wrapper | Covers |
 |---------|--------|
 | `prepare_gfs.sh` | GFS atmosphere and/or SST — `--product atm` (default), `sst`, or `both` (one download, both products). |
-| `prepare_era5.sh` | ERA5 atmosphere (download pl+sl → `GFS:` intermediate). |
+| `prepare_era5.sh` | ERA5 atmosphere (download pl+sl → `ERA5:` intermediate). |
 | `prepare_oisst.sh` | OISST SST over a date range — observed (default) or `--climatology` (→ daily `SST:` intermediates). |
 
 ## Dependencies
@@ -131,16 +132,21 @@ Key options:
   `--fields {full,sst}` (`sst` = small surface-only download for SST update files),
   `--outdir`.
 - `download_era5.py`: `--date YYYY-MM-DD`, `--time HH`, `--area N W S E` (default
-  global), `--outdir`.
+  global), `--outdir`. **`--area` subsets only the *download*, not the simulated
+  region** (that is the mesh): use it for a regional run to keep a long LBC series
+  small, sized as the regional mesh extent **+ a margin** (order `N W S E`, degrees,
+  lon −180..180; W/E negative over Brazil). For a global mesh, omit it. See
+  [hindcast_regional](recipes/hindcast_regional.md) for a worked example (single
+  month to multiple years).
 - `download_oisst.py` / `oisst_to_intermediate.py`: `--start`/`--end` date range,
   `--hour HH`, `--outdir`.
 - `oisst_clim_to_intermediate.py` (climatology): `--start`/`--end` target dates,
   `--hour HH`, `--outdir` (no download — reads the LTM over OPeNDAP).
-- converters: `--outdir`, `--prefix` (default `GFS` for atmosphere, `SST` for SST);
-  GFS uses `--grib`, ERA5 uses `--pl`/`--sl`.
+- converters: `--outdir`, `--prefix` (default `GFS` for the GFS converter, `ERA5`
+  for the ERA5 converter, `SST` for SST); GFS uses `--grib`, ERA5 uses `--pl`/`--sl`.
 
-The converter sets the date tag from the data's valid time, so the output name
-(`GFS:YYYY-MM-DD_HH`) is what `config_start_time` must match.
+The converter sets the date tag from the data's valid time, so the date part of the
+output name (`<prefix>:YYYY-MM-DD_HH`) is what `config_start_time` must match.
 
 ## Verify the init.nc (recommended after every run)
 
diff --git a/usp-utils/pre_proc/real_data/era5_to_intermediate.py b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
index 76d73273..08d05367 100755
--- a/usp-utils/pre_proc/real_data/era5_to_intermediate.py
+++ b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
@@ -1,11 +1,13 @@
 #!/usr/bin/env python3
 """
 era5_to_intermediate.py — convert ERA5 GRIB files (from download_era5.py) into a
-WPS intermediate-format file (e.g. GFS:2014-09-10_00) that MPAS init_atmosphere
-reads for real-data initialization (config_met_prefix='GFS').
+WPS intermediate-format file (e.g. ERA5:2014-09-10_00) that MPAS init_atmosphere
+reads for real-data initialization (config_met_prefix='ERA5').
 
 This is the ERA5 counterpart of gfs_to_intermediate.py (same pywinter output,
-same intermediate field names), for dates the GFS archive does not cover.
+same intermediate field names), for dates the GFS archive does not cover. The
+only difference is the file prefix: ERA5 writes 'ERA5:', GFS writes 'GFS:', so
+set config_met_prefix to match the source you used.
 
 Run inside the conda env that has cfgrib + pywinter:
     conda run -n cgfd-usp-mpas python era5_to_intermediate.py \
@@ -80,8 +82,8 @@ def main() -> int:
     ap.add_argument("--sl", required=True, help="ERA5 single-level GRIB file")
     ap.add_argument("--outdir", default=None,
                     help="Output dir (default: the --pl file's directory)")
-    ap.add_argument("--prefix", default="GFS",
-                    help="Intermediate prefix (config_met_prefix; keep 'GFS')")
+    ap.add_argument("--prefix", default="ERA5",
+                    help="Intermediate prefix (config_met_prefix; default 'ERA5')")
     args = ap.parse_args()
 
     pl, sl = Path(args.pl), Path(args.sl)
diff --git a/usp-utils/pre_proc/real_data/prepare_era5.sh b/usp-utils/pre_proc/real_data/prepare_era5.sh
index 09a4af2a..d3f3cff5 100755
--- a/usp-utils/pre_proc/real_data/prepare_era5.sh
+++ b/usp-utils/pre_proc/real_data/prepare_era5.sh
@@ -2,7 +2,7 @@
 #
 # prepare_era5.sh — one-shot ERA5 atmospheric data for MPAS real-data init:
 # download pressure- and single-level GRIB (cdsapi) + convert to the WPS
-# intermediate file (GFS:YYYY-MM-DD_HH) that init_atmosphere reads.
+# intermediate file (ERA5:YYYY-MM-DD_HH) that init_atmosphere reads.
 #
 # Use for the scientific/downscaling pipeline (ERA5 atmosphere; pair with OISST
 # for SST via prepare_oisst.sh). ERA5 covers 1940 -> present.
@@ -47,9 +47,9 @@ pl="$era5dir/era5_pl_${ymd}${hh}.grib"
 sl="$era5dir/era5_sl_${ymd}${hh}.grib"
 [ -f "$pl" ] && [ -f "$sl" ] || { echo "ERROR: ERA5 GRIB not found in $era5dir" >&2; exit 1; }
 
-echo "[prepare_era5] convert -> WPS intermediate (GFS:)"
+echo "[prepare_era5] convert -> WPS intermediate (ERA5:)"
 run era5_to_intermediate.py --pl "$pl" --sl "$sl"
 
-echo "[prepare_era5] done. Set config_met_prefix='GFS' and the printed"
+echo "[prepare_era5] done. Set config_met_prefix='ERA5' and the printed"
 echo "               config_nfglevels in namelist.init_atmosphere; link the"
-echo "               GFS:* intermediate (and the *.static.nc) into the run dir."
+echo "               ERA5:* intermediate (and the *.static.nc) into the run dir."
diff --git a/usp-utils/pre_proc/real_data/recipes/hindcast_global.md b/usp-utils/pre_proc/real_data/recipes/hindcast_global.md
index 38fe8322..101be516 100644
--- a/usp-utils/pre_proc/real_data/recipes/hindcast_global.md
+++ b/usp-utils/pre_proc/real_data/recipes/hindcast_global.md
@@ -14,7 +14,7 @@ daily analysis for SST. (For a recent date, 2021-03-23+, GFS analysis works too.
 
 ## 1. Atmosphere initial conditions
 ```sh
-./prepare_era5.sh --date 2010-01-01 --time 00          # -> GFS:2010-01-01_00
+./prepare_era5.sh --date 2010-01-01 --time 00          # -> ERA5:2010-01-01_00
 #   recent date alternative:
 # ./prepare_gfs.sh --date 2023-09-10 --cycle 00 --product atm
 ```
@@ -27,7 +27,8 @@ For a "typical" period rather than a specific year, use `--climatology` instead
 observed range. `config_fg_interval` (case-8 namelist) sets how often MPAS reads SST.
 
 ## Hand-off to the model
-- `GFS:` → `init_atmosphere` **case 7** (global `*.static.nc`) → `x1.*.init.nc`.
+- `ERA5:` (or `GFS:` for the recent-date alternative) → `init_atmosphere` **case 7**
+  (global `*.static.nc`) → `x1.*.init.nc`.
 - `SST:` → `init_atmosphere` **case 8** → `x1.*.sfc_update.nc`; model `config_sst_update`
   + the `sfc_update` stream. See README *SST / sea-ice update*.
 - **Verify the `init.nc`** (README *Verify the init.nc*) — this catches unit/interpolation
diff --git a/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
index 4a8e52d2..4ccf96bd 100644
--- a/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
+++ b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
@@ -1,9 +1,9 @@
 # Recipe — Hindcast / downscaling, regional (limited-area) mesh
 
 **Goal:** prepare the input *data* for a **downscaling / case study** on a **regional**
-mesh over a **past** period (the most common scientific use). Like the global hindcast
-but with **lateral boundary conditions (LBCs)** from ERA5 through the run. Scope here is
-data preparation only.
+(limited-area) mesh over a **past** period (the most common scientific use). Like the
+global hindcast but with **lateral boundary conditions (LBCs)** from ERA5 through the
+run. Scope here is data preparation only.
 
 | What | Source | When |
 |------|--------|------|
@@ -11,35 +11,85 @@ data preparation only.
 | Lateral boundary data | **ERA5** at the LBC cadence | every `N` h across the run |
 | SST / sea-ice | **OISST observed** daily | at the update cadence |
 
-## 1. Atmosphere initial conditions
+## The `--area` download box
+
+`--area "N W S E"` restricts the **ERA5 download**, not the simulated region — the
+region is the regional mesh (`create_region`, out of scope here). It keeps a long LBC
+series small. Size it to the **regional mesh extent (including its relaxation zone)
+plus a few degrees of margin** on every side; a box that does not cover the mesh makes
+`init_atmosphere` extrapolate or fail at the boundary.
+
+- Order **N W S E**, degrees. Latitude −90..90 (N > S); longitude −180..180 (W/E
+  negative over Brazil).
+- Brazilian equatorial margin (~5°N–5°S, 52°W–34°W) with buffer: `--area "12 -60 -12 -28"`.
+- Global mesh: omit `--area` (download global).
+
+## The run window (month, year, or multiple years)
+
+The same block drives any span — set `START`/`END` and the loops fill it. A single
+month, a single year, or several years are just different endpoints:
+
+| Span | `START` | `END` |
+|------|---------|-------|
+| One month | `2014-09-01` | `2014-09-30` |
+| One year | `2014-01-01` | `2014-12-31` |
+| Multiple years | `2014-01-01` | `2016-12-31` |
+
+ERA5 is taken 6-hourly (the usual LBC cadence). Downloads skip files already present,
+so the loops are resumable and the init time overlapping the first LBC time is free.
+
+```sh
+# --- Window: edit these three lines ------------------------------------------
+START=2014-09-01                                      # run start (init at START 00z)
+END=2014-09-30                                         # last day of the run
+AREA="12 -60 -12 -28"                                 # N W S E, mesh + margin
+
+STOP=$(date -u -d "$END +1 day" +%Y-%m-%d)            # day after END (closing LBC)
+```
+
+### 1. Atmosphere initial conditions (run start)
 ```sh
-./prepare_era5.sh --date 2014-09-10 --time 00          # -> GFS:2014-09-10_00
+./prepare_era5.sh --date "$START" --time 00 --area "$AREA"     # -> ERA5:<START>_00
 ```
 
-## 2. Boundary data through the run (the LBC step consumes these)
-One `GFS:` intermediate per boundary time = ERA5 at the LBC cadence (here 6-hourly):
+### 2. Lateral boundary data, 6-hourly across the window
+One `ERA5:` intermediate per boundary time (00/06/12/18 every day), plus the **closing
+boundary** at `STOP` 00z that the run's stop time needs:
 ```sh
-for d in 2014-09-10 2014-09-11 2014-09-12; do
+d="$START"
+while [ "$d" != "$STOP" ]; do
     for t in 00 06 12 18; do
-        ./prepare_era5.sh --date $d --time $t          # GFS:<d>_<t>
+        ./prepare_era5.sh --date "$d" --time "$t" --area "$AREA"
     done
+    d=$(date -u -d "$d +1 day" +%Y-%m-%d)
 done
+./prepare_era5.sh --date "$STOP" --time 00 --area "$AREA"      # closing boundary
 ```
-(Optionally pass `--area "N W S E"` to ERA5 to download only a box around the regional
-domain and keep the files small.)
+Each day at 6-hourly is 4 boundary times (a small pl+sl GRIB each over the box), so a
+30-day month is `30×4 + 1 = 121` and a year is `365×4 + 1`. A finer cadence
+(3-hourly/hourly) raises boundary fidelity at more cost and storage; it must match
+`config_fg_interval` in the case-9 step (21600 s for 6 h).
 
-## 3. SST / sea-ice over the run (observed, daily)
+### 3. SST / sea-ice over the run (observed, daily)
 ```sh
-./prepare_oisst.sh --start 2014-09-10 --end 2014-09-15
+./prepare_oisst.sh --start "$START" --end "$END"
 ```
 
+For long windows, keep CDS requests **serial** (the queue throttles parallel jobs) and
+budget disk space — a regional box is small per file, but years × 6-hourly accumulates.
+The skip-existing behavior makes the whole window resumable. Download the geog/static
+bundle once (see `../../static_fields/`) and reuse it across all years.
+
 ## Hand-off to the model
-- Initial conditions: `GFS:<init>` → `init_atmosphere` **case 7** (regional `*.static.nc`).
-- LBCs: the `GFS:` series at the boundary times → `init_atmosphere` **case 9** →
-  `lbc.*.nc`; model `config_apply_lbcs=.true.` + the `lbc` stream.
-- SST: `SST:` series → case 8 + `config_sst_update`.
-- **Verify the `init.nc`** (README *Verify the init.nc*).
+- Initial conditions: `ERA5:<START>_00` → `init_atmosphere` **case 7** (regional `*.static.nc`).
+- LBCs: the `ERA5:` series at the boundary times → `init_atmosphere` **case 9** →
+  `lbc.*.nc`; model `config_apply_lbcs=.true.` + the `lbc` stream
+  (`config_fg_interval` = LBC cadence in seconds).
+- SST: `SST:` series → case 8 + `config_sst_update` (see the README *SST / sea-ice
+  update files*).
+- **Verify the `init.nc`** (README *Verify the init.nc*) — the ERA5 path is not yet
+  validated end-to-end here, so this check matters.
 
 > **Out of scope here (separate tutorial / branch):** creating the limited-area mesh
 > (`create_region` from MPAS-Limited-Area) and running `init_atmosphere` case 7/9 and the
-> model. This directory only prepares the `GFS:`/`SST:` intermediates those steps read.
+> model. This directory only prepares the `ERA5:`/`SST:` intermediates those steps read.

From 085753ebd0d4e6a6c69f16d33054e7fd7af379d1 Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Sun, 21 Jun 2026 15:18:38 +0000
Subject: [PATCH 6/7] pre_proc/real_data: fix ERA5 converter dropping PMSL and
 SEAICE

era5_to_intermediate.py silently emitted only 12 of 14 fields because two
single-level lookups did not match the CDS GRIB:
- mean sea level pressure (msl) is at typeOfLevel 'surface' here, not 'meanSea';
- sea-ice cover has GRIB shortName 'ci', which cfgrib loads as variable 'siconc',
  so the 'siconc'-keyed lookup found nothing.

Make the lookups robust: _open_var accepts shortName aliases and returns the
single filtered variable regardless of how cfgrib names it; sea ice tries
'ci'/'siconc'; msl tries 'surface' then 'meanSea'. The converter now writes all
14 fields (adds PMSL and SEAICE). README field-mapping notes the ci/siconc alias.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 usp-utils/pre_proc/real_data/README.md        |  2 +-
 .../real_data/era5_to_intermediate.py         | 59 ++++++++++++-------
 2 files changed, 40 insertions(+), 21 deletions(-)

diff --git a/usp-utils/pre_proc/real_data/README.md b/usp-utils/pre_proc/real_data/README.md
index bca71f83..280cd7e9 100644
--- a/usp-utils/pre_proc/real_data/README.md
+++ b/usp-utils/pre_proc/real_data/README.md
@@ -183,7 +183,7 @@ Both converters emit the same intermediate field names (WPS GFS Vtable):
 | terrain height      | `SOILHGT` | 2D            | HGT sfc   | z(sfc) / g0 |
 | skin temperature    | `SKINTEMP`| 2D            | TMP sfc   | skt |
 | land-sea mask       | `LANDSEA` | 2D            | LAND      | lsm |
-| sea-ice fraction    | `SEAICE`  | 2D            | ICEC      | siconc |
+| sea-ice fraction    | `SEAICE`  | 2D            | ICEC      | ci (a.k.a. siconc) |
 | snow water equiv.   | `SNOW`    | 2D            | WEASD     | sd × 1000 |
 | soil temperature    | `ST<tag>` | 2D (4 layers) | TSOIL     | stl1–4 |
 | soil moisture       | `SM<tag>` | 2D (4 layers) | SOILW     | swvl1–4 |
diff --git a/usp-utils/pre_proc/real_data/era5_to_intermediate.py b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
index 08d05367..43b0bff8 100755
--- a/usp-utils/pre_proc/real_data/era5_to_intermediate.py
+++ b/usp-utils/pre_proc/real_data/era5_to_intermediate.py
@@ -44,14 +44,16 @@
     "v": ("VV", None),
     "r": ("RH", None),
 }
-# 2D surface (typeOfLevel='surface'): shortName -> (intermediate name, transform)
+# 2D surface (typeOfLevel='surface'): intermediate name -> (shortName aliases, transform).
+# Aliases cover ERA5 GRIB naming differences across CDS deliveries (e.g. sea-ice
+# cover comes as 'ci' or 'siconc').
 SURFACE = {
-    "sp": ("PSFC", None),
-    "skt": ("SKINTEMP", None),
-    "z": ("SOILHGT", lambda a: a / G0),
-    "lsm": ("LANDSEA", None),
-    "siconc": ("SEAICE", None),
-    "sd": ("SNOW", lambda a: a * 1000.0),  # m w.e. -> kg m-2
+    "PSFC": (["sp"], None),
+    "SKINTEMP": (["skt"], None),
+    "SOILHGT": (["z"], lambda a: a / G0),
+    "LANDSEA": (["lsm"], None),
+    "SEAICE": (["ci", "siconc"], None),
+    "SNOW": (["sd"], lambda a: a * 1000.0),  # m w.e. -> kg m-2
 }
 # Soil layers: ERA5 shortNames and the matching WPS depth tags (top-bottom, cm)
 SOIL_T = ["stl1", "stl2", "stl3", "stl4"]
@@ -66,14 +68,27 @@ def _open(path, **keys):
     )
 
 
-def _open_var(path, sn, **keys):
-    """Open a single variable by shortName; return the DataArray or None."""
-    try:
-        ds = _open(path, shortName=sn, **keys)
-    except Exception as e:  # noqa: BLE001
-        print(f"[WARN] could not read {sn}: {e}")
-        return None
-    return ds[sn] if sn in ds else None
+def _open_var(path, shortnames, **keys):
+    """Open the first available variable among shortnames; return DataArray or None.
+
+    shortnames may be a single shortName or a list of aliases — ERA5 GRIB names
+    vary across CDS deliveries (e.g. sea-ice cover is 'ci' here but 'siconc'
+    elsewhere), so try each and return the first that opens and is present.
+    """
+    for sn in ([shortnames] if isinstance(shortnames, str) else shortnames):
+        try:
+            ds = _open(path, shortName=sn, **keys)
+        except Exception:  # noqa: BLE001
+            continue
+        if sn in ds:
+            return ds[sn]
+        # cfgrib may name the variable by its CF name, not the GRIB shortName
+        # (e.g. shortName 'ci' -> variable 'siconc'); the shortName+typeOfLevel
+        # filter isolates a single field, so return it.
+        dvars = list(ds.data_vars)
+        if len(dvars) == 1:
+            return ds[dvars[0]]
+    return None
 
 
 def main() -> int:
@@ -129,18 +144,22 @@ def arr3d(a):
 
     # --- 2D surface (open each var individually to avoid cfgrib hypercube clashes) ---
     print(f"[INFO] Reading {sl} (single levels)")
-    for sn, (name, fn) in SURFACE.items():
-        da = _open_var(sl, sn, typeOfLevel="surface")
+    for name, (snames, fn) in SURFACE.items():
+        da = _open_var(sl, snames, typeOfLevel="surface")
         if da is not None:
             a = arr2d(da.values)
             variables.append(pw.V2d(name, fn(a) if fn else a))
         else:
-            print(f"[WARN] surface field missing: {sn}")
-    msl = _open_var(sl, "msl", typeOfLevel="meanSea")
+            print(f"[WARN] surface field missing: {name} ({'/'.join(snames)})")
+    # ERA5 mean sea level pressure: typeOfLevel is 'surface' in the CDS GRIB, but
+    # some deliveries label it 'meanSea' — try both.
+    msl = _open_var(sl, "msl", typeOfLevel="surface")
+    if msl is None:
+        msl = _open_var(sl, "msl", typeOfLevel="meanSea")
     if msl is not None:
         variables.append(pw.V2d("PMSL", arr2d(msl.values)))
     else:
-        print("[WARN] surface field missing: msl")
+        print("[WARN] surface field missing: msl (PMSL)")
 
     # --- Soil layers (stack the 4 ERA5 layers; pywinter Vsl -> ST/SM<tag>) ---
     st = [_open_var(sl, s, typeOfLevel="depthBelowLandLayer") for s in SOIL_T]

From 135a043ee3b57017e1d85bd2ac5e178d599d2ac1 Mon Sep 17 00:00:00 2001
From: daniloceano <danilo.oceano@usp.br>
Date: Sun, 21 Jun 2026 16:34:24 +0000
Subject: [PATCH 7/7] pre_proc/real_data: add resumable ERA5 series downloader
 for long runs
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

prepare_era5_series.sh builds the ERA5 (ERA5:) intermediate series over a date
window — the LBC series of a month / year / multiple years — by wrapping
prepare_era5.sh in a window loop built for unattended runs:

- Resumable: skips every time whose ERA5:<date>_<hh> intermediate already exists,
  so re-running only fetches the gaps.
- Fault-tolerant: a failed CDS request logs FAIL and the series continues.
- nohup-friendly: timestamped log lines; redirect and detach.
- Optional parallelism (--jobs N via xargs -P), with the CDS throttle documented
  (keep N small).
- Timestamps are parsed as UTC and iterated by epoch seconds, so the dates do not
  shift with the server's local timezone.

hindcast_regional.md now recommends it (nohup example) for windows longer than a
few days, keeping the equivalent manual loop for reference; README lists it.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 usp-utils/pre_proc/real_data/README.md        |   1 +
 .../pre_proc/real_data/prepare_era5_series.sh | 108 ++++++++++++++++++
 .../real_data/recipes/hindcast_regional.md    |  43 ++++---
 3 files changed, 137 insertions(+), 15 deletions(-)
 create mode 100755 usp-utils/pre_proc/real_data/prepare_era5_series.sh

diff --git a/usp-utils/pre_proc/real_data/README.md b/usp-utils/pre_proc/real_data/README.md
index 280cd7e9..9197a7d6 100644
--- a/usp-utils/pre_proc/real_data/README.md
+++ b/usp-utils/pre_proc/real_data/README.md
@@ -80,6 +80,7 @@ One-shot wrappers (download + convert in one call, inside the conda env):
 |---------|--------|
 | `prepare_gfs.sh` | GFS atmosphere and/or SST — `--product atm` (default), `sst`, or `both` (one download, both products). |
 | `prepare_era5.sh` | ERA5 atmosphere (download pl+sl → `ERA5:` intermediate). |
+| `prepare_era5_series.sh` | ERA5 `ERA5:` series over a date window (the LBC series of a month/year/years). Resumable, fault-tolerant, `nohup`-friendly, optional `--jobs`. |
 | `prepare_oisst.sh` | OISST SST over a date range — observed (default) or `--climatology` (→ daily `SST:` intermediates). |
 
 ## Dependencies
diff --git a/usp-utils/pre_proc/real_data/prepare_era5_series.sh b/usp-utils/pre_proc/real_data/prepare_era5_series.sh
new file mode 100755
index 00000000..1bb18a0d
--- /dev/null
+++ b/usp-utils/pre_proc/real_data/prepare_era5_series.sh
@@ -0,0 +1,108 @@
+#!/usr/bin/env bash
+#
+# prepare_era5_series.sh — build a *series* of ERA5 WPS intermediates over a date
+# window, for long downscaling runs (the LBC series of a month / year / years).
+# Wraps prepare_era5.sh in a window loop. Designed for unattended runs:
+#
+#   - Resumable: skips every time whose ERA5:<date>_<hh> intermediate already
+#     exists, so re-running only fetches what is missing.
+#   - Fault-tolerant: a failed CDS request logs FAIL and the series continues;
+#     re-run later to pick up the gaps.
+#   - nohup-friendly: every line is timestamped; redirect to a log and detach.
+#   - Optional parallelism (--jobs N). NOTE: the CDS queue throttles concurrent
+#     requests per user, so keep N small (2-3); large N just queues or is rejected.
+#
+# Usage:
+#   nohup ./prepare_era5_series.sh --start 2020-01-01 --end 2020-02-01 \
+#         --area "10 -57 0 -30" --cadence 6 --jobs 2 > era5_2020.log 2>&1 &
+#   tail -f era5_2020.log
+#
+#   --start YYYY-MM-DD   first day (series starts at START 00z)            [required]
+#   --end   YYYY-MM-DD   last day to include, at 00z (set to the run stop) [required]
+#   --area  "N W S E"    ERA5 download box (mesh + margin; W/E negative over Brazil)
+#   --cadence H          LBC cadence in hours (default 6 -> 00/06/12/18)
+#   --jobs N             concurrent downloads (default 1; keep <=3 for CDS)
+#   --env ENV            conda env (default cgfd-usp-mpas)
+#   --outbase DIR        base output dir (default <repo>/met_data/era5)
+#
+# Requires a CDS account + ~/.cdsapirc. SST is a single separate call
+# (prepare_oisst.sh --start --end), not part of this series.
+
+set -uo pipefail   # NOT -e: one failed CDS request must not abort the whole series
+SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+
+START="" END="" AREA="" CAD=6 JOBS=1 ENV="cgfd-usp-mpas" OUTBASE=""
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --start)   START="$2"; shift 2 ;;
+        --end)     END="$2"; shift 2 ;;
+        --area)    AREA="$2"; shift 2 ;;
+        --cadence) CAD="$2"; shift 2 ;;
+        --jobs)    JOBS="$2"; shift 2 ;;
+        --env)     ENV="$2"; shift 2 ;;
+        --outbase) OUTBASE="$2"; shift 2 ;;
+        -h|--help) sed -n '2,30p' "$0" | sed 's/^# \?//'; exit 0 ;;
+        *) echo "Unknown argument: $1" >&2; exit 2 ;;
+    esac
+done
+[ -n "$START" ] && [ -n "$END" ] || { echo "ERROR: --start and --end (YYYY-MM-DD) are required" >&2; exit 2; }
+
+REPO_ROOT=$( cd "$SCRIPT_DIR/../../.." && pwd )
+MET="${OUTBASE:-$REPO_ROOT/met_data/era5}"
+
+ts() { date -u +%Y-%m-%dT%H:%M:%SZ; }
+
+# Download+convert one (date, hour), skipping if the intermediate already exists.
+do_one() {
+    local d="$1" hh="$2"
+    local ymd="${d//-/}"
+    local intf="$MET/${ymd}${hh}/ERA5:${d}_${hh}"
+    if [ -s "$intf" ]; then
+        echo "[$(ts)] SKIP $d ${hh}z (have ERA5:${d}_${hh})"
+        return 0
+    fi
+    echo "[$(ts)] GET  $d ${hh}z"
+    local args=(--date "$d" --time "$hh" --env "$ENV")
+    [ -n "$AREA" ] && args+=(--area "$AREA")
+    if "$SCRIPT_DIR/prepare_era5.sh" "${args[@]}"; then
+        echo "[$(ts)] DONE $d ${hh}z"
+    else
+        echo "[$(ts)] FAIL $d ${hh}z (re-run to retry)"
+        return 1
+    fi
+}
+export -f do_one ts
+export SCRIPT_DIR MET AREA ENV
+
+# Build the task list: timestamps from START 00z to END 00z (inclusive), step CAD h.
+# Parse as UTC (the '... UTC' suffix) and iterate by epoch seconds so the result is
+# independent of the server's local timezone and of date(1)'s output locale.
+start_ts=$( date -u -d "$START 00:00 UTC" +%s ) || { echo "ERROR: bad --start date" >&2; exit 2; }
+end_ts=$(   date -u -d "$END 00:00 UTC"   +%s ) || { echo "ERROR: bad --end date" >&2; exit 2; }
+[ "$start_ts" -le "$end_ts" ] || { echo "ERROR: --start is after --end" >&2; exit 2; }
+step=$(( CAD * 3600 ))
+tasks=()
+t_ts="$start_ts"
+while [ "$t_ts" -le "$end_ts" ]; do
+    tasks+=("$( date -u -d "@$t_ts" +%Y-%m-%d ) $( date -u -d "@$t_ts" +%H )")
+    t_ts=$(( t_ts + step ))
+done
+
+echo "[$(ts)] ERA5 series: ${#tasks[@]} times | ${START}..${END} | cadence ${CAD}h | jobs ${JOBS} | area '${AREA:-global}'"
+
+if [ "$JOBS" -gt 1 ]; then
+    printf '%s\n' "${tasks[@]}" | xargs -P "$JOBS" -n 2 bash -c 'do_one "$1" "$2"' _
+else
+    for task in "${tasks[@]}"; do do_one $task; done
+fi
+
+# Summary: how many intermediates are present now (resumability check).
+have=0
+for task in "${tasks[@]}"; do
+    set -- $task; d="$1"; hh="$2"; ymd="${d//-/}"
+    [ -s "$MET/${ymd}${hh}/ERA5:${d}_${hh}" ] && have=$((have + 1))
+done
+miss=$(( ${#tasks[@]} - have ))
+echo "[$(ts)] summary: ${have}/${#tasks[@]} intermediates present, ${miss} missing"
+[ "$miss" -gt 0 ] && echo "[$(ts)] re-run the same command to retry the ${miss} missing time(s)."
+exit 0
diff --git a/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
index 4ccf96bd..dcedd3cf 100644
--- a/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
+++ b/usp-utils/pre_proc/real_data/recipes/hindcast_regional.md
@@ -52,33 +52,46 @@ STOP=$(date -u -d "$END +1 day" +%Y-%m-%d)            # day after END (closing L
 ./prepare_era5.sh --date "$START" --time 00 --area "$AREA"     # -> ERA5:<START>_00
 ```
 
-### 2. Lateral boundary data, 6-hourly across the window
-One `ERA5:` intermediate per boundary time (00/06/12/18 every day), plus the **closing
-boundary** at `STOP` 00z that the run's stop time needs:
+### 2. Lateral boundary data across the window (6-hourly)
+One `ERA5:` intermediate per boundary time. For anything longer than a few days use
+the **series helper**, which is built for unattended runs — it skips times already
+fetched (resumable), keeps going if a single CDS request fails, and is meant to run
+under `nohup`:
 ```sh
-d="$START"
-while [ "$d" != "$STOP" ]; do
-    for t in 00 06 12 18; do
-        ./prepare_era5.sh --date "$d" --time "$t" --area "$AREA"
-    done
-    d=$(date -u -d "$d +1 day" +%Y-%m-%d)
-done
-./prepare_era5.sh --date "$STOP" --time 00 --area "$AREA"      # closing boundary
+nohup ./prepare_era5_series.sh --start "$START" --end "$STOP" \
+      --area "$AREA" --cadence 6 --jobs 2 > era5_series.log 2>&1 &
+tail -f era5_series.log
 ```
+`--end "$STOP"` makes the series include the **closing boundary** at the next-day 00z
+that the run's stop time needs. `--jobs` adds parallelism, but the CDS queue throttles
+concurrent requests, so keep it small (2-3). Re-run the same command to retry any
+times that failed.
+
 Each day at 6-hourly is 4 boundary times (a small pl+sl GRIB each over the box), so a
 30-day month is `30×4 + 1 = 121` and a year is `365×4 + 1`. A finer cadence
 (3-hourly/hourly) raises boundary fidelity at more cost and storage; it must match
 `config_fg_interval` in the case-9 step (21600 s for 6 h).
 
+> Under the hood it just loops `prepare_era5.sh` over the window; by hand that is:
+> ```sh
+> d="$START"
+> while [ "$d" != "$STOP" ]; do
+>     for t in 00 06 12 18; do ./prepare_era5.sh --date "$d" --time "$t" --area "$AREA"; done
+>     d=$(date -u -d "$d +1 day" +%Y-%m-%d)
+> done
+> ./prepare_era5.sh --date "$STOP" --time 00 --area "$AREA"   # closing boundary
+> ```
+
 ### 3. SST / sea-ice over the run (observed, daily)
 ```sh
 ./prepare_oisst.sh --start "$START" --end "$END"
 ```
 
-For long windows, keep CDS requests **serial** (the queue throttles parallel jobs) and
-budget disk space — a regional box is small per file, but years × 6-hourly accumulates.
-The skip-existing behavior makes the whole window resumable. Download the geog/static
-bundle once (see `../../static_fields/`) and reuse it across all years.
+For long windows, keep `--jobs` low (2-3) — the CDS queue throttles concurrent
+requests, so more workers just queue or get rejected — and budget disk space (a
+regional box is small per file, but years × 6-hourly accumulates). The series helper
+is resumable, so an interrupted multi-year pull picks up where it stopped. Download
+the geog/static bundle once (see `../../static_fields/`) and reuse it across all years.
 
 ## Hand-off to the model
 - Initial conditions: `ERA5:<START>_00` → `init_atmosphere` **case 7** (regional `*.static.nc`).