prop/rust/prop_grid_rs/src/pipeline.rs

//! End-to-end f01..f18 chain step: fetch surface + pressure GRIBs, decode,
//! build per-cell `Conditions`, score every band, write score-grid files
//! (`<band>/<iso>.prop`) to the shared scores directory.
//!
//! This is the Rust equivalent of `PropagationGridWorker.process_forecast_hour/4`
//! minus the f00-only enrichment (native duct, NEXRAD, commercial link
//! boost, ProfilesFile write) — Elixir retains f00.

use std::path::Path;
use std::sync::Arc;

use chrono::{DateTime, Datelike, Timelike, Utc};

use crate::band_config;
use crate::commercial::{self, LinkLookupEntry};
use crate::decoder::{self, CellValues, PointGrid};
use crate::fetcher::{self, HrrrClient, Product};
use crate::grid::wgrib2_grid_spec;
use crate::native_duct;
use crate::nexrad::{self, NexradObservation};
use crate::profiles_file::{self, CellEntry};
use crate::scorer::{self, Conditions};
use crate::scores_file::{self, ScorePoint};
use crate::sounding_params::{self, Level};
use crate::weather_scalar_file::{self, ScalarRow};

#[derive(Debug, thiserror::Error)]
pub enum PipelineError {
    #[error("fetch: {0}")]
    Fetch(#[from] fetcher::FetchError),
    #[error("decode: {0}")]
    Decode(#[from] decoder::DecodeError),
    #[error("write: {0}")]
    Write(#[from] scores_file::WriteError),
    #[error("profile write: {0}")]
    ProfileWrite(#[from] profiles_file::WriteError),
    #[error("scalar write: {0}")]
    ScalarWrite(#[from] weather_scalar_file::WriteError),
    #[error("native duct: {0}")]
    NativeDuct(#[from] native_duct::NativeDuctError),
    #[error("nexrad: {0}")]
    Nexrad(#[from] nexrad::NexradError),
    #[error("commercial: {0}")]
    Commercial(#[from] commercial::CommercialError),
    #[error("forecast hour 0 is reserved for Elixir")]
    F00Reserved,
    #[error("surface grib was empty")]
    EmptyGrib,
}

#[derive(Debug, Clone)]
pub struct ChainStepInput {
    pub run_time: DateTime<Utc>,
    pub forecast_hour: u8,
}

impl ChainStepInput {
    pub fn valid_time(&self) -> DateTime<Utc> {
        self.run_time + chrono::Duration::hours(self.forecast_hour as i64)
    }
}

#[derive(Debug, Clone)]
pub struct ChainStepStats {
    pub score_files_written: u32,
    pub point_count: u32,
    pub band_count: u32,
    /// Number of cells written to the per-valid_time profile file.
    /// Always `point_count` after the f01..f18 profile-file write was
    /// added on 2026-04-25 — kept as a separate field so the log line
    /// is symmetric with the analysis step's `profile_cells_written`
    /// in `RunStepStats`, and so a future failure mode that writes
    /// scores but skips the profile file shows up as a divergence.
    pub profile_cells_written: u32,
}

#[tracing::instrument(
    name = "pipeline.run_chain_step",
    skip(client, scores_dir),
    fields(
        run_time = %step.run_time,
        forecast_hour = step.forecast_hour,
    )
)]
pub async fn run_chain_step(
    client: &HrrrClient,
    scores_dir: &Path,
    step: &ChainStepInput,
) -> Result<ChainStepStats, PipelineError> {
    if step.forecast_hour == 0 {
        return Err(PipelineError::F00Reserved);
    }
    let valid_time = step.valid_time();
    let date = step.run_time.date_naive();
    let hour = step.run_time.hour() as u8;

    let sfc_wanted = fetcher::surface_messages_owned();
    let prs_wanted = fetcher::pressure_messages_grid();
    let sfc_fut = client.fetch_product_blob(
        date,
        hour,
        Product::Surface,
        step.forecast_hour,
        &sfc_wanted,
    );
    let prs_fut = client.fetch_product_blob(
        date,
        hour,
        Product::Pressure,
        step.forecast_hour,
        &prs_wanted,
    );
    let (sfc_blob, prs_blob) = tokio::try_join!(sfc_fut, prs_fut)?;
    if sfc_blob.is_empty() {
        return Err(PipelineError::EmptyGrib);
    }

    let grid_spec = wgrib2_grid_spec();
    let sfc_pattern = match_pattern(&sfc_wanted);
    let prs_pattern = match_pattern(&prs_wanted);

    // Decoder runs wgrib2 subprocess — block so we don't hog the tokio
    // reactor. `spawn_blocking` lifts it onto the dedicated pool.
    let grid_spec_sfc = grid_spec;
    let sfc_grid = tokio::task::spawn_blocking(move || {
        decoder::extract_grid(&sfc_blob, &sfc_pattern, grid_spec_sfc)
    })
    .await
    .expect("blocking join")?;

    let grid_spec_prs = grid_spec;
    let prs_grid = tokio::task::spawn_blocking(move || {
        decoder::extract_grid(&prs_blob, &prs_pattern, grid_spec_prs)
    })
    .await
    .expect("blocking join")?;

    let merged = merge_grids(sfc_grid, prs_grid);
    let point_count = merged.len() as u32;

    // Build `prepared` (for scoring) and `profile_entries` (for the
    // per-valid_time profile file consumed by /skewt and the LiveView
    // point-detail panels) in a single pass over `merged`. We need a
    // borrowing iterator (not `into_iter`) because `cell_to_profile_entry`
    // also reads from the cell — `merged` is dropped immediately after
    // this loop finishes, so the ~200 MB string-keyed HashMaps don't
    // outlive the scoring stage. The price is briefly retaining both
    // vecs at the same time as `merged`, which is well inside the
    // 512 MiB cgroup ceiling on talos5.
    let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> =
        Vec::with_capacity(merged.len());
    let mut profile_entries: Vec<CellEntry> = Vec::with_capacity(merged.len());
    let mut scalar_rows: Vec<ScalarRow> = Vec::with_capacity(merged.len());
    for (key, cell) in merged.iter() {
        let (lat, lon) = decoder::key_to_latlon(*key);
        if let Some(conditions) = cell_to_conditions(cell, lat, lon, &valid_time) {
            let invariants = scorer::precompute_band_invariants(&conditions);
            prepared.push((lat, lon, conditions, invariants));
            profile_entries.push(cell_to_profile_entry(lat, lon, cell));
            if let Some(row) = weather_scalar_file::derive_row(lat, lon, valid_time, cell) {
                scalar_rows.push(row);
            }
        }
    }
    drop(merged);
    let profile_cells_written = profile_entries.len() as u32;

    // Spawn the profile-file write on the blocking pool so it overlaps
    // with the scoring loop instead of serialising on the hot path.
    let scores_dir_owned = scores_dir.to_path_buf();
    let scores_dir_for_profiles = scores_dir_owned.clone();
    let profile_future = {
        let profiles = profile_entries;
        tokio::task::spawn_blocking(move || {
            profiles_file::write_atomic(&scores_dir_for_profiles, valid_time, &profiles)
                .map(|_| 0u32)
        })
    };

    // Persist the derived scalar artifact alongside the raw profile so
    // Elixir's `/weather` reads land on a pre-derived chunk file
    // instead of falling back to the cold ProfilesFile decode + per-cell
    // SoundingParams + WeatherLayers derivation. Identical wire format
    // (gzipped MessagePack chunked by 5°×5°) to what the Elixir writer
    // produces — the reader doesn't care which side wrote the bytes.
    let scores_dir_for_scalars = scores_dir_owned.clone();
    let scalar_future = {
        let rows = scalar_rows;
        tokio::task::spawn_blocking(move || {
            weather_scalar_file::write_atomic(&scores_dir_for_scalars, valid_time, &rows)
                .map(|_| 0u32)
        })
    };

    // Step 2: score and write every band in parallel.
    //
    // Scoring: rayon's thread pool over the 23 bands. Each band is
    // independent pure math across the ~92k-cell `prepared` slice;
    // sharing is read-only so no locking. This cuts per-step scoring
    // wall from ~20s serial to ~5-8s on a 4-core box.
    //
    // Writing: each band file is an independent NFS fsync. Launching
    // them as parallel `spawn_blocking` tasks lets them overlap on the
    // blocking pool instead of serializing one-at-a-time on the hot
    // path. Memory stays bounded — each Vec<ScorePoint> is ~1.8 MB so
    // 23 concurrent ≈ ~40 MB peak.
    let bands = band_config::all_bands();
    let prepared = Arc::new(prepared);
    let scored: Vec<(u32, Vec<ScorePoint>)> = {
        use rayon::prelude::*;
        bands
            .par_iter()
            .map(|band| {
                let mut scores: Vec<ScorePoint> = Vec::with_capacity(prepared.len());
                for (lat, lon, conditions, invariants) in prepared.iter() {
                    let r = scorer::composite_score_with(conditions, band, Some(*invariants));
                    scores.push(ScorePoint {
                        lat: *lat,
                        lon: *lon,
                        score: r.score,
                    });
                }
                (band.freq_mhz, scores)
            })
            .collect()
    };

    let write_handles: Vec<_> = scored
        .into_iter()
        .map(|(band_mhz, scores)| {
            let dir = scores_dir_owned.clone();
            tokio::task::spawn_blocking(move || {
                scores_file::write_atomic(&dir, band_mhz, valid_time, &scores)
            })
        })
        .collect();
    for h in write_handles {
        h.await.expect("blocking join")?;
    }
    profile_future
        .await
        .expect("blocking join")
        .map_err(PipelineError::ProfileWrite)?;
    scalar_future
        .await
        .expect("blocking join")
        .map_err(PipelineError::ScalarWrite)?;
    let files_written = bands.len() as u32;

    Ok(ChainStepStats {
        score_files_written: files_written,
        point_count,
        band_count: bands.len() as u32,
        profile_cells_written,
    })
}

/// wgrib2 `-match` regex: `":(A|B|C):"`. Var names only — levels are
/// post-filtered when we read the binary back.
fn match_pattern(wanted: &[(String, String)]) -> String {
    let mut vars: Vec<&str> = wanted.iter().map(|(v, _)| v.as_str()).collect();
    vars.sort();
    vars.dedup();
    format!(":({}):", vars.join("|"))
}

fn merge_grids(mut sfc: PointGrid, prs: PointGrid) -> PointGrid {
    for (k, v) in prs {
        sfc.entry(k).or_default().extend(v);
    }
    sfc
}

fn cell_to_conditions(
    cell: &CellValues,
    lat: f64,
    lon: f64,
    valid_time: &DateTime<Utc>,
) -> Option<Conditions> {
    let tmp_k = cell.get("TMP:2 m above ground").copied()?;
    let dpt_k = cell.get("DPT:2 m above ground").copied()?;
    let temp_c = (tmp_k - 273.15) as f64;
    let dewpoint_c = (dpt_k - 273.15) as f64;
    let temp_f = scorer::c_to_f(temp_c);
    let dewpoint_f = scorer::c_to_f(dewpoint_c);
    let abs_humidity = scorer::absolute_humidity(temp_c, dewpoint_c);

    let wind_speed_kts = match (
        cell.get("UGRD:10 m above ground").copied(),
        cell.get("VGRD:10 m above ground").copied(),
    ) {
        (Some(u), Some(v)) => Some(scorer::wind_speed_kts(u as f64, v as f64)),
        _ => None,
    };

    let sky_cover_pct = cell
        .get("TCDC:entire atmosphere")
        .copied()
        .map(|v| v as f64);
    let pwat_mm = cell
        .get("PWAT:entire atmosphere (considered as a single layer)")
        .copied()
        .map(|v| v as f64);
    let pressure_mb = cell
        .get("PRES:surface")
        .copied()
        .map(|pa| pa as f64 / 100.0);
    // Mirror Elixir's merged_rain_rate: pick the heavier of HRRR's
    // accumulation-derived rate and NEXRAD's reflectivity-derived rate
    // so a fast convective cell that hasn't yet shown up in the hourly
    // APCP still triggers the rain penalty.
    let hrrr_rate = cell
        .get("APCP:surface")
        .copied()
        .map(|mm| mm as f64)
        .unwrap_or(0.0);
    let nexrad_rate = scorer::dbz_to_rain_rate_mmhr(
        cell.get("nexrad_max_reflectivity_dbz")
            .copied()
            .map(|v| v as f64),
    );
    let merged_rate = hrrr_rate.max(nexrad_rate);
    let rain_rate_mmhr = if merged_rate > 0.0 {
        Some(merged_rate)
    } else {
        None
    };
    let bl_depth_m = cell.get("HPBL:surface").copied().map(|v| v as f64);
    let best_duct_band_ghz = cell.get("best_duct_freq_ghz").copied().map(|v| v as f64);

    let levels: Vec<Level> = fetcher::grid_level_keys()
        .iter()
        .filter_map(|k| {
            let t = cell.get(k.tmp.as_str()).copied()?;
            let d = cell.get(k.dpt.as_str()).copied();
            let h = cell.get(k.hgt.as_str()).copied()?;
            Some(Level {
                pres_mb: k.pres_mb,
                hght_m: h as f64,
                tmpc: (t - 273.15) as f64,
                dwpc: d.map(|v| (v - 273.15) as f64),
            })
        })
        .collect();
    let min_refractivity_gradient = sounding_params::min_refractivity_gradient(levels);

    Some(Conditions {
        abs_humidity,
        temp_f,
        dewpoint_f,
        wind_speed_kts,
        sky_cover_pct,
        utc_hour: valid_time.hour() as u8,
        utc_minute: valid_time.minute() as u8,
        month: valid_time.month() as u8,
        longitude: lon,
        latitude: Some(lat),
        pressure_mb,
        prev_pressure_mb: None,
        rain_rate_mmhr,
        min_refractivity_gradient,
        bl_depth_m,
        pwat_mm,
        best_duct_band_ghz,
        bulk_richardson: None,
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn match_pattern_dedups_vars() {
        let wanted = vec![
            ("TMP".into(), "2 m above ground".into()),
            ("DPT".into(), "2 m above ground".into()),
            ("TMP".into(), "1000 mb".into()),
        ];
        assert_eq!(match_pattern(&wanted), ":(DPT|TMP):");
    }

    #[test]
    fn cell_to_conditions_maps_fields() {
        use chrono::TimeZone;
        let mut cell = CellValues::new();
        cell.insert("TMP:2 m above ground".into(), 298.15); // 25 °C
        cell.insert("DPT:2 m above ground".into(), 293.15); // 20 °C
        cell.insert("UGRD:10 m above ground".into(), 3.0);
        cell.insert("VGRD:10 m above ground".into(), 4.0);
        cell.insert("TCDC:entire atmosphere".into(), 30.0);
        cell.insert("PRES:surface".into(), 101_000.0);
        cell.insert("HPBL:surface".into(), 400.0);
        cell.insert(
            "PWAT:entire atmosphere (considered as a single layer)".into(),
            25.0,
        );

        let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
        let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
        assert!((c.temp_f - 77.0).abs() < 0.1);
        assert!((c.dewpoint_f - 68.0).abs() < 0.1);
        assert!((c.wind_speed_kts.unwrap() - 9.72).abs() < 0.1);
        assert!((c.pressure_mb.unwrap() - 1010.0).abs() < 0.01);
        assert_eq!(c.sky_cover_pct, Some(30.0));
        assert_eq!(c.bl_depth_m, Some(400.0));
        assert_eq!(c.pwat_mm, Some(25.0));
        assert_eq!(c.month, 6);
        assert_eq!(c.longitude, -97.0);
    }

    #[test]
    fn cell_missing_surface_returns_none() {
        use chrono::TimeZone;
        let cell = CellValues::new();
        let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
        assert!(cell_to_conditions(&cell, 32.0, -97.0, &vt).is_none());
    }

    #[test]
    fn nexrad_reflectivity_lifts_rain_rate_above_hrrr_precip() {
        // Mirrors the Elixir scorer's rain merge: when HRRR's hourly
        // accumulation hasn't caught a fast convective cell but NEXRAD
        // sees 45 dBZ overhead, take whichever rate is higher. Without
        // this the Rust pipeline silently underestimates rain attenuation.
        use chrono::TimeZone;
        let mut cell = CellValues::new();
        cell.insert("TMP:2 m above ground".into(), 295.0);
        cell.insert("DPT:2 m above ground".into(), 285.0);
        cell.insert("APCP:surface".into(), 0.0);
        cell.insert("nexrad_max_reflectivity_dbz".into(), 45.0);

        let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
        let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
        let rate = c.rain_rate_mmhr.expect("nexrad-derived rain rate present");
        assert!(rate > 5.0, "expected >5 mm/hr from 45 dBZ, got {rate}");
    }

    #[test]
    fn duct_freq_threads_into_best_duct_band_ghz() {
        // After native_duct::merge_duct_grid runs, every cell with a
        // detected trapping layer carries `best_duct_freq_ghz`. The
        // scorer reads `best_duct_band_ghz` from Conditions for the 15%
        // Native Duct Boost. cell_to_conditions must wire one to the
        // other or the boost never fires from the Rust pipeline.
        use chrono::TimeZone;
        let mut cell = CellValues::new();
        cell.insert("TMP:2 m above ground".into(), 295.0);
        cell.insert("DPT:2 m above ground".into(), 285.0);
        cell.insert("best_duct_freq_ghz".into(), 24.0);

        let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
        let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
        assert_eq!(c.best_duct_band_ghz, Some(24.0));
    }

    #[test]
    fn rejects_forecast_hour_zero() {
        use chrono::TimeZone;
        let client = HrrrClient::default_base().unwrap();
        let dir = tempfile::tempdir().unwrap();
        let step = ChainStepInput {
            run_time: Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(),
            forecast_hour: 0,
        };
        let rt = tokio::runtime::Runtime::new().unwrap();
        let err = rt
            .block_on(run_chain_step(&client, dir.path(), &step))
            .unwrap_err();
        assert!(matches!(err, PipelineError::F00Reserved));
    }

    /// Integration-ish: hand-build a 2-cell surface + pressure grid,
    /// merge → cell_to_conditions → scorer → scores_file::write_atomic
    /// → scores_file::read, asserting the full chain agrees with
    /// itself. Covers the post-fetch hot path (the one that scales
    /// with grid size) without standing up wgrib2.
    #[test]
    fn merge_score_write_read_roundtrip() {
        use chrono::TimeZone;

        // decoder encodes lat/lon at millidegree precision in the key.
        fn latlon_to_key(lat: f64, lon: f64) -> (i32, i32) {
            ((lat * 1000.0).round() as i32, (lon * 1000.0).round() as i32)
        }

        let mut sfc = PointGrid::new();
        let mut prs = PointGrid::new();

        let cells = [(32.0_f64, -97.0_f64), (32.03_f64, -97.03_f64)];
        for &(lat, lon) in &cells {
            let mut s = CellValues::new();
            s.insert("TMP:2 m above ground".into(), 295.0);
            s.insert("DPT:2 m above ground".into(), 285.0);
            s.insert("UGRD:10 m above ground".into(), 2.0);
            s.insert("VGRD:10 m above ground".into(), 1.0);
            s.insert("TCDC:entire atmosphere".into(), 20.0);
            s.insert("PRES:surface".into(), 101_000.0);
            s.insert("HPBL:surface".into(), 500.0);
            s.insert(
                "PWAT:entire atmosphere (considered as a single layer)".into(),
                20.0,
            );
            sfc.insert(latlon_to_key(lat, lon), s);

            let mut p = CellValues::new();
            p.insert("TMP:1000 mb".into(), 295.0);
            p.insert("DPT:1000 mb".into(), 285.0);
            p.insert("HGT:1000 mb".into(), 100.0);
            p.insert("TMP:975 mb".into(), 293.0);
            p.insert("DPT:975 mb".into(), 280.0);
            p.insert("HGT:975 mb".into(), 300.0);
            prs.insert(latlon_to_key(lat, lon), p);
        }

        let merged = merge_grids(sfc, prs);
        assert_eq!(merged.len(), 2);

        let valid_time = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();

        let prepared: Vec<_> = merged
            .into_iter()
            .filter_map(|(key, cell)| {
                let (lat, lon) = crate::decoder::key_to_latlon(key);
                let c = cell_to_conditions(&cell, lat, lon, &valid_time)?;
                let inv = scorer::precompute_band_invariants(&c);
                Some((lat, lon, c, inv))
            })
            .collect();
        assert_eq!(prepared.len(), 2);

        let band = &band_config::all_bands()[0];
        let scores: Vec<ScorePoint> = prepared
            .iter()
            .map(|(lat, lon, c, inv)| {
                let r = scorer::composite_score_with(c, band, Some(*inv));
                ScorePoint {
                    lat: *lat,
                    lon: *lon,
                    score: r.score,
                }
            })
            .collect();
        assert_eq!(scores.len(), 2);
        for s in &scores {
            assert!(s.score <= 100);
        }

        // Write to an atomic .prop file and read it back.
        let dir = tempfile::tempdir().unwrap();
        scores_file::write_atomic(dir.path(), band.freq_mhz, valid_time, &scores).unwrap();

        let path = scores_file::path_for(dir.path(), band.freq_mhz, valid_time);
        let bytes = std::fs::read(&path).expect("read back score file");
        let decoded = scores_file::decode(&bytes).expect("decode");
        // Decoded grid covers the full CONUS bounding box as a dense
        // byte array; our 2 score points land at their respective
        // cells and the rest are sentinel no-data. Spot-check: header
        // fields survived the round trip and the body is the expected
        // size for a row_major n_rows × n_cols grid.
        assert_eq!(decoded.band_mhz, band.freq_mhz);
        assert_eq!(
            decoded.body.len(),
            decoded.n_rows as usize * decoded.n_cols as usize
        );
    }
}

// =====================================================================
// Analysis pipeline (f00): fetch surface+pressure + native duct, merge
// with NEXRAD + commercial-link degradation, score all bands, write the
// ProfilesFile (MessagePack) and the band score files.
// =====================================================================

#[derive(Debug, Clone)]
pub struct AnalysisStepInput {
    pub run_time: DateTime<Utc>,
}

impl AnalysisStepInput {
    pub fn valid_time(&self) -> DateTime<Utc> {
        // f00 analysis valid_time == run_time.
        self.run_time
    }
}

#[derive(Debug, Clone)]
pub struct AnalysisStepStats {
    pub score_files_written: u32,
    pub point_count: u32,
    pub band_count: u32,
    pub profile_cells_written: u32,
    pub duct_cells: u32,
    pub nexrad_cells_with_echo: u32,
    pub commercial_cells_boosted: u32,
}

#[tracing::instrument(
    name = "pipeline.run_analysis_step",
    skip(client, pool, scores_dir),
    fields(run_time = %step.run_time)
)]
pub async fn run_analysis_step(
    client: &HrrrClient,
    pool: &sqlx::PgPool,
    scores_dir: &Path,
    step: &AnalysisStepInput,
) -> Result<AnalysisStepStats, PipelineError> {
    let valid_time = step.valid_time();
    let date = step.run_time.date_naive();
    let hour = step.run_time.hour() as u8;
    let grid_spec = wgrib2_grid_spec();

    // Parallel fetch + decode of surface + pressure + native duct. All
    // three are independent HRRR products; overlap the HTTP legs and
    // the subsequent wgrib2 decodes so the long pole is a single
    // decode, not three serial decodes.
    let sfc_wanted = fetcher::surface_messages_owned();
    let prs_wanted = fetcher::pressure_messages_grid();
    let sfc_pattern = match_pattern(&sfc_wanted);
    let prs_pattern = match_pattern(&prs_wanted);

    // All three futures return different error types; lift into
    // PipelineError inside each branch so try_join! stays happy.
    let sfc_fut = async {
        client
            .fetch_product_blob(date, hour, Product::Surface, 0, &sfc_wanted)
            .await
            .map_err(PipelineError::Fetch)
    };
    let prs_fut = async {
        client
            .fetch_product_blob(date, hour, Product::Pressure, 0, &prs_wanted)
            .await
            .map_err(PipelineError::Fetch)
    };
    let duct_fut = async {
        native_duct::fetch_native_duct_grid(client, date, hour, 0, grid_spec)
            .await
            .map_err(PipelineError::NativeDuct)
    };

    let (sfc_blob, prs_blob, duct_map) = tokio::try_join!(sfc_fut, prs_fut, duct_fut)?;

    let duct_count = duct_map.len() as u32;

    let sfc_grid = tokio::task::spawn_blocking(move || {
        decoder::extract_grid(&sfc_blob, &sfc_pattern, grid_spec)
    })
    .await
    .expect("blocking join")?;

    let prs_grid = tokio::task::spawn_blocking(move || {
        decoder::extract_grid(&prs_blob, &prs_pattern, grid_spec)
    })
    .await
    .expect("blocking join")?;

    // Merge pressure into surface, then fold duct metrics into every
    // cell that has both. The duct values live alongside the raw
    // CellValues strings so build_grid_cache_rows can pick them up
    // under their atom keys after the Elixir reader normalises.
    let mut merged = merge_grids(sfc_grid, prs_grid);
    merged = native_duct::merge_duct_grid(merged, &duct_map);

    // NEXRAD composite reflectivity overlay — only valid for f00
    // because forecast hours can't see the future radar image.
    let nexrad_http = reqwest::Client::builder()
        .user_agent("prop-grid-rs/0.1")
        .build()
        .expect("reqwest client");
    let nexrad_points: Vec<(f64, f64)> =
        merged.keys().map(|&k| decoder::key_to_latlon(k)).collect();
    let nexrad_obs: Vec<NexradObservation> = match nexrad::fetch_frame(
        &nexrad_http,
        valid_time,
        &nexrad_points,
    )
    .await
    {
        Ok(obs) => obs,
        Err(e) => {
            tracing::warn!(error = %e, "NEXRAD fetch failed — continuing without radar overlay");
            Vec::new()
        }
    };
    let nexrad_cells_with_echo = apply_nexrad(&mut merged, &nexrad_obs);

    // Commercial-link degradation — precompute once, then fold into
    // every in-range cell. Only ~7 links cluster around DFW so most
    // cells stay untouched.
    let commercial_lookup: Vec<LinkLookupEntry> = commercial::build_link_lookup(pool, valid_time)
        .await
        .map_err(PipelineError::Commercial)?;
    let commercial_cells_boosted = apply_commercial(&mut merged, &commercial_lookup);

    let point_count = merged.len() as u32;

    // Build the ProfilesFile cell list (the raw enriched grid) and
    // the Conditions vector (for scoring) in one pass. Profile keeps
    // atmospheric + duct + NEXRAD + commercial fields under
    // string-keyed msgpack; Rust doesn't need to mutate those further.
    let mut profile_entries: Vec<CellEntry> = Vec::with_capacity(merged.len());
    let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> =
        Vec::with_capacity(merged.len());
    let mut scalar_rows: Vec<ScalarRow> = Vec::with_capacity(merged.len());
    for (key, cell) in merged.iter() {
        let (lat, lon) = decoder::key_to_latlon(*key);
        if let Some(conditions) = cell_to_conditions(cell, lat, lon, &valid_time) {
            let invariants = scorer::precompute_band_invariants(&conditions);
            prepared.push((lat, lon, conditions, invariants));
            profile_entries.push(cell_to_profile_entry(lat, lon, cell));
            if let Some(row) = weather_scalar_file::derive_row(lat, lon, valid_time, cell) {
                scalar_rows.push(row);
            }
        }
    }

    // Write the profile file on the blocking pool — msgpack encode +
    // gzip + atomic rename. One file per run, ~2 MB compressed, so the
    // write runs well inside the scoring loop's wall time.
    let scores_dir_owned = scores_dir.to_path_buf();
    let scores_dir_for_profiles = scores_dir_owned.clone();
    let profile_future = {
        let profiles = profile_entries;
        tokio::task::spawn_blocking(move || {
            profiles_file::write_atomic(&scores_dir_for_profiles, valid_time, &profiles)
                .map(|_| 0u32)
        })
    };

    // Same wire format as Elixir's `ScalarFile`. See run_chain_step
    // (forecast hours) for the full rationale — both code paths go
    // through the same writer so f00 and f01..f18 land identical
    // chunked artifacts on NFS.
    let scores_dir_for_scalars = scores_dir_owned.clone();
    let scalar_future = {
        let rows = scalar_rows;
        tokio::task::spawn_blocking(move || {
            weather_scalar_file::write_atomic(&scores_dir_for_scalars, valid_time, &rows)
                .map(|_| 0u32)
        })
    };

    // Bands and score writes run in parallel, same shape as
    // run_chain_step. Rayon across bands, try_join_all across writes.
    let bands = band_config::all_bands();
    let prepared_arc = Arc::new(prepared);
    let scored: Vec<(u32, Vec<ScorePoint>)> = {
        use rayon::prelude::*;
        bands
            .par_iter()
            .map(|band| {
                let mut scores: Vec<ScorePoint> = Vec::with_capacity(prepared_arc.len());
                for (lat, lon, conditions, invariants) in prepared_arc.iter() {
                    let r = scorer::composite_score_with(conditions, band, Some(*invariants));
                    scores.push(ScorePoint {
                        lat: *lat,
                        lon: *lon,
                        score: r.score,
                    });
                }
                (band.freq_mhz, scores)
            })
            .collect()
    };

    let write_handles: Vec<_> = scored
        .into_iter()
        .map(|(band_mhz, scores)| {
            let dir = scores_dir_owned.clone();
            tokio::task::spawn_blocking(move || {
                scores_file::write_atomic(&dir, band_mhz, valid_time, &scores)
            })
        })
        .collect();
    for h in write_handles {
        h.await.expect("blocking join")?;
    }
    profile_future
        .await
        .expect("blocking join")
        .map_err(PipelineError::ProfileWrite)?;
    scalar_future
        .await
        .expect("blocking join")
        .map_err(PipelineError::ScalarWrite)?;

    let files_written = bands.len() as u32;
    Ok(AnalysisStepStats {
        score_files_written: files_written,
        point_count,
        band_count: bands.len() as u32,
        profile_cells_written: prepared_arc.len() as u32,
        duct_cells: duct_count,
        nexrad_cells_with_echo,
        commercial_cells_boosted,
    })
}

fn apply_nexrad(grid: &mut PointGrid, obs: &[NexradObservation]) -> u32 {
    let mut with_echo = 0u32;
    for o in obs {
        let key = (
            (o.lat * 1000.0).round() as i32,
            (o.lon * 1000.0).round() as i32,
        );
        // merged's keys are whatever extract_grid produced; use direct
        // lat/lon match by scanning. This is O(N*M) worst-case but N is
        // ~92k and obs mirrors the same point list, so in practice each
        // observation matches exactly one cell by int-key.
        if let Some(cell) = grid.get_mut(&key) {
            if o.max_reflectivity_dbz > 0.0 {
                with_echo += 1;
            }
            cell.insert(
                "nexrad_max_reflectivity_dbz".into(),
                o.max_reflectivity_dbz as f32,
            );
        }
    }
    with_echo
}

fn apply_commercial(grid: &mut PointGrid, lookup: &[LinkLookupEntry]) -> u32 {
    let mut boosted = 0u32;
    for (key, cell) in grid.iter_mut() {
        let (lat, lon) = decoder::key_to_latlon(*key);
        if let Some(d) = commercial::degradation_at((lat, lon), lookup) {
            cell.insert("commercial_degradation_db".into(), d.degradation_db as f32);
            cell.insert("commercial_baseline_dbm".into(), d.baseline_dbm as f32);
            cell.insert("commercial_current_dbm".into(), d.current_dbm as f32);
            cell.insert("commercial_n_links".into(), d.n_links as f32);
            boosted += 1;
        }
    }
    boosted
}

/// Turn one enriched cell into a ProfilesFile entry. The msgpack
/// `profile` map keeps string keys; the Elixir reader atomizes the
/// whitelisted subset (surface_*, hpbl_m, pwat_mm, duct metrics, …).
fn cell_to_profile_entry(lat: f64, lon: f64, cell: &CellValues) -> CellEntry {
    use profiles_file::value_map;
    use rmpv::Value as V;

    let mut kv: Vec<(&'static str, V)> = Vec::with_capacity(16);

    if let Some(&v) = cell.get("TMP:2 m above ground") {
        kv.push(("surface_temp_c", V::F64((v as f64) - 273.15)));
    }
    if let Some(&v) = cell.get("DPT:2 m above ground") {
        kv.push(("surface_dewpoint_c", V::F64((v as f64) - 273.15)));
    }
    if let Some(&v) = cell.get("PRES:surface") {
        kv.push(("surface_pressure_mb", V::F64((v as f64) / 100.0)));
    }
    if let Some(&v) = cell.get("HPBL:surface") {
        kv.push(("hpbl_m", V::F64(v as f64)));
    }
    if let Some(&v) = cell.get("PWAT:entire atmosphere (considered as a single layer)") {
        kv.push(("pwat_mm", V::F64(v as f64)));
    }
    // Wind / cloud / precip are read by Elixir's per-QSO scorer
    // (`Propagation.factors_for` in propagation.ex) when a user clicks a
    // grid cell on /map. Without these, wind/sky/rain factors silently
    // score 0 for Rust-produced ProfilesFile cells.
    if let Some(&u) = cell.get("UGRD:10 m above ground") {
        kv.push(("wind_u", V::F64(u as f64)));
    }
    if let Some(&vv) = cell.get("VGRD:10 m above ground") {
        kv.push(("wind_v", V::F64(vv as f64)));
    }
    if let Some(&c) = cell.get("TCDC:entire atmosphere") {
        kv.push(("cloud_cover_pct", V::F64(c as f64)));
    }
    if let Some(&p) = cell.get("APCP:surface") {
        kv.push(("precip_mm", V::F64(p as f64)));
    }
    if let Some(&v) = cell.get("native_min_gradient") {
        kv.push(("native_min_gradient", V::F64(v as f64)));
    }
    if let Some(&v) = cell.get("best_duct_freq_ghz") {
        kv.push(("best_duct_freq_ghz", V::F64(v as f64)));
    }
    if let Some(&v) = cell.get("max_duct_thickness_m") {
        kv.push(("max_duct_thickness_m", V::F64(v as f64)));
    }
    if let Some(&v) = cell.get("duct_count") {
        kv.push(("duct_count", V::Integer((v as i64).into())));
    }
    if let Some(&v) = cell.get("nexrad_max_reflectivity_dbz") {
        kv.push(("nexrad_max_reflectivity_dbz", V::F64(v as f64)));
    }
    if let Some(&d) = cell.get("commercial_degradation_db") {
        let links = cell.get("commercial_n_links").copied().unwrap_or(0.0);
        let baseline = cell.get("commercial_baseline_dbm").copied().unwrap_or(0.0);
        let current = cell.get("commercial_current_dbm").copied().unwrap_or(0.0);
        kv.push((
            "commercial_link_degradation",
            value_map([
                ("degradation_db", V::F64(d as f64)),
                ("baseline_dbm", V::F64(baseline as f64)),
                ("current_dbm", V::F64(current as f64)),
                ("n_links", V::Integer((links as i64).into())),
            ]),
        ));
    }

    // Pressure-level profile list for SoundingParams.derive on the
    // Elixir side. Mirrors the fields derive expects.
    let levels: Vec<V> = fetcher::grid_level_keys()
        .iter()
        .filter_map(|k| {
            let t = cell.get(k.tmp.as_str()).copied()?;
            let d = cell.get(k.dpt.as_str()).copied();
            let h = cell.get(k.hgt.as_str()).copied()?;
            let mut pairs: Vec<(V, V)> = Vec::with_capacity(4);
            pairs.push((V::String("pres_mb".into()), V::F64(k.pres_mb)));
            pairs.push((V::String("hght_m".into()), V::F64(h as f64)));
            pairs.push((V::String("tmpc".into()), V::F64((t as f64) - 273.15)));
            if let Some(dv) = d {
                pairs.push((V::String("dwpc".into()), V::F64((dv as f64) - 273.15)));
            }
            Some(V::Map(pairs))
        })
        .collect();
    if !levels.is_empty() {
        kv.push(("profile", V::Array(levels)));
    }

    let profile = value_map(kv);

    CellEntry { lat, lon, profile }
}