#include "tsunami_core.h"

#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void tsunami_fatal(const char *message) {
    fprintf(stderr, "Error: %s\n", message);
    exit(EXIT_FAILURE);
}

void tsunami_fatal_errno(const char *message) {
    fprintf(stderr, "Error: %s: %s\n", message, strerror(errno));
    exit(EXIT_FAILURE);
}

static int read_exact(FILE *fp, void *dst, size_t bytes) {
    return fread(dst, 1, bytes, fp) == bytes;
}

void *tsunami_xcalloc(size_t count, size_t size) {
    void *ptr = calloc(count, size);
    if (!ptr) tsunami_fatal("out of memory");
    return ptr;
}

static float g_source_ramp_seconds = SOURCE_RAMP_DEFAULT_SECONDS;

void set_source_ramp_seconds(float seconds) {
    if (seconds < 1.0f) seconds = 1.0f;
    if (seconds > 3600.0f) seconds = 3600.0f;
    g_source_ramp_seconds = seconds;
}

float source_ramp_seconds(void) {
    return g_source_ramp_seconds;
}

static float clampf_local(float x, float lo, float hi) {
    if (x < lo) return lo;
    if (x > hi) return hi;
    return x;
}

float visible_effective_depth_m(const GridData *g, int visible_idx);

/*
   Original PenguinMod lists:
     관측점   -> station grid index list
     station_id    -> JMA tsunami forecast zone id for each station point
     예보구 명 -> JMA forecast zone names, ids 1..66

   In the original project, station_id value 67 is used for points outside the
   official Japanese forecast-zone list. Most of these points are kept for GUI
   display only and are not exported to station or forecast-zone CSV files.

   The vertical north-south line at the far right edge of the exported grid is
   an artificial boundary observation line, not a real coastal observation set.
   It is hidden/ignored through station_enabled[] without deleting rows, so the
   station_index <-> station_id mapping never shifts.
*/
static const char *JMA_FORECAST_ZONE_NAMES[66] = {
    "オホーツク海沿岸",
    "北海道太平洋沿岸東部",
    "北海道太平洋沿岸中部",
    "北海道太平洋沿岸西部",
    "北海道日本海沿岸北部",
    "北海道日本海沿岸南部",
    "陸奥湾",
    "青森県日本海沿岸",
    "青森県太平洋沿岸",
    "岩手県",
    "宮城県",
    "福島県",
    "秋田県",
    "山形県",
    "茨城県",
    "千葉県九十九里・外房",
    "千葉県内房",
    "東京湾内湾",
    "伊豆諸島",
    "小笠原諸島",
    "相模湾・三浦半島",
    "静岡県",
    "愛知県外海",
    "伊勢・三河湾",
    "三重県南部",
    "新潟県上中下越",
    "佐渡",
    "富山県",
    "石川県能登",
    "石川県加賀",
    "福井県",
    "京都府",
    "兵庫県北部",
    "兵庫県瀬戸内海沿岸",
    "淡路島南部",
    "大阪府",
    "和歌山県",
    "鳥取県",
    "島根県出雲・石見",
    "隠岐",
    "岡山県",
    "広島県",
    "香川県",
    "愛媛県瀬戸内海沿岸",
    "愛媛県宇和海沿岸",
    "徳島県",
    "高知県",
    "山口県瀬戸内海沿岸",
    "山口県日本海沿岸",
    "福岡県瀬戸内海沿岸",
    "福岡県日本海沿岸",
    "佐賀県北部",
    "長崎県西方",
    "壱岐・対馬",
    "有明・八代海",
    "熊本県天草灘沿岸",
    "大分県瀬戸内海沿岸",
    "大分県豊後水道沿岸",
    "宮崎県",
    "鹿児島県東部",
    "鹿児島県西部",
    "種子島・屋久島地方",
    "奄美群島・トカラ列島",
    "沖縄本島地方",
    "大東島地方",
    "宮古島・八重山地方",
};


const char *forecast_zone_name_from_id(int region_id) {
    if (region_id >= 1 && region_id <= 66) return JMA_FORECAST_ZONE_NAMES[region_id - 1];
    if (region_id == 67) return "outside_forecast_zone";
    if (region_id == 68) return "강원 북부 해안";
    if (region_id == 69) return "강원 중부 해안";
    if (region_id == 70) return "강원 남부 해안";
    if (region_id == 71) return "경북 북부 해안";
    if (region_id == 72) return "경북 남부 해안";
    if (region_id == 73) return "울산 해안";
    if (region_id == 74) return "부산 해안";
    if (region_id == 75) return "거제시 동부 해안";
    if (region_id == 76) return "경남 중부 남해 해안";
    if (region_id == 77) return "경남 서부 남해 해안";
    if (region_id == 78) return "제주도 북부 해안";
    if (region_id == 79) return "제주도 동부 해안";
    if (region_id == 80) return "제주도 남부 해안";
    if (region_id == 81) return "제주도 서부 해안";
    if (region_id == 82) return "전남 동부 남해 해안";
    if (region_id == 83) return "전남 서부 남해 해안";
    if (region_id == 84) return "전남 남부 서해 해안";
    if (region_id == 85) return "전남 중부 서해 해안";
    if (region_id == 86) return "전남 북부 서해 해안";
    if (region_id == 87) return "전북 남부 서해 해안";
    if (region_id == 88) return "전북 북부 서해 해안";
    if (region_id == 89) return "충남 남부 서해 해안";
    if (region_id == 90) return "충남 북부 서해 해안";
    if (region_id == 91) return "인천.경기 남부 해안";
    if (region_id == 92) return "인천.경기 북부 해안";
    if (region_id == 93) return "울릉도 해안";
    return "unknown";
}

const char *tsunami_height_color_bin(float height_m) {
    /* Station / forecast-zone warning colors use only positive water-level rise.
       Negative drawdown from the subsidence side of a fault is not counted as
       tsunami height. */
    float h = height_m > 0.0f ? height_m : 0.0f;
    if (h < 0.05f) return "0-5cm";
    if (h < 0.10f) return "5-10cm";
    if (h < 0.20f) return "10-20cm";
    if (h < 0.50f) return "20-50cm";
    if (h < 1.00f) return "50-100cm";
    if (h < 2.00f) return "100-200cm";
    if (h < 3.00f) return "200-300cm";
    if (h < 4.00f) return "300-400cm";
    if (h < 5.00f) return "400-500cm";
    if (h < 8.00f) return "500-800cm";
    if (h < 10.0f) return "800-1000cm";
    return "1000cm+";
}

const char *tsunami_height_color_hex(float height_m) {
    float h = height_m > 0.0f ? height_m : 0.0f;
    /* These hex colors are copied from the original station/warning SVG costumes. */
    if (h < 0.05f) return "none";
    if (h < 0.10f) return "#99ff80";
    if (h < 0.20f) return "#00ff00";
    if (h < 0.50f) return "#fff280";
    if (h < 1.00f) return "#ffe600";
    if (h < 2.00f) return "#ff8080";
    if (h < 3.00f) return "#ff0000";
    if (h < 4.00f) return "#e580ff";
    if (h < 5.00f) return "#cc00ff";
    if (h < 8.00f) return "#8b00ad";
    if (h < 10.0f) return "#460057";
    return "#000000";
}

static void csv_write_escaped(FILE *fp, const char *text) {
    fputc('"', fp);
    for (const char *p = text; *p; ++p) {
        if (*p == '"') fputc('"', fp);
        fputc(*p, fp);
    }
    fputc('"', fp);
}

float shoaling_factor(float depth_m) {
    float d = clampf_local(depth_m, 1.0f, 7000.0f);
    /* Green's-law shoreline-equivalent multiplier.
       Keep this function as the physical depth^-1/4 conversion used by
       station/JESIM 1 m gauge equivalence.  Do not call it directly for the
       displayed wave field, otherwise the solver's bathymetry response is
       counted twice. */
    return 9.145f * powf(d, -0.25f);
}

float output_shoaling_factor(float depth_m) {
#if WAVE_OUTPUT_SHOALING_ENABLED
    return shoaling_factor(depth_m);
#else
    (void)depth_m;
    return 1.0f;
#endif
}

float station_shoreline_equiv_factor(float source_depth_m) {
#if STATION_SHORELINE_EQUIV_SHOALING_ENABLED
    float local = shoaling_factor(source_depth_m);
    if (local < 0.0001f) local = 0.0001f;
    return shoaling_factor(STATION_OBSERVATION_DEPTH_M) / local;
#else
    (void)source_depth_m;
    return 1.0f;
#endif
}

static float wave_velocity_gain_from_depth(float depth_m) {
    float d = clampf_local(depth_m, 1.0f, 7000.0f);
    /*
       The height solver uses this value as a wave-equation coefficient, i.e.
       roughly c^2.  Therefore the actual propagation speed is proportional to
       sqrt(velocity_gain).

       Old: velocity_gain = sqrt(depth) / 230
            => actual speed ~ depth^(1/4)

       v85: velocity_gain = depth / 18000
            => correct depth exponent, but the whole field was about 2x slow.

       v86 tried depth / 4500 to make speed ~2x faster, but that also injected
       about 4x more wave-equation acceleration and destabilized/over-amplified
       wave height.

       v87 restored the stable v85 coefficient and used two internal numerical
       updates per 30-second displayed tick.

       v125 keeps this stable coefficient but makes the actual model tick 15 s,
       with one numerical update per tick.  This keeps the same effective
       numerical cadence while making delay/source/observation timing explicit.
    */
    return d / 18000.0f;
}

static float wave_face_depth_m(float center_depth_m, float neighbor_depth_m) {
    float d0 = clampf_local(center_depth_m, 1.0f, 7000.0f);
    float d1 = clampf_local(neighbor_depth_m, 1.0f, 7000.0f);
#if WAVE_FACE_DEPTH_MODE_GEOMETRIC
    /* v137: use the interface depth for the flux coefficient.
       This is not a shallow-depth floor; a uniformly shallow cell still uses
       its real shallow depth.  Only abrupt deep-to-shallow interfaces get an
       intermediate coupling, reducing the artificial wall/reflection produced
       by the old center-depth Laplacian form. */
    return sqrtf(d0 * d1);
#else
    return 0.5f * (d0 + d1);
#endif
}

static float wave_velocity_gain_from_face_depth(float center_depth_m, float neighbor_depth_m) {
    return wave_velocity_gain_from_depth(wave_face_depth_m(center_depth_m, neighbor_depth_m));
}

float global_height_decay_per_step(void) {
    /* v125: global decay is time-based.  With DEFAULT_DT_SECONDS=15 and one
       numerical update per tick this is the exact 15-second keep factor.  If
       internal substeps are re-enabled later, divide the displayed dt by the
       substep count so the product over the displayed tick remains correct. */
    float substep_dt_s = (float)DEFAULT_DT_SECONDS / (float)WAVE_INTERNAL_SUBSTEPS;
    return powf(GLOBAL_HEIGHT_DECAY_KEEP_6H, substep_dt_s / GLOBAL_HEIGHT_DECAY_PERIOD_S);
}

static int g_coast_reflection_absorb_enabled = COAST_REFLECTION_DEFAULT_ENABLED;
static float g_coast_reflection_energy_loss = COAST_REFLECTION_DEFAULT_ENERGY_LOSS;

void set_coastal_reflection_absorption(int enabled, float energy_loss) {
    g_coast_reflection_absorb_enabled = enabled ? 1 : 0;
    if (energy_loss < 0.0f) energy_loss = 0.0f;
    if (energy_loss > 0.95f) energy_loss = 0.95f;
    g_coast_reflection_energy_loss = energy_loss;
}

int coastal_reflection_absorption_enabled(void) {
    return g_coast_reflection_absorb_enabled;
}

float coastal_reflection_absorption_energy_loss(void) {
    return g_coast_reflection_energy_loss;
}

float coastal_reflection_keep_per_step(void) {
    /*
       In this stable boundary version, the coast absorber is not applied to
       the falling limb of a crest.  Instead, solid-land neighbor cells are
       treated as a partially absorbing ghost boundary.

       energy_loss is an approximate reflected-wave energy loss at the coast.
       The corresponding free-surface amplitude reflection coefficient is R:

           R = sqrt(1 - energy_loss)

       R = 1.0 reproduces the old hard-wall coastline.  Smaller R absorbs more
       energy and prevents the coast cell from holding water by damping its
       return-to-zero motion.
    */
    float r = sqrtf(1.0f - g_coast_reflection_energy_loss);
    if (r < 0.05f) r = 0.05f;
    if (r > 1.0f) r = 1.0f;
    return r;
}

static int g_narrow_sea_damping_enabled = NARROW_SEA_DEFAULT_ENABLED;
static float g_narrow_sea_energy_loss = NARROW_SEA_DEFAULT_ENERGY_LOSS;

void set_narrow_sea_damping(int enabled, float energy_loss) {
    g_narrow_sea_damping_enabled = enabled ? 1 : 0;
    if (energy_loss < 0.0f) energy_loss = 0.0f;
    if (energy_loss > 0.90f) energy_loss = 0.90f;
    g_narrow_sea_energy_loss = energy_loss;
}

int narrow_sea_damping_enabled(void) {
    return g_narrow_sea_damping_enabled;
}

float narrow_sea_damping_energy_loss(void) {
    return g_narrow_sea_energy_loss;
}

float narrow_sea_keep_per_step(void) {
    /*
       One-cell channels are below the practical resolution of the 897x540 grid.
       Instead of switching to the slow 2x grid, treat those cells as a short
       unresolved inlet/strait and lose part of the transmitted energy there.
       The user-facing value is energy loss; the solver applies the equivalent
       free-surface amplitude keep factor.
    */
    float r = sqrtf(1.0f - g_narrow_sea_energy_loss);
    if (r < 0.10f) r = 0.10f;
    if (r > 1.0f) r = 1.0f;
    return r;
}

/* Manual one-cell choke for the Shimonoseki/Kanmon strait repair.
   Important: this value is a KEEP ratio, not a loss ratio.
   0.10 means only about 10% of height/velocity remains through the requested
   under-resolved cell, not a 10% loss leaving 90%.

   v125 uses a 15-second model tick.  The calibration target remains the old
   30-second displayed-tick result: about 10% remains after 30 seconds of
   model-time passage through this under-resolved cell.  Therefore each
   numerical substep applies target^(substep_dt/30s). */
#define SHIMONOSEKI_CHOKE_X 332
#define SHIMONOSEKI_CHOKE_Y 225
#define SHIMONOSEKI_CHOKE_TARGET_KEEP_PER_30S 0.10f
#define SHIMONOSEKI_CHOKE_REFERENCE_PERIOD_S 30.0f

static int visible_cell_is_shimonoseki_choke(const GridData *g, int visible_idx) {
    if (!g || visible_idx < 0) return 0;
    int x = visible_idx % (int)g->width;
    int y = visible_idx / (int)g->width;
    return (x == SHIMONOSEKI_CHOKE_X && y == SHIMONOSEKI_CHOKE_Y);
}

static float shimonoseki_choke_keep_per_internal_substep(void) {
    float target = SHIMONOSEKI_CHOKE_TARGET_KEEP_PER_30S;
    if (target < 0.0f) target = 0.0f;
    if (target > 1.0f) target = 1.0f;
    if (target <= 0.0f) return 0.0f;
    float substep_dt_s = (float)DEFAULT_DT_SECONDS / (float)WAVE_INTERNAL_SUBSTEPS;
    float exponent = substep_dt_s / SHIMONOSEKI_CHOKE_REFERENCE_PERIOD_S;
    if (exponent < 0.0f) exponent = 0.0f;
    return powf(target, exponent);
}

static float station_cell_shoreline_equiv_height(const GridData *g, const Simulation *s, int gx, int gy) {
    if (!g || !s) return 0.0f;
    if (gx < 0 || gy < 0 || gx >= (int)g->width || gy >= (int)g->height) return 0.0f;

    size_t idx = (size_t)gy * g->width + (uint32_t)gx;
    if (g->edge_absorb_gate && g->edge_absorb_gate[idx]) return 0.0f;
    if (g->solid[idx]) return 0.0f;

    int sx = gx + (int)s->outer_pad;
    int sy = gy + (int)s->outer_pad;
    if (sx < 0 || sy < 0 || sx >= (int)s->solver_width || sy >= (int)s->solver_height) return 0.0f;

    size_t p = (size_t)sy * (size_t)s->solver_width + (size_t)sx;
    float solver_h = s->height_norm[p];
    if (solver_h >= SCRATCH_LAND_VALUE * 0.5f) return 0.0f;

    /*
       The displayed/max-wave field now reads the solver height directly to avoid
       double-counting shoaling.  Stations are different: they represent a 1 m
       shoreline gauge, while the finite grid may sample a 100 m+ offshore cell.
       Convert each nearby ocean-cell height to the common 1 m shoreline-
       equivalent reference before comparing.
    */
    float local_depth_m = visible_effective_depth_m(g, idx);
    float physical_h = solver_h * output_shoaling_factor(local_depth_m);
    return physical_h * station_shoreline_equiv_factor(local_depth_m);
}

float station_observed_height_m(const GridData *g, const Simulation *s, uint32_t station_i) {
    if (!g || !s || station_i >= g->station_count) return 0.0f;
    if (!station_is_enabled(g, station_i)) return 0.0f;

    int32_t idx = g->station_sample_index ? g->station_sample_index[station_i] : g->station_index[station_i];
    if (idx < 0 || (size_t)idx >= (size_t)g->width * (size_t)g->height) return 0.0f;

    /*
       v150: station gauges now read only their assigned sample cell.
       Do not search a radius and do not take the nearby maximum, because that
       can borrow a higher wave from an adjacent bay/coast cell and overstate the
       observation.  Keep the 1 m shoreline-equivalent conversion for the selected
       cell, since the finite grid sample may still sit in a 50-200 m offshore
       cell because of map resolution.
    */
    int gx = idx % (int)g->width;
    int gy = idx / (int)g->width;
    return station_cell_shoreline_equiv_height(g, s, gx, gy);
}

static int is_right_edge_artifact_station(const GridData *g, uint32_t station_i) {
    if (!g || !g->station_index || !g->station_id || station_i >= g->station_count) return 1;
    int32_t idx = g->station_index[station_i];
    size_t cells = (size_t)g->width * (size_t)g->height;
    if (idx < 0 || (size_t)idx >= cells) return 1;
    int x = idx % (int)g->width;
    /* The exported project has 539 display-only points at x=895 on a 897-wide
       grid.  Treat only that far-right station_id=67 line as an artifact. */
    return (g->station_id[station_i] == 67 && x >= (int)g->width - 2);
}

int station_is_enabled(const GridData *g, uint32_t station_i) {
    if (!g || station_i >= g->station_count) return 0;
    if (!g->station_enabled) return 1;
    return g->station_enabled[station_i] != 0;
}

uint32_t count_enabled_stations(const GridData *g) {
    if (!g) return 0;
    uint32_t count = 0;
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (station_is_enabled(g, i)) count++;
    }
    return count;
}


static void fill_isolated_japan_coast_outside_zones(GridData *g) {
    /*
       Some coastal observation points in the source data are station_id=67
       (outside_forecast_zone) even though they sit alone inside a Japanese
       coastal forecast-zone chain.  Treat only small isolated 67 holes as
       missing data and fill them from nearby JMA zones.

       Conservative rule:
       - candidate must be station_id=67 and not the far-east artificial line
       - within 5 cells, at least 4 nearby JMA-zone stations must exist
       - one JMA zone must be at least 70% of those nearby JMA stations
       - only a small number of nearby 67 points are allowed
    */
    if (!g || !g->station_id || !g->station_index) return;

    enum {
        RADIUS_CELLS = 5,
        RADIUS2 = RADIUS_CELLS * RADIUS_CELLS,
        MIN_JMA_NEIGHBORS = 4,
        MAX_NEARBY_67 = 2
    };

    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (g->station_id[i] != 67) continue;
        if (is_right_edge_artifact_station(g, i)) continue;

        int32_t idx_i = g->station_index[i];
        if (idx_i < 0) continue;
        int xi = idx_i % (int)g->width;
        int yi = idx_i / (int)g->width;

        int counts[67];
        memset(counts, 0, sizeof(counts));
        int jma_neighbors = 0;
        int nearby_67 = 0;
        int nearby_foreign = 0;

        for (uint32_t j = 0; j < g->station_count; ++j) {
            if (i == j) continue;
            int32_t idx_j = g->station_index[j];
            if (idx_j < 0) continue;

            int xj = idx_j % (int)g->width;
            int yj = idx_j / (int)g->width;
            int dx = xi - xj;
            int dy = yi - yj;
            int d2 = dx * dx + dy * dy;
            if (d2 > RADIUS2) continue;

            int zid = g->station_id[j];
            if (zid >= 1 && zid <= 66) {
                counts[zid]++;
                jma_neighbors++;
            } else if (zid == 67) {
                nearby_67++;
            } else if (zid >= 68 && zid <= 100) {
                nearby_foreign++;
            }
        }

        if (jma_neighbors < MIN_JMA_NEIGHBORS) continue;
        if (nearby_67 > MAX_NEARBY_67) continue;
        if (nearby_foreign > jma_neighbors) continue;

        int best_zone = 0;
        int best_count = 0;
        for (int z = 1; z <= 66; ++z) {
            if (counts[z] > best_count) {
                best_count = counts[z];
                best_zone = z;
            }
        }

        if (best_zone > 0 && best_count * 100 >= jma_neighbors * 70) {
            g->station_id[i] = best_zone;
        }
    }
}


static void fill_isolated_chugoku_45_48_inside_hiroshima(GridData *g) {
    /*
       Some Seto/Chugoku boundary points are still encoded as 45 even
       though they are isolated inside a continuous zone-42 Hiroshima cluster.

       Conservative cleanup:
       - candidate must currently be 45
       - within 5 cells, zone 42 must strongly dominate
       - same-zone support for the candidate must be weak
       - then reassign to 42 広島県
    */
    if (!g || !g->station_id || !g->station_index) return;

    enum {
        RADIUS_CELLS = 5,
        RADIUS2 = RADIUS_CELLS * RADIUS_CELLS,
        MIN_NEARBY_42 = 4,
        MAX_SAME_ZONE_SUPPORT = 1
    };

    for (uint32_t i = 0; i < g->station_count; ++i) {
        int zid = g->station_id[i];
        if (zid != 45) continue;

        int32_t idx_i = g->station_index[i];
        if (idx_i < 0) continue;
        int xi = idx_i % (int)g->width;
        int yi = idx_i / (int)g->width;

        int count42 = 0;
        int same_zone = 0;
        int other_jma = 0;

        for (uint32_t j = 0; j < g->station_count; ++j) {
            if (i == j) continue;
            int32_t idx_j = g->station_index[j];
            if (idx_j < 0) continue;

            int xj = idx_j % (int)g->width;
            int yj = idx_j / (int)g->width;
            int dx = xi - xj;
            int dy = yi - yj;
            int d2 = dx * dx + dy * dy;
            if (d2 > RADIUS2) continue;

            int nz = g->station_id[j];
            if (nz == 42) count42++;
            else if (nz == zid) same_zone++;
            else if (nz >= 1 && nz <= 66 && nz != 67) other_jma++;
        }

        if (count42 >= MIN_NEARBY_42 &&
            same_zone <= MAX_SAME_ZONE_SUPPORT &&
            count42 >= other_jma + same_zone) {
            g->station_id[i] = 42;
        }
    }
}


static void apply_station_id_overrides(GridData *g) {
    if (!g || !g->station_id) return;

    /* Manual corrections requested from the original observation-point mapping.
       These station numbers are 1-based, matching the station_no used in CSV
       export and the GUI hover popup.  Do not delete or reorder rows; only
       correct the station_id/forecast-zone id in-place so station_index <-> station_id
       alignment stays unchanged. */
    static const uint32_t yamaguchi_seto_station_no[] = {
        1463u, 1483u, 1500u, 1501u, 1502u, 1524u
    };

    struct StationZoneOverride {
        uint32_t station_no;
        int32_t zone_id;
    };

    static const struct StationZoneOverride single_station_overrides[] = {
        { 2683u, 28 }, /* 富山県 */
        { 1141u, 58 }, /* 大分県豊後水道沿岸 */
        { 1031u, 59 }, /* 宮崎県 */

        /* Boundary/outlier corrections.
           These station numbers are near forecast-zone borders, but their
           original id was an isolated wrong zone. */
        { 1047u, 56 }, /* 熊本県天草灘沿岸: was isolated 壱岐・対馬 near Nagasaki/Amakusa boundary */
        { 2378u, 33 }, /* 兵庫県北部: was 24 */
        { 2649u, 12 }, /* 福島県: was isolated 秋田県 near Fukushima/Ibaraki boundary */

        /* Boundary isolated wrong id near Tokyo Bay/Izu area. */
        { 2209u, 18 }, /* 東京湾内湾: was isolated 小笠原諸島 */
        { 1610u, 35 },
        { 1637u, 35 }, /* 淡路島南部 */
        /* Remaining 67 outside-forecast holes manually filled by local zone. */
        { 1247u, 45 }, /* 愛媛県宇和海沿岸 */
        { 1245u, 44 }, /* 愛媛県瀬戸内海沿岸 */
        { 1266u, 44 }, /* 愛媛県瀬戸内海沿岸 */
        { 3861u, 5 },  /* 北海道日本海沿岸北部 */
        { 3862u, 4 },  /* 北海道太平洋沿岸西部 */

        { 1647u, 83 }, /* 전남 서부 남해 해안 */ /* 淡路島南部: was isolated 隠岐 near Awaji/Tokushima boundary */
        { 2169u, 17 }, /* 千葉県内房: was isolated 伊豆諸島 near Tokyo Bay/Chiba boundary */

        /* Chugoku / Seto Inland Sea split.
           Some points that should continue zone 45 or 48 were encoded as
           zone 42 広島県.  Keep the main Hiroshima cluster as 42, but pull
           the leakage points back into the adjoining official zones. */
        { 1543u, 42 }, /* 広島県: corrected back from 45 */
        { 1559u, 42 }, /* 広島県: corrected back from 45 */
        { 1560u, 42 }, /* 広島県: corrected back from 45 */
        { 1561u, 42 }, /* 広島県: corrected back from 45 */
        { 1562u, 42 }, /* 広島県: corrected back from 45 */
        { 1582u, 42 }, /* 広島県: corrected back from 45 */

        { 1603u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1604u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1628u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1629u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */

        /* Station 4891 is not zone 45.  It sits on the Hokkaido 1/6 boundary;
           use the local surrounding zone instead of the isolated original id. */
        { 4891u, 6 },  /* 北海道日本海沿岸南部 */

        /* Zone 45 愛媛県宇和海沿岸 was missing a few outside-zone(67) points
           on its eastern side.  Extend zone 45 into those adjacent 67 points.
           Do not include the Chugoku/Hiroshima blob here. */
        { 1068u, 45 },
        { 1084u, 45 },
        { 1095u, 45 },

        /* Chugoku/Hiroshima continuous 45-looking blob plus confirmed 48 exceptions.
           These are not isolated single points; they form a small continuous
           45-side patch inside/along the 42 広島県 cluster.  Lock confirmed points to 42.
           For zone 48, only stations 1558 and 1581 are confirmed exceptions. */
        { 1431u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1432u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1433u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1450u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1469u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1489u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1490u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1491u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1503u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1504u, 42 },
        { 1505u, 42 },
        { 1525u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1526u, 42 },
        { 1527u, 42 },
        { 1541u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1542u, 42 },
        { 1557u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1558u, 42 },
        { 1580u, 48 }, /* 山口県瀬戸内海沿岸: restored; only 1558/1581 should become 42 */
        { 1581u, 42 }
    };


    /* 1st application of Korea Meteorological Administration official coastal
       forecast-zone assignments provided via zone_assignments*.csv.
       Station numbers are 1-based and match GUI/CSV station_no. */
    static const struct StationZoneOverride kma_official_zone_assignments[] = {
        { 1182u, 81 }, /* 제주도 서부 해안 */
        { 1183u, 80 }, /* 제주도 남부 해안 */
        { 1184u, 80 }, /* 제주도 남부 해안 */
        { 1185u, 80 }, /* 제주도 남부 해안 */
        { 1186u, 80 }, /* 제주도 남부 해안 */
        { 1193u, 81 }, /* 제주도 서부 해안 */
        { 1194u, 80 }, /* 제주도 남부 해안 */
        { 1195u, 80 }, /* 제주도 남부 해안 */
        { 1196u, 80 }, /* 제주도 남부 해안 */
        { 1197u, 80 }, /* 제주도 남부 해안 */
        { 1198u, 80 }, /* 제주도 남부 해안 */
        { 1207u, 81 }, /* 제주도 서부 해안 */
        { 1208u, 80 }, /* 제주도 남부 해안 */
        { 1209u, 79 }, /* 제주도 동부 해안 */
        { 1221u, 81 }, /* 제주도 서부 해안 */
        { 1222u, 79 }, /* 제주도 동부 해안 */
        { 1238u, 81 }, /* 제주도 서부 해안 */
        { 1239u, 79 }, /* 제주도 동부 해안 */
        { 1258u, 81 }, /* 제주도 서부 해안 */
        { 1259u, 79 }, /* 제주도 동부 해안 */
        { 1274u, 81 }, /* 제주도 서부 해안 */
        { 1275u, 79 }, /* 제주도 동부 해안 */
        { 1290u, 81 }, /* 제주도 서부 해안 */
        { 1291u, 81 }, /* 제주도 서부 해안 */
        { 1292u, 79 }, /* 제주도 동부 해안 */
        { 1311u, 81 }, /* 제주도 서부 해안 */
        { 1312u, 81 }, /* 제주도 서부 해안 */
        { 1313u, 79 }, /* 제주도 동부 해안 */
        { 1324u, 78 }, /* 제주도 북부 해안 */
        { 1325u, 78 }, /* 제주도 북부 해안 */
        { 1326u, 78 }, /* 제주도 북부 해안 */
        { 1327u, 78 }, /* 제주도 북부 해안 */
        { 1328u, 79 }, /* 제주도 동부 해안 */
        { 1342u, 78 }, /* 제주도 북부 해안 */
        { 1343u, 78 }, /* 제주도 북부 해안 */
        { 1344u, 78 }, /* 제주도 북부 해안 */
        { 1345u, 78 }, /* 제주도 북부 해안 */
        { 1519u, 83 }, /* 전남 서부 남해 해안 */
        { 1520u, 83 }, /* 전남 서부 남해 해안 */
        { 1536u, 83 }, /* 전남 서부 남해 해안 */
        { 1537u, 83 }, /* 전남 서부 남해 해안 */
        { 1552u, 83 }, /* 전남 서부 남해 해안 */
        { 1553u, 83 }, /* 전남 서부 남해 해안 */
        { 1574u, 83 }, /* 전남 서부 남해 해안 */
        { 1575u, 83 }, /* 전남 서부 남해 해안 */
        { 1576u, 83 }, /* 전남 서부 남해 해안 */
        { 1596u, 83 }, /* 전남 서부 남해 해안 */
        { 1597u, 83 }, /* 전남 서부 남해 해안 */
        { 1598u, 83 }, /* 전남 서부 남해 해안 */
        { 1599u, 83 }, /* 전남 서부 남해 해안 */
        { 1622u, 83 }, /* 전남 서부 남해 해안 */
        { 1623u, 83 }, /* 전남 서부 남해 해안 */
        { 1624u, 83 }, /* 전남 서부 남해 해안 */
        { 1643u, 83 }, /* 전남 서부 남해 해안 */
        { 1644u, 83 }, /* 전남 서부 남해 해안 */
        { 1645u, 83 }, /* 전남 서부 남해 해안 */
        { 1646u, 83 }, /* 전남 서부 남해 해안 */
        { 1669u, 84 }, /* 전남 남부 서해 해안 */
        { 1670u, 83 }, /* 전남 서부 남해 해안 */
        { 1671u, 83 }, /* 전남 서부 남해 해안 */
        { 1672u, 83 }, /* 전남 서부 남해 해안 */
        { 1673u, 83 }, /* 전남 서부 남해 해안 */
        { 1674u, 82 }, /* 전남 동부 남해 해안 */
        { 1675u, 82 }, /* 전남 동부 남해 해안 */
        { 1694u, 84 }, /* 전남 남부 서해 해안 */
        { 1695u, 83 }, /* 전남 서부 남해 해안 */
        { 1696u, 83 }, /* 전남 서부 남해 해안 */
        { 1697u, 83 }, /* 전남 서부 남해 해안 */
        { 1698u, 83 }, /* 전남 서부 남해 해안 */
        { 1699u, 83 }, /* 전남 서부 남해 해안 */
        { 1700u, 83 }, /* 전남 서부 남해 해안 */
        { 1701u, 82 }, /* 전남 동부 남해 해안 */
        { 1702u, 82 }, /* 전남 동부 남해 해안 */
        { 1703u, 82 }, /* 전남 동부 남해 해안 */
        { 1704u, 82 }, /* 전남 동부 남해 해안 */
        { 1725u, 84 }, /* 전남 남부 서해 해안 */
        { 1726u, 83 }, /* 전남 서부 남해 해안 */
        { 1727u, 83 }, /* 전남 서부 남해 해안 */
        { 1728u, 83 }, /* 전남 서부 남해 해안 */
        { 1729u, 82 }, /* 전남 동부 남해 해안 */
        { 1730u, 82 }, /* 전남 동부 남해 해안 */
        { 1757u, 84 }, /* 전남 남부 서해 해안 */
        { 1758u, 84 }, /* 전남 남부 서해 해안 */
        { 1759u, 83 }, /* 전남 서부 남해 해안 */
        { 1760u, 82 }, /* 전남 동부 남해 해안 */
        { 1761u, 82 }, /* 전남 동부 남해 해안 */
        { 1762u, 82 }, /* 전남 동부 남해 해안 */
        { 1763u, 82 }, /* 전남 동부 남해 해안 */
        { 1764u, 82 }, /* 전남 동부 남해 해안 */
        { 1777u, 84 }, /* 전남 남부 서해 해안 */
        { 1778u, 82 }, /* 전남 동부 남해 해안 */
        { 1779u, 82 }, /* 전남 동부 남해 해안 */
        { 1780u, 82 }, /* 전남 동부 남해 해안 */
        { 1781u, 82 }, /* 전남 동부 남해 해안 */
        { 1782u, 77 }, /* 경남 서부 남해 해안 */
        { 1800u, 84 }, /* 전남 남부 서해 해안 */
        { 1801u, 82 }, /* 전남 동부 남해 해안 */
        { 1802u, 82 }, /* 전남 동부 남해 해안 */
        { 1803u, 82 }, /* 전남 동부 남해 해안 */
        { 1804u, 82 }, /* 전남 동부 남해 해안 */
        { 1805u, 82 }, /* 전남 동부 남해 해안 */
        { 1806u, 82 }, /* 전남 동부 남해 해안 */
        { 1807u, 77 }, /* 경남 서부 남해 해안 */
        { 1840u, 84 }, /* 전남 남부 서해 해안 */
        { 1841u, 84 }, /* 전남 남부 서해 해안 */
        { 1842u, 84 }, /* 전남 남부 서해 해안 */
        { 1843u, 84 }, /* 전남 남부 서해 해안 */
        { 1844u, 82 }, /* 전남 동부 남해 해안 */
        { 1845u, 82 }, /* 전남 동부 남해 해안 */
        { 1846u, 82 }, /* 전남 동부 남해 해안 */
        { 1847u, 82 }, /* 전남 동부 남해 해안 */
        { 1848u, 77 }, /* 경남 서부 남해 해안 */
        { 1876u, 84 }, /* 전남 남부 서해 해안 */
        { 1877u, 84 }, /* 전남 남부 서해 해안 */
        { 1878u, 84 }, /* 전남 남부 서해 해안 */
        { 1879u, 84 }, /* 전남 남부 서해 해안 */
        { 1880u, 84 }, /* 전남 남부 서해 해안 */
        { 1881u, 84 }, /* 전남 남부 서해 해안 */
        { 1882u, 82 }, /* 전남 동부 남해 해안 */
        { 1883u, 82 }, /* 전남 동부 남해 해안 */
        { 1884u, 77 }, /* 경남 서부 남해 해안 */
        { 1885u, 77 }, /* 경남 서부 남해 해안 */
        { 1886u, 77 }, /* 경남 서부 남해 해안 */
        { 1888u, 75 }, /* 거제시 동부 해안 */
        { 1908u, 85 }, /* 전남 중부 서해 해안 */
        { 1909u, 85 }, /* 전남 중부 서해 해안 */
        { 1910u, 85 }, /* 전남 중부 서해 해안 */
        { 1911u, 85 }, /* 전남 중부 서해 해안 */
        { 1912u, 85 }, /* 전남 중부 서해 해안 */
        { 1913u, 82 }, /* 전남 동부 남해 해안 */
        { 1914u, 82 }, /* 전남 동부 남해 해안 */
        { 1915u, 82 }, /* 전남 동부 남해 해안 */
        { 1916u, 82 }, /* 전남 동부 남해 해안 */
        { 1917u, 77 }, /* 경남 서부 남해 해안 */
        { 1918u, 77 }, /* 경남 서부 남해 해안 */
        { 1919u, 77 }, /* 경남 서부 남해 해안 */
        { 1920u, 76 }, /* 경남 중부 남해 해안 */
        { 1921u, 75 }, /* 거제시 동부 해안 */
        { 1948u, 85 }, /* 전남 중부 서해 해안 */
        { 1949u, 85 }, /* 전남 중부 서해 해안 */
        { 1950u, 85 }, /* 전남 중부 서해 해안 */
        { 1951u, 85 }, /* 전남 중부 서해 해안 */
        { 1952u, 85 }, /* 전남 중부 서해 해안 */
        { 1953u, 82 }, /* 전남 동부 남해 해안 */
        { 1954u, 82 }, /* 전남 동부 남해 해안 */
        { 1955u, 77 }, /* 경남 서부 남해 해안 */
        { 1956u, 77 }, /* 경남 서부 남해 해안 */
        { 1957u, 76 }, /* 경남 중부 남해 해안 */
        { 1958u, 75 }, /* 거제시 동부 해안 */
        { 1980u, 85 }, /* 전남 중부 서해 해안 */
        { 1981u, 85 }, /* 전남 중부 서해 해안 */
        { 1982u, 85 }, /* 전남 중부 서해 해안 */
        { 1983u, 85 }, /* 전남 중부 서해 해안 */
        { 1984u, 85 }, /* 전남 중부 서해 해안 */
        { 1985u, 82 }, /* 전남 동부 남해 해안 */
        { 1986u, 77 }, /* 경남 서부 남해 해안 */
        { 1987u, 77 }, /* 경남 서부 남해 해안 */
        { 1988u, 76 }, /* 경남 중부 남해 해안 */
        { 1989u, 75 }, /* 거제시 동부 해안 */
        { 2001u, 85 }, /* 전남 중부 서해 해안 */
        { 2002u, 85 }, /* 전남 중부 서해 해안 */
        { 2003u, 85 }, /* 전남 중부 서해 해안 */
        { 2004u, 77 }, /* 경남 서부 남해 해안 */
        { 2005u, 76 }, /* 경남 중부 남해 해안 */
        { 2006u, 76 }, /* 경남 중부 남해 해안 */
        { 2007u, 76 }, /* 경남 중부 남해 해안 */
        { 2008u, 76 }, /* 경남 중부 남해 해안 */
        { 2009u, 75 }, /* 거제시 동부 해안 */
        { 2021u, 85 }, /* 전남 중부 서해 해안 */
        { 2022u, 85 }, /* 전남 중부 서해 해안 */
        { 2023u, 85 }, /* 전남 중부 서해 해안 */
        { 2029u, 77 }, /* 경남 서부 남해 해안 */
        { 2030u, 77 }, /* 경남 서부 남해 해안 */
        { 2031u, 76 }, /* 경남 중부 남해 해안 */
        { 2032u, 75 }, /* 거제시 동부 해안 */
        { 2044u, 85 }, /* 전남 중부 서해 해안 */
        { 2045u, 85 }, /* 전남 중부 서해 해안 */
        { 2046u, 85 }, /* 전남 중부 서해 해안 */
        { 2051u, 76 }, /* 경남 중부 남해 해안 */
        { 2052u, 75 }, /* 거제시 동부 해안 */
        { 2065u, 85 }, /* 전남 중부 서해 해안 */
        { 2066u, 85 }, /* 전남 중부 서해 해안 */
        { 2067u, 85 }, /* 전남 중부 서해 해안 */
        { 2070u, 76 }, /* 경남 중부 남해 해안 */
        { 2071u, 76 }, /* 경남 중부 남해 해안 */
        { 2072u, 75 }, /* 거제시 동부 해안 */
        { 2084u, 86 }, /* 전남 북부 서해 해안 */
        { 2085u, 86 }, /* 전남 북부 서해 해안 */
        { 2086u, 76 }, /* 경남 중부 남해 해안 */
        { 2087u, 76 }, /* 경남 중부 남해 해안 */
        { 2088u, 76 }, /* 경남 중부 남해 해안 */
        { 2089u, 75 }, /* 거제시 동부 해안 */
        { 2090u, 74 }, /* 부산 해안 */
        { 2091u, 74 }, /* 부산 해안 */
        { 2103u, 86 }, /* 전남 북부 서해 해안 */
        { 2104u, 86 }, /* 전남 북부 서해 해안 */
        { 2105u, 86 }, /* 전남 북부 서해 해안 */
        { 2106u, 86 }, /* 전남 북부 서해 해안 */
        { 2107u, 76 }, /* 경남 중부 남해 해안 */
        { 2108u, 76 }, /* 경남 중부 남해 해안 */
        { 2111u, 75 }, /* 거제시 동부 해안 */
        { 2112u, 75 }, /* 거제시 동부 해안 */
        { 2113u, 74 }, /* 부산 해안 */
        { 2114u, 74 }, /* 부산 해안 */
        { 2115u, 74 }, /* 부산 해안 */
        { 2127u, 86 }, /* 전남 북부 서해 해안 */
        { 2128u, 86 }, /* 전남 북부 서해 해안 */
        { 2129u, 86 }, /* 전남 북부 서해 해안 */
        { 2130u, 74 }, /* 부산 해안 */
        { 2131u, 74 }, /* 부산 해안 */
        { 2132u, 74 }, /* 부산 해안 */
        { 2144u, 86 }, /* 전남 북부 서해 해안 */
        { 2145u, 86 }, /* 전남 북부 서해 해안 */
        { 2146u, 86 }, /* 전남 북부 서해 해안 */
        { 2147u, 74 }, /* 부산 해안 */
        { 2148u, 74 }, /* 부산 해안 */
        { 2159u, 86 }, /* 전남 북부 서해 해안 */
        { 2160u, 86 }, /* 전남 북부 서해 해안 */
        { 2161u, 86 }, /* 전남 북부 서해 해안 */
        { 2162u, 74 }, /* 부산 해안 */
        { 2173u, 86 }, /* 전남 북부 서해 해안 */
        { 2174u, 74 }, /* 부산 해안 */
        { 2187u, 86 }, /* 전남 북부 서해 해안 */
        { 2188u, 74 }, /* 부산 해안 */
        { 2204u, 86 }, /* 전남 북부 서해 해안 */
        { 2205u, 73 }, /* 울산 해안 */
        { 2214u, 87 }, /* 전북 남부 서해 해안 */
        { 2215u, 73 }, /* 울산 해안 */
        { 2224u, 87 }, /* 전북 남부 서해 해안 */
        { 2225u, 73 }, /* 울산 해안 */
        { 2234u, 87 }, /* 전북 남부 서해 해안 */
        { 2235u, 73 }, /* 울산 해안 */
        { 2236u, 73 }, /* 울산 해안 */
        { 2245u, 87 }, /* 전북 남부 서해 해안 */
        { 2246u, 73 }, /* 울산 해안 */
        { 2255u, 87 }, /* 전북 남부 서해 해안 */
        { 2256u, 73 }, /* 울산 해안 */
        { 2268u, 87 }, /* 전북 남부 서해 해안 */
        { 2269u, 73 }, /* 울산 해안 */
        { 2293u, 87 }, /* 전북 남부 서해 해안 */
        { 2294u, 73 }, /* 울산 해안 */
        { 2320u, 87 }, /* 전북 남부 서해 해안 */
        { 2321u, 72 }, /* 경북 남부 해안 */
        { 2338u, 87 }, /* 전북 남부 서해 해안 */
        { 2339u, 72 }, /* 경북 남부 해안 */
        { 2349u, 87 }, /* 전북 남부 서해 해안 */
        { 2350u, 87 }, /* 전북 남부 서해 해안 */
        { 2351u, 72 }, /* 경북 남부 해안 */
        { 2368u, 88 }, /* 전북 북부 서해 해안 */
        { 2369u, 72 }, /* 경북 남부 해안 */
        { 2376u, 88 }, /* 전북 북부 서해 해안 */
        { 2377u, 72 }, /* 경북 남부 해안 */
        { 2386u, 88 }, /* 전북 북부 서해 해안 */
        { 2387u, 72 }, /* 경북 남부 해안 */
        { 2393u, 88 }, /* 전북 북부 서해 해안 */
        { 2394u, 88 }, /* 전북 북부 서해 해안 */
        { 2395u, 72 }, /* 경북 남부 해안 */
        { 2400u, 88 }, /* 전북 북부 서해 해안 */
        { 2401u, 72 }, /* 경북 남부 해안 */
        { 2406u, 88 }, /* 전북 북부 서해 해안 */
        { 2407u, 72 }, /* 경북 남부 해안 */
        { 2415u, 88 }, /* 전북 북부 서해 해안 */
        { 2416u, 88 }, /* 전북 북부 서해 해안 */
        { 2417u, 88 }, /* 전북 북부 서해 해안 */
        { 2418u, 72 }, /* 경북 남부 해안 */
        { 2428u, 88 }, /* 전북 북부 서해 해안 */
        { 2429u, 72 }, /* 경북 남부 해안 */
        { 2437u, 89 }, /* 충남 남부 서해 해안 */
        { 2438u, 88 }, /* 전북 북부 서해 해안 */
        { 2439u, 88 }, /* 전북 북부 서해 해안 */
        { 2440u, 72 }, /* 경북 남부 해안 */
        { 2441u, 72 }, /* 경북 남부 해안 */
        { 2442u, 72 }, /* 경북 남부 해안 */
        { 2449u, 89 }, /* 충남 남부 서해 해안 */
        { 2450u, 72 }, /* 경북 남부 해안 */
        { 2451u, 72 }, /* 경북 남부 해안 */
        { 2462u, 89 }, /* 충남 남부 서해 해안 */
        { 2463u, 72 }, /* 경북 남부 해안 */
        { 2471u, 89 }, /* 충남 남부 서해 해안 */
        { 2472u, 72 }, /* 경북 남부 해안 */
        { 2479u, 89 }, /* 충남 남부 서해 해안 */
        { 2480u, 72 }, /* 경북 남부 해안 */
        { 2487u, 89 }, /* 충남 남부 서해 해안 */
        { 2488u, 71 }, /* 경북 북부 해안 */
        { 2498u, 90 }, /* 충남 북부 서해 해안 */
        { 2499u, 71 }, /* 경북 북부 해안 */
        { 2509u, 90 }, /* 충남 북부 서해 해안 */
        { 2510u, 89 }, /* 충남 남부 서해 해안 */
        { 2511u, 90 }, /* 충남 북부 서해 해안 */
        { 2512u, 71 }, /* 경북 북부 해안 */
        { 2519u, 90 }, /* 충남 북부 서해 해안 */
        { 2520u, 90 }, /* 충남 북부 서해 해안 */
        { 2521u, 71 }, /* 경북 북부 해안 */
        { 2526u, 90 }, /* 충남 북부 서해 해안 */
        { 2527u, 90 }, /* 충남 북부 서해 해안 */
        { 2528u, 71 }, /* 경북 북부 해안 */
        { 2533u, 90 }, /* 충남 북부 서해 해안 */
        { 2534u, 90 }, /* 충남 북부 서해 해안 */
        { 2535u, 71 }, /* 경북 북부 해안 */
        { 2541u, 90 }, /* 충남 북부 서해 해안 */
        { 2542u, 90 }, /* 충남 북부 서해 해안 */
        { 2543u, 71 }, /* 경북 북부 해안 */
        { 2549u, 90 }, /* 충남 북부 서해 해안 */
        { 2550u, 90 }, /* 충남 북부 서해 해안 */
        { 2551u, 71 }, /* 경북 북부 해안 */
        { 2557u, 90 }, /* 충남 북부 서해 해안 */
        { 2558u, 90 }, /* 충남 북부 서해 해안 */
        { 2559u, 90 }, /* 충남 북부 서해 해안 */
        { 2560u, 90 }, /* 충남 북부 서해 해안 */
        { 2561u, 71 }, /* 경북 북부 해안 */
        { 2567u, 90 }, /* 충남 북부 서해 해안 */
        { 2568u, 90 }, /* 충남 북부 서해 해안 */
        { 2569u, 90 }, /* 충남 북부 서해 해안 */
        { 2570u, 90 }, /* 충남 북부 서해 해안 */
        { 2571u, 71 }, /* 경북 북부 해안 */
        { 2580u, 90 }, /* 충남 북부 서해 해안 */
        { 2581u, 71 }, /* 경북 북부 해안 */
        { 2586u, 90 }, /* 충남 북부 서해 해안 */
        { 2587u, 71 }, /* 경북 북부 해안 */
        { 2597u, 90 }, /* 충남 북부 서해 해안 */
        { 2598u, 71 }, /* 경북 북부 해안 */
        { 2608u, 90 }, /* 충남 북부 서해 해안 */
        { 2609u, 71 }, /* 경북 북부 해안 */
        { 2623u, 90 }, /* 충남 북부 서해 해안 */
        { 2624u, 90 }, /* 충남 북부 서해 해안 */
        { 2625u, 90 }, /* 충남 북부 서해 해안 */
        { 2626u, 71 }, /* 경북 북부 해안 */
        { 2633u, 90 }, /* 충남 북부 서해 해안 */
        { 2634u, 90 }, /* 충남 북부 서해 해안 */
        { 2635u, 90 }, /* 충남 북부 서해 해안 */
        { 2636u, 71 }, /* 경북 북부 해안 */
        { 2644u, 90 }, /* 충남 북부 서해 해안 */
        { 2645u, 71 }, /* 경북 북부 해안 */
        { 2652u, 90 }, /* 충남 북부 서해 해안 */
        { 2653u, 71 }, /* 경북 북부 해안 */
        { 2662u, 91 }, /* 인천.경기 남부 해안 */
        { 2663u, 91 }, /* 인천.경기 남부 해안 */
        { 2664u, 91 }, /* 인천.경기 남부 해안 */
        { 2665u, 91 }, /* 인천.경기 남부 해안 */
        { 2666u, 91 }, /* 인천.경기 남부 해안 */
        { 2667u, 71 }, /* 경북 북부 해안 */
        { 2676u, 91 }, /* 인천.경기 남부 해안 */
        { 2677u, 91 }, /* 인천.경기 남부 해안 */
        { 2678u, 91 }, /* 인천.경기 남부 해안 */
        { 2679u, 91 }, /* 인천.경기 남부 해안 */
        { 2680u, 71 }, /* 경북 북부 해안 */
        { 2690u, 91 }, /* 인천.경기 남부 해안 */
        { 2691u, 70 }, /* 강원 남부 해안 */
        { 2708u, 91 }, /* 인천.경기 남부 해안 */
        { 2709u, 91 }, /* 인천.경기 남부 해안 */
        { 2710u, 91 }, /* 인천.경기 남부 해안 */
        { 2711u, 70 }, /* 강원 남부 해안 */
        { 2722u, 91 }, /* 인천.경기 남부 해안 */
        { 2723u, 70 }, /* 강원 남부 해안 */
        { 2734u, 91 }, /* 인천.경기 남부 해안 */
        { 2735u, 91 }, /* 인천.경기 남부 해안 */
        { 2736u, 70 }, /* 강원 남부 해안 */
        { 2747u, 91 }, /* 인천.경기 남부 해안 */
        { 2748u, 91 }, /* 인천.경기 남부 해안 */
        { 2749u, 91 }, /* 인천.경기 남부 해안 */
        { 2750u, 91 }, /* 인천.경기 남부 해안 */
        { 2751u, 70 }, /* 강원 남부 해안 */
        { 2762u, 91 }, /* 인천.경기 남부 해안 */
        { 2763u, 91 }, /* 인천.경기 남부 해안 */
        { 2764u, 70 }, /* 강원 남부 해안 */
        { 2773u, 91 }, /* 인천.경기 남부 해안 */
        { 2774u, 91 }, /* 인천.경기 남부 해안 */
        { 2775u, 70 }, /* 강원 남부 해안 */
        { 2786u, 92 }, /* 인천.경기 북부 해안 */
        { 2787u, 91 }, /* 인천.경기 남부 해안 */
        { 2788u, 91 }, /* 인천.경기 남부 해안 */
        { 2789u, 91 }, /* 인천.경기 남부 해안 */
        { 2790u, 70 }, /* 강원 남부 해안 */
        { 2791u, 93 }, /* 울릉도 해안 */
        { 2792u, 93 }, /* 울릉도 해안 */
        { 2808u, 92 }, /* 인천.경기 북부 해안 */
        { 2809u, 91 }, /* 인천.경기 남부 해안 */
        { 2810u, 69 }, /* 강원 중부 해안 */
        { 2811u, 93 }, /* 울릉도 해안 */
        { 2812u, 93 }, /* 울릉도 해안 */
        { 2833u, 92 }, /* 인천.경기 북부 해안 */
        { 2834u, 69 }, /* 강원 중부 해안 */
        { 2835u, 93 }, /* 울릉도 해안 */
        { 2836u, 93 }, /* 울릉도 해안 */
        { 2850u, 92 }, /* 인천.경기 북부 해안 */
        { 2851u, 92 }, /* 인천.경기 북부 해안 */
        { 2852u, 92 }, /* 인천.경기 북부 해안 */
        { 2853u, 92 }, /* 인천.경기 북부 해안 */
        { 2854u, 92 }, /* 인천.경기 북부 해안 */
        { 2855u, 69 }, /* 강원 중부 해안 */
        { 2856u, 93 }, /* 울릉도 해안 */
        { 2857u, 93 }, /* 울릉도 해안 */
        { 2870u, 92 }, /* 인천.경기 북부 해안 */
        { 2871u, 92 }, /* 인천.경기 북부 해안 */
        { 2872u, 92 }, /* 인천.경기 북부 해안 */
        { 2873u, 92 }, /* 인천.경기 북부 해안 */
        { 2874u, 69 }, /* 강원 중부 해안 */
        { 2875u, 93 }, /* 울릉도 해안 */
        { 2876u, 93 }, /* 울릉도 해안 */
        { 2877u, 93 }, /* 울릉도 해안 */
        { 2887u, 92 }, /* 인천.경기 북부 해안 */
        { 2888u, 69 }, /* 강원 중부 해안 */
        { 2898u, 92 }, /* 인천.경기 북부 해안 */
        { 2899u, 92 }, /* 인천.경기 북부 해안 */
        { 2900u, 69 }, /* 강원 중부 해안 */
        { 2914u, 92 }, /* 인천.경기 북부 해안 */
        { 2915u, 69 }, /* 강원 중부 해안 */
        { 2926u, 92 }, /* 인천.경기 북부 해안 */
        { 2927u, 69 }, /* 강원 중부 해안 */
        { 2943u, 69 }, /* 강원 중부 해안 */
        { 2964u, 69 }, /* 강원 중부 해안 */
        { 2981u, 69 }, /* 강원 중부 해안 */
        { 2996u, 69 }, /* 강원 중부 해안 */
        { 3007u, 68 }, /* 강원 북부 해안 */
        { 3019u, 68 }, /* 강원 북부 해안 */
        { 3031u, 68 }, /* 강원 북부 해안 */
        { 3045u, 68 }, /* 강원 북부 해안 */
        { 3058u, 68 }, /* 강원 북부 해안 */
        { 3072u, 68 }, /* 강원 북부 해안 */
        { 3087u, 68 }, /* 강원 북부 해안 */
        { 3095u, 68 }, /* 강원 북부 해안 */
        { 3104u, 68 }, /* 강원 북부 해안 */
        { 3114u, 68 }, /* 강원 북부 해안 */
        { 3126u, 68 }, /* 강원 북부 해안 */
        { 3140u, 68 }, /* 강원 북부 해안 */
        { 3147u, 68 }, /* 강원 북부 해안 */
        { 3155u, 68 }, /* 강원 북부 해안 */
        { 3163u, 68 }, /* 강원 북부 해안 */
        { 3171u, 68 }, /* 강원 북부 해안 */
        { 3179u, 68 }, /* 강원 북부 해안 */
        { 3187u, 68 }, /* 강원 북부 해안 */
    };

    static const struct StationZoneOverride final_user_station_overrides[] = {
        { 1500u, 49 }, /* user correction: station_id / forecast-zone id 49 */
        { 1501u, 49 }, /* user correction: station_id / forecast-zone id 49 */
        { 1502u, 49 }, /* user correction: station_id / forecast-zone id 49 */
        { 1524u, 49 }  /* user correction: station_id / forecast-zone id 49 */
    };

    const int32_t yamaguchi_seto_zone_id = 48; /* 山口県瀬戸内海沿岸 */

    for (size_t k = 0; k < sizeof(yamaguchi_seto_station_no) / sizeof(yamaguchi_seto_station_no[0]); ++k) {
        uint32_t station_no = yamaguchi_seto_station_no[k];
        if (station_no >= 1u && station_no <= g->station_count) {
            g->station_id[station_no - 1u] = yamaguchi_seto_zone_id;
        }
    }

    for (size_t k = 0; k < sizeof(single_station_overrides) / sizeof(single_station_overrides[0]); ++k) {
        uint32_t station_no = single_station_overrides[k].station_no;
        if (station_no >= 1u && station_no <= g->station_count) {
            g->station_id[station_no - 1u] = single_station_overrides[k].zone_id;
        }
    }

    for (size_t k = 0; k < sizeof(kma_official_zone_assignments) / sizeof(kma_official_zone_assignments[0]); ++k) {
        uint32_t station_no = kma_official_zone_assignments[k].station_no;
        if (station_no >= 1u && station_no <= g->station_count) {
            g->station_id[station_no - 1u] = kma_official_zone_assignments[k].zone_id;
        }
    }

    fill_isolated_japan_coast_outside_zones(g);
    fill_isolated_chugoku_45_48_inside_hiroshima(g);

    /* Apply user-requested station-id fixes last so automatic hole-filling cannot
       pull these specific observation points back into a neighboring zone. */
    for (size_t k = 0; k < sizeof(final_user_station_overrides) / sizeof(final_user_station_overrides[0]); ++k) {
        uint32_t station_no = final_user_station_overrides[k].station_no;
        if (station_no >= 1u && station_no <= g->station_count) {
            g->station_id[station_no - 1u] = final_user_station_overrides[k].zone_id;
        }
    }
}

static int is_land_value(float v) {
    return v >= 50.0f;
}

GridData load_grid_data(const char *path) {
    FILE *fp = fopen(path, "rb");
    if (!fp) tsunami_fatal_errno("could not open data file");

    char magic[TSUNAMI_MAGIC_SIZE];
    if (!read_exact(fp, magic, sizeof(magic))) tsunami_fatal("data file is truncated");
    if (memcmp(magic, TSUNAMI_MAGIC, strlen(TSUNAMI_MAGIC)) != 0) {
        tsunami_fatal("data file magic does not match TSUNAMI_SIM_V1");
    }

    GridData g;
    memset(&g, 0, sizeof(g));
    if (!read_exact(fp, &g.width, sizeof(g.width)) ||
        !read_exact(fp, &g.height, sizeof(g.height)) ||
        !read_exact(fp, &g.station_count, sizeof(g.station_count))) {
        tsunami_fatal("data file header is truncated");
    }

    if (g.width == 0 || g.height == 0 || g.width > 10000 || g.height > 10000) {
        tsunami_fatal("invalid grid dimensions in data file");
    }

    size_t cells = (size_t)g.width * (size_t)g.height;
    g.depth_m = (float *)tsunami_xcalloc(cells, sizeof(float));
    g.effective_depth_m = (float *)tsunami_xcalloc(cells, sizeof(float));
    g.base_land = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    g.land = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    g.solid = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    g.user_override = (int8_t *)tsunami_xcalloc(cells, sizeof(int8_t));
    g.station_index = (int32_t *)tsunami_xcalloc(g.station_count, sizeof(int32_t));
    g.station_sample_index = (int32_t *)tsunami_xcalloc(g.station_count, sizeof(int32_t));
    g.station_id = (int32_t *)tsunami_xcalloc(g.station_count, sizeof(int32_t));
    g.station_enabled = (uint8_t *)tsunami_xcalloc(g.station_count, sizeof(uint8_t));
    g.edge_absorb_gate = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    g.edge_absorb_depth_m = (float *)tsunami_xcalloc(cells, sizeof(float));

    /*
       Older exported .bin files contained a saved Scratch buffer as the
       "solid" mask. That buffer can include a large artificial horizontal
       band, which makes the southern side of the map look clipped and still
       leaves real coastline cells unblocked. Read the byte block for binary
       compatibility, but rebuild the actual land/coast mask from DEP below.
    */
    uint8_t *embedded_solid = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));

    if (!read_exact(fp, g.depth_m, cells * sizeof(float)) ||
        !read_exact(fp, embedded_solid, cells * sizeof(uint8_t)) ||
        !read_exact(fp, g.station_index, g.station_count * sizeof(int32_t)) ||
        !read_exact(fp, g.station_id, g.station_count * sizeof(int32_t))) {
        free(embedded_solid);
        tsunami_fatal("data file body is truncated");
    }
    free(embedded_solid);

    apply_station_id_overrides(&g);

    for (uint32_t i = 0; i < g.station_count; ++i) {
        int is_edge_artifact = is_right_edge_artifact_station(&g, i);
        g.station_enabled[i] = is_edge_artifact ? 0u : 1u;

        /* Do not export or draw the artificial north-south observation line at
           the far east edge.  However, keep it as a solver-only "open boundary
           gate" so waves can enter the hidden outer absorbing layer there.
           This avoids the east edge turning into a hard wall after those
           display-only stations are disabled.  The station/station_id arrays are
           not modified, so their index matching never shifts. */
        if (is_edge_artifact) {
            int32_t idx = g.station_index[i];
            if (idx >= 0 && (size_t)idx < cells) {
                int y = idx / (int)g.width;
                /* Mark only the last visible edge cell for this row as an
                   open-boundary gate. Earlier versions used a 10-cell visible
                   handoff strip, but that made the east side behave differently
                   from the south side. Now the actual absorber is outside the
                   visible grid; the visible map itself is not damped. */
                int start_x = (int)g.width - EAST_VISIBLE_GATE_CELLS;
                if (start_x < 0) start_x = 0;
                for (int xx = start_x; xx < (int)g.width; ++xx) {
                    g.edge_absorb_gate[(size_t)y * g.width + (uint32_t)xx] = 1u;
                }
            }
        }
    }

    /* The removed east observation line covers almost the full map height.
       Make the east handoff continuous row-by-row so one missing row does not
       become a hard reflector. This only marks the final edge column as the
       outside-absorber gate; it does not create a 10-cell visible absorber. */
    {
        int start_x = (int)g.width - EAST_VISIBLE_GATE_CELLS;
        if (start_x < 0) start_x = 0;
        for (uint32_t y = 0; y < g.height; ++y) {
            int has_gate_near_row = 0;
            for (uint32_t yy = (y > 2u ? y - 2u : 0u); yy <= y + 2u && yy < g.height; ++yy) {
                if (g.edge_absorb_gate[(size_t)yy * g.width + (g.width - 1u)]) {
                    has_gate_near_row = 1;
                    break;
                }
            }
            if (has_gate_near_row || y + 2u >= g.height) {
                for (int xx = start_x; xx < (int)g.width; ++xx) {
                    g.edge_absorb_gate[(size_t)y * g.width + (uint32_t)xx] = 1u;
                }
            }
        }
    }

    for (size_t i = 0; i < cells; ++i) {
        /* DEP comes from a color heightmap, so it is not perfectly identical to the real coastline. */
        g.base_land[i] = (g.depth_m[i] > 0.0f) ? 0u : 1u;
        g.land[i] = g.base_land[i];
        g.user_override[i] = -1;
    }

    /* Correct the visual DEP mask using the original shoreline/forecast-point list.
       Cells inside the station coastline are land, even if the heightmap color left
       a small positive depth there. The station points themselves are kept as ocean,
       matching the original PenguinMod setup.

       Additional manual corrections can be layered on top later through user_override.
    */
    rebuild_land_mask_from_coast_stations(&g);
    rebuild_solid_mask(&g, DEFAULT_COAST_GUARD_CELLS);

    fclose(fp);
    return g;
}


static void mark_mask_disc(uint8_t *mask, uint32_t width, uint32_t height, int cx, int cy, int radius) {
    if (!mask || radius < 0) return;
    int x0 = cx - radius;
    int x1 = cx + radius;
    int y0 = cy - radius;
    int y1 = cy + radius;
    if (x0 < 0) x0 = 0;
    if (y0 < 0) y0 = 0;
    if (x1 >= (int)width) x1 = (int)width - 1;
    if (y1 >= (int)height) y1 = (int)height - 1;

    for (int y = y0; y <= y1; ++y) {
        for (int x = x0; x <= x1; ++x) {
            int dx = x - cx;
            int dy = y - cy;
            if (dx * dx + dy * dy <= radius * radius) {
                mask[(size_t)y * width + (uint32_t)x] = 1u;
            }
        }
    }
}

void rebuild_land_mask_from_coast_stations(GridData *g) {
    if (!g || !g->land || !g->station_index || !g->base_land) return;

    const uint32_t W = g->width;
    const uint32_t H = g->height;
    const size_t cells = (size_t)W * (size_t)H;
    const int station_barrier_radius = 0;

    uint8_t *blocked = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    uint8_t *visited = (uint8_t *)tsunami_xcalloc(cells, sizeof(uint8_t));
    uint32_t *queue = (uint32_t *)tsunami_xcalloc(cells, sizeof(uint32_t));

    /* Start from the DEP-derived base mask, then apply user overrides before
       the shoreline flood-fill so manual land cells behave as barriers and
       manual water cells can reopen narrow channels. */
    memcpy(blocked, g->base_land, cells * sizeof(uint8_t));
    if (g->user_override) {
        for (size_t i = 0; i < cells; ++i) {
            if (g->user_override[i] > 0) blocked[i] = 1u;
            else if (g->user_override[i] == 0) blocked[i] = 0u;
        }
    }

    /* The original 관측점/station_id coastline is a stronger coastline reference than DEP.
       Use station cells as a TEMPORARY barrier only while deciding what is inland.
       After the flood fill, the station cells are explicitly put back into the ocean.
       station_id=67 is included too because those points represent non-JMA/outside coast
       or map-edge coastline in the original project. */
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        /* station_enabled==2 marks a synthetic/virtual coastal observation
           point.  It is a real observation station for warnings/display, but
           must NOT become a coastline barrier while rebuilding the land mask.
           This lets small-island gauges such as the Izu Islands observe the
           nearest shallow ocean cell without adding new land pixels. */
        if (g->station_enabled && g->station_enabled[i] == 2u) continue;
        int32_t idx = g->station_index[i];
        if (idx < 0 || (size_t)idx >= cells) continue;
        int x = idx % (int)W;
        int y = idx / (int)W;
        mark_mask_disc(blocked, W, H, x, y, station_barrier_radius);
    }

    /* Flood-fill ocean from the map boundary. Anything not reached is inland/land.
       Four-neighbor fill intentionally avoids diagonal leaks through one-pixel
       coastline corners. */
    size_t head = 0, tail = 0;
#define PUSH_IF_WATER(INDEX_EXPR) do {         size_t qi__ = (size_t)(INDEX_EXPR);         if (!blocked[qi__] && !visited[qi__]) {             visited[qi__] = 1u;             queue[tail++] = (uint32_t)qi__;         }     } while (0)

    for (uint32_t x = 0; x < W; ++x) {
        PUSH_IF_WATER(x);
        PUSH_IF_WATER((size_t)(H - 1u) * W + x);
    }
    for (uint32_t y = 0; y < H; ++y) {
        PUSH_IF_WATER((size_t)y * W);
        PUSH_IF_WATER((size_t)y * W + (W - 1u));
    }

    while (head < tail) {
        uint32_t q = queue[head++];
        int x = (int)(q % W);
        int y = (int)(q / W);
        if (x > 0) PUSH_IF_WATER((size_t)y * W + (uint32_t)(x - 1));
        if (x + 1 < (int)W) PUSH_IF_WATER((size_t)y * W + (uint32_t)(x + 1));
        if (y > 0) PUSH_IF_WATER((size_t)(y - 1) * W + (uint32_t)x);
        if (y + 1 < (int)H) PUSH_IF_WATER((size_t)(y + 1) * W + (uint32_t)x);
    }
#undef PUSH_IF_WATER

    for (size_t i = 0; i < cells; ++i) {
        g->land[i] = visited[i] ? 0u : 1u;
    }

    /* Important: 관측점 was used as an ocean-side observation list in the
       original project.  The temporary station barrier above is only for
       separating outside ocean from inland pixels; station cells themselves
       must remain water for source placement, display, and sampling. */
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        int32_t idx = g->station_index[i];
        if (idx < 0 || (size_t)idx >= cells) continue;
        g->land[idx] = 0u;
    }

    /* Re-apply manual overrides last so the operator can fix narrow straits,
       remove false islands, or force tiny embayments to land/water exactly
       as shown on a reference map. */
    if (g->user_override) {
        for (size_t i = 0; i < cells; ++i) {
            if (g->user_override[i] > 0) g->land[i] = 1u;
            else if (g->user_override[i] == 0) g->land[i] = 0u;
        }
    }

    /* Solver-only correction for the removed far-east observation line.
       Those points are hidden from display/export, but the east edge must stay
       connected to the hidden absorber instead of becoming a reflective wall. */
    if (g->edge_absorb_gate) {
        for (size_t i = 0; i < cells; ++i) {
            if (g->edge_absorb_gate[i] && (!g->user_override || g->user_override[i] <= 0)) {
                g->land[i] = 0u;
            }
        }
    }

    free(queue);
    free(visited);
    free(blocked);
}

static void rebuild_station_sample_indices(GridData *g) {
    if (!g || !g->station_index || !g->station_sample_index || !g->solid) return;

    const int max_radius = 16;
    const uint32_t W = g->width;
    const uint32_t H = g->height;
    const size_t cells = (size_t)W * (size_t)H;

    for (uint32_t i = 0; i < g->station_count; ++i) {
        int32_t station_idx = g->station_index[i];
        g->station_sample_index[i] = station_idx;
        if (!station_is_enabled(g, i)) continue;
        if (station_idx < 0 || (size_t)station_idx >= cells) continue;

        int sx = station_idx % (int)W;
        int sy = station_idx / (int)W;
        int best_idx = -1;
        int best_d2 = 1000000000;

        /* Station cells are ocean by definition in the original project.
           Prefer sampling the station cell itself whenever the solver mask
           leaves it open, even if DEP came from a heightmap color that says 0. */
        if (!g->solid[(size_t)station_idx]) {
            best_idx = station_idx;
            best_d2 = 0;
        }

        for (int r = 0; r <= max_radius && best_d2 != 0; ++r) {
            int x0 = sx - r;
            int x1 = sx + r;
            int y0 = sy - r;
            int y1 = sy + r;
            if (x0 < 0) x0 = 0;
            if (y0 < 0) y0 = 0;
            if (x1 >= (int)W) x1 = (int)W - 1;
            if (y1 >= (int)H) y1 = (int)H - 1;

            for (int y = y0; y <= y1; ++y) {
                for (int x = x0; x <= x1; ++x) {
                    if (x != x0 && x != x1 && y != y0 && y != y1 && r > 0) continue;
                    size_t idx = (size_t)y * W + (uint32_t)x;
                    if (g->solid[idx]) continue;
                    if (g->land[idx]) continue;
                    int dx = x - sx;
                    int dy = y - sy;
                    int d2 = dx * dx + dy * dy;
                    if (d2 < best_d2) {
                        best_d2 = d2;
                        best_idx = (int)idx;
                    }
                }
            }
            if (best_idx >= 0) break;
        }

        if (best_idx >= 0) g->station_sample_index[i] = best_idx;
    }
}


static float east_gate_reference_depth_for_row(const GridData *g, int row_y) {
    if (!g || !g->depth_m || !g->solid) return OUTER_ABSORB_DEPTH_M;

    const int W = (int)g->width;
    const int H = (int)g->height;
    if (W <= 0 || H <= 0) return OUTER_ABSORB_DEPTH_M;

    int gate_start_x = W - EAST_VISIBLE_GATE_CELLS;
    if (gate_start_x < 0) gate_start_x = 0;

    /*
       The hidden absorber should connect to the NORMAL offshore depth, not to
       the green/shallow DEP pixels left by the removed artificial station line.
       Scan westward from the gate and average the first usable open-water DEP
       cells from nearby rows.  The 4000 m floor intentionally keeps this east
       open boundary deep enough to avoid a shallow artificial shelf/reflector.
    */
    float sum = 0.0f;
    int count = 0;
    const int max_count = 48;

    for (int row_radius = 0; row_radius <= 3 && count < max_count; ++row_radius) {
        int y0 = row_y - row_radius;
        int y1 = row_y + row_radius;
        if (y0 < 0) y0 = 0;
        if (y1 >= H) y1 = H - 1;

        for (int y = y0; y <= y1 && count < max_count; ++y) {
            for (int x = gate_start_x - 1; x >= 0 && x >= gate_start_x - EAST_GATE_DEPTH_SCAN_CELLS && count < max_count; --x) {
                size_t idx = (size_t)y * g->width + (uint32_t)x;
                if (g->solid[idx]) continue;
                if (g->edge_absorb_gate && g->edge_absorb_gate[idx]) continue;

                float d = g->depth_m[idx];
                if (d <= 0.0f) continue;
                sum += clampf_local(d, 1.0f, 7000.0f);
                count++;
            }
        }
    }

    if (count > 0) {
        float avg = sum / (float)count;
        return clampf_local(avg, EAST_GATE_MIN_DEPTH_M, 7000.0f);
    }

    return OUTER_ABSORB_DEPTH_M;
}


static int is_force_water_repair_cell(const GridData *g, int visible_idx) {
    if (!g || !g->user_override || visible_idx < 0) return 0;
    size_t cells = (size_t)g->width * (size_t)g->height;
    if ((size_t)visible_idx >= cells) return 0;
    return g->user_override[(size_t)visible_idx] == 0;
}

static float reference_depth_for_force_water_cell(const GridData *g, int fx, int fy) {
    if (!g || !g->depth_m || !g->land || !g->solid) return FORCE_WATER_MIN_DEPTH_M;

    const int W = (int)g->width;
    const int H = (int)g->height;
    if (fx < 0 || fy < 0 || fx >= W || fy >= H) return FORCE_WATER_MIN_DEPTH_M;

    /*
       Manual force-water pixels usually repair DEP/heightmap mistakes such as
       a narrow strait, bay mouth, or false island.  Their original DEP value is
       often 0 or very shallow, which made the corrected water behave like a
       1 m-deep shelf.  Copy a nearby NORMAL ocean depth instead.

       Exclude other force-water cells from the reference set so a whole drawn
       channel does not keep copying its own bad 0-depth pixels.  Use the first
       usable ring and average several cells to avoid a noisy one-pixel value.
    */
    const int max_count = 32;
    for (int r = 1; r <= FORCE_WATER_DEPTH_SCAN_CELLS; ++r) {
        float sum = 0.0f;
        int count = 0;
        int x0 = fx - r;
        int x1 = fx + r;
        int y0 = fy - r;
        int y1 = fy + r;
        if (x0 < 0) x0 = 0;
        if (y0 < 0) y0 = 0;
        if (x1 >= W) x1 = W - 1;
        if (y1 >= H) y1 = H - 1;

        for (int y = y0; y <= y1 && count < max_count; ++y) {
            for (int x = x0; x <= x1 && count < max_count; ++x) {
                if (x != x0 && x != x1 && y != y0 && y != y1) continue;
                size_t idx = (size_t)y * g->width + (uint32_t)x;
                if (g->land[idx] || g->solid[idx]) continue;
                if (g->user_override && g->user_override[idx] == 0) continue;
                if (g->edge_absorb_gate && g->edge_absorb_gate[idx]) continue;
                float d = g->depth_m[idx];
                if (d <= 0.0f) continue;
                sum += clampf_local(d, FORCE_WATER_MIN_DEPTH_M, 7000.0f);
                count++;
            }
        }
        if (count > 0) {
            return clampf_local(sum / (float)count, FORCE_WATER_MIN_DEPTH_M, 7000.0f);
        }
    }

    return FORCE_WATER_MIN_DEPTH_M;
}

static void rebuild_effective_depths(GridData *g) {
    if (!g || !g->effective_depth_m || !g->depth_m) return;

    const size_t cells = (size_t)g->width * (size_t)g->height;
    for (size_t i = 0; i < cells; ++i) {
        g->effective_depth_m[i] = clampf_local(g->depth_m[i], 1.0f, 7000.0f);
    }

    if (!g->user_override) return;

    for (uint32_t y = 0; y < g->height; ++y) {
        for (uint32_t x = 0; x < g->width; ++x) {
            size_t idx = (size_t)y * g->width + x;
            if (g->user_override[idx] == 0) {
                g->effective_depth_m[idx] = reference_depth_for_force_water_cell(g, (int)x, (int)y);
            }
        }
    }
}

float visible_effective_depth_m(const GridData *g, int visible_idx) {
    if (!g || visible_idx < 0) return 1.0f;
    size_t cells = (size_t)g->width * (size_t)g->height;
    if ((size_t)visible_idx >= cells) return 1.0f;
    if (g->effective_depth_m) return clampf_local(g->effective_depth_m[(size_t)visible_idx], 1.0f, 7000.0f);
    return clampf_local(g->depth_m[(size_t)visible_idx], 1.0f, 7000.0f);
}

static void rebuild_edge_absorb_depths(GridData *g) {
    if (!g || !g->edge_absorb_depth_m || !g->depth_m) return;

    const size_t cells = (size_t)g->width * (size_t)g->height;
    for (size_t i = 0; i < cells; ++i) {
        g->edge_absorb_depth_m[i] = visible_effective_depth_m(g, (int)i);
    }

    if (!g->edge_absorb_gate) return;

    int gate_start_x = (int)g->width - EAST_VISIBLE_GATE_CELLS;
    if (gate_start_x < 0) gate_start_x = 0;

    for (uint32_t y = 0; y < g->height; ++y) {
        float ref_depth = east_gate_reference_depth_for_row(g, (int)y);
        for (int x = gate_start_x; x < (int)g->width; ++x) {
            size_t idx = (size_t)y * g->width + (uint32_t)x;
            if (g->edge_absorb_gate[idx]) {
                g->edge_absorb_depth_m[idx] = ref_depth;
            }
        }
    }
}

void rebuild_solid_mask(GridData *g, int coast_guard_cells) {
    if (!g || !g->land || !g->solid) return;
    if (coast_guard_cells < 0) coast_guard_cells = 0;
    if (coast_guard_cells > 8) coast_guard_cells = 8;

    size_t cells = (size_t)g->width * (size_t)g->height;

    /*
       Make the wave field match the displayed land/sea boundary.

       Older builds expanded each land cell by coast_guard_cells into the ocean
       side.  That was useful for preventing one-pixel leaks, but it also moved
       the solver wall one or more cells offshore while the map still showed
       those cells as sea.  The result was the observed mismatch: a strong wave
       appeared one cell offshore, while the displayed coastline/observation
       cell stayed weak or unreachable.

       Now solid is primarily the same as land.  Narrow-channel/leak control is
       handled by the existing coastal absorption / strait transmission logic
       rather than by moving the solver coastline offshore.
    */
    for (size_t i = 0; i < cells; ++i) {
        g->solid[i] = g->land[i] ? 1u : 0u;
    }

    /* Manual overrides must affect the solver boundary exactly as displayed. */
    if (g->user_override) {
        for (size_t i = 0; i < cells; ++i) {
            if (g->user_override[i] == 0) g->solid[i] = 0u;
            else if (g->user_override[i] > 0) g->solid[i] = 1u;
        }
    }

    /*
       Observation points are ocean-side gauges in the original project and
       should remain reachable.  Reopen station cells and their immediate
       water-classified neighbors so the wave field reaches the displayed
       coastline/marker instead of stopping offshore.
    */
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        int32_t station_idx = g->station_index[i];
        if (station_idx < 0 || (size_t)station_idx >= cells) continue;

        int sx = station_idx % (int)g->width;
        int sy = station_idx / (int)g->width;

        for (int yy = sy - 1; yy <= sy + 1; ++yy) {
            for (int xx = sx - 1; xx <= sx + 1; ++xx) {
                if (xx < 0 || yy < 0 || xx >= (int)g->width || yy >= (int)g->height) continue;
                size_t idx = (size_t)yy * g->width + (uint32_t)xx;
                if (!g->land[idx]) g->solid[idx] = 0u;
            }
        }
    }

    if (g->edge_absorb_gate) {
        for (size_t i = 0; i < cells; ++i) {
            if (g->edge_absorb_gate[i] && (!g->user_override || g->user_override[i] <= 0)) {
                g->solid[i] = 0u;
            }
        }
    }

    rebuild_effective_depths(g);
    rebuild_station_sample_indices(g);
}


void clear_land_overrides(GridData *g) {
    if (!g || !g->user_override) return;
    size_t cells = (size_t)g->width * (size_t)g->height;
    for (size_t i = 0; i < cells; ++i) g->user_override[i] = -1;
}

void set_land_override(GridData *g, int grid_x, int grid_y, int mode) {
    if (!g || !g->user_override) return;
    if (grid_x < 0 || grid_y < 0 || grid_x >= (int)g->width || grid_y >= (int)g->height) return;
    size_t idx = tsunami_index(g, grid_x, grid_y);
    if (mode > 0) g->user_override[idx] = 1;
    else if (mode == 0) g->user_override[idx] = 0;
    else g->user_override[idx] = -1;
}

size_t count_land_overrides(const GridData *g) {
    if (!g || !g->user_override) return 0;
    size_t cells = (size_t)g->width * (size_t)g->height;
    size_t count = 0;
    for (size_t i = 0; i < cells; ++i) {
        if (g->user_override[i] >= 0) count++;
    }
    return count;
}

int save_land_override_csv(const char *path, const GridData *g) {
    if (!g || !g->user_override) return 0;
    FILE *fp = fopen(path, "wb");
    if (!fp) return 0;
    fputs("grid_x,grid_y,mode\n", fp);
    for (uint32_t y = 0; y < g->height; ++y) {
        for (uint32_t x = 0; x < g->width; ++x) {
            size_t idx = (size_t)y * g->width + x;
            int v = g->user_override[idx];
            if (v < 0) continue;
            fprintf(fp, "%u,%u,%s\n", x, y, v ? "land" : "water");
        }
    }
    fclose(fp);
    return 1;
}

static int parse_land_override_line(const char *text, int *x, int *y, int *mode_out) {
    char mode_text[64] = {0};
    char extra = '\0';
    if (sscanf(text, " %d , %d , %63[^,\r\n ] %c", x, y, mode_text, &extra) < 3) return 0;
    if (strcmp(mode_text, "1") == 0 || strcmp(mode_text, "land") == 0 || strcmp(mode_text, "LAND") == 0 || strcmp(mode_text, "l") == 0 || strcmp(mode_text, "L") == 0) {
        *mode_out = 1;
        return 1;
    }
    if (strcmp(mode_text, "0") == 0 || strcmp(mode_text, "water") == 0 || strcmp(mode_text, "WATER") == 0 || strcmp(mode_text, "w") == 0 || strcmp(mode_text, "W") == 0) {
        *mode_out = 0;
        return 1;
    }
    if (strcmp(mode_text, "-1") == 0 || strcmp(mode_text, "clear") == 0 || strcmp(mode_text, "CLEAR") == 0 || strcmp(mode_text, "none") == 0) {
        *mode_out = -1;
        return 1;
    }
    return 0;
}

int load_land_override_csv(GridData *g, const char *path) {
    if (!g || !g->user_override) return 0;
    FILE *fp = fopen(path, "r");
    if (!fp) return 0;
    clear_land_overrides(g);
    char line[256];
    int line_no = 0;
    while (fgets(line, sizeof(line), fp)) {
        line_no++;
        char *p = line;
        while (*p == ' ' || *p == '\t') p++;
        if (*p == '\0' || *p == '\n' || *p == '\r' || *p == '#') continue;
        if (line_no == 1 && strstr(p, "grid_x") != NULL) continue;
        int x = 0, y = 0, mode = -1;
        if (parse_land_override_line(p, &x, &y, &mode)) {
            set_land_override(g, x, y, mode);
        }
    }
    fclose(fp);
    rebuild_land_mask_from_coast_stations(g);
    return 1;
}

void free_grid_data(GridData *g) {
    free(g->depth_m);
    free(g->effective_depth_m);
    free(g->base_land);
    free(g->land);
    free(g->solid);
    free(g->user_override);
    free(g->station_index);
    free(g->station_sample_index);
    free(g->station_id);
    free(g->station_enabled);
    free(g->edge_absorb_gate);
    free(g->edge_absorb_depth_m);
    memset(g, 0, sizeof(*g));
}

static int clampi_core(int x, int lo, int hi) {
    if (x < lo) return lo;
    if (x > hi) return hi;
    return x;
}

static size_t solver_index_core(const Simulation *s, int sx, int sy) {
    return (size_t)sy * (size_t)s->solver_width + (uint32_t)sx;
}

static int solver_to_visible_index(const GridData *g, const Simulation *s, int sx, int sy) {
    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;
    if (vx < 0 || vy < 0 || vx >= (int)g->width || vy >= (int)g->height) return -1;
    return (int)((size_t)vy * g->width + (uint32_t)vx);
}

static int nearest_visible_index_for_solver(const GridData *g, const Simulation *s, int sx, int sy) {
    int vx = clampi_core(sx - (int)s->outer_pad, 0, (int)g->width - 1);
    int vy = clampi_core(sy - (int)s->outer_pad, 0, (int)g->height - 1);
    return (int)((size_t)vy * g->width + (uint32_t)vx);
}

static int visible_cell_is_coastal_water(const GridData *g, const Simulation *s, int sx, int sy);
static int visible_cell_is_one_cell_narrow_sea(const GridData *g, const Simulation *s, int sx, int sy);
static int visible_cell_is_narrow_reflection_bypass_zone(const GridData *g, const Simulation *s, int sx, int sy);
static int visible_cell_is_strait_transmission_zone(const GridData *g, const Simulation *s, int sx, int sy);
static int diagonal_move_crosses_strait_transmission_zone(const GridData *g, const Simulation *s,
                                                         int center_sx, int center_sy, int dx, int dy);


static int is_east_absorber_gate_solver_cell(const GridData *g, const Simulation *s, int sx, int sy) {
    if (!g || !s || !g->edge_absorb_gate) return 0;
    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;
    if (vx < (int)g->width) return 0;          /* only the east outside pad */
    if (vy < 0 || vy >= (int)g->height) return 0;
    int edge_idx = (int)((size_t)vy * g->width + (g->width - 1u));
    return g->edge_absorb_gate[(size_t)edge_idx] != 0;
}

static int solver_solid_at_or_outside(const GridData *g, const Simulation *s, int sx, int sy) {
    if (sx < 0 || sy < 0 || sx >= (int)s->solver_width || sy >= (int)s->solver_height) return 1;

    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx >= 0) {
        if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx]) return 0;
        return g->solid[(size_t)visible_idx] != 0;
    }

    if (is_east_absorber_gate_solver_cell(g, s, sx, sy)) return 0;

    /* Hidden outside absorber cells do not have their own coastline data.
       Use the nearest visible edge cell as a gate: if that edge cell is land,
       this outside cell is also treated as solid. That prevents waves from
       escaping through land edges and then sneaking around behind the map. */
    int nearest_idx = nearest_visible_index_for_solver(g, s, sx, sy);
    return g->solid[(size_t)nearest_idx] != 0;
}

static float solver_depth_m(const GridData *g, const Simulation *s, int sx, int sy) {
    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx >= 0) {
        if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx] && g->edge_absorb_depth_m) {
            return clampf_local(g->edge_absorb_depth_m[(size_t)visible_idx], EAST_GATE_MIN_DEPTH_M, 7000.0f);
        }
        return visible_effective_depth_m(g, visible_idx);
    }

    if (is_east_absorber_gate_solver_cell(g, s, sx, sy)) {
        int vx = sx - (int)s->outer_pad;
        int vy = sy - (int)s->outer_pad;
        (void)vx;
        int edge_idx = (int)((size_t)vy * g->width + (g->width - 1u));
        if (g->edge_absorb_depth_m) {
            return clampf_local(g->edge_absorb_depth_m[(size_t)edge_idx], EAST_GATE_MIN_DEPTH_M, 7000.0f);
        }
        return east_gate_reference_depth_for_row(g, vy);
    }

    int nearest_idx = nearest_visible_index_for_solver(g, s, sx, sy);
    if (!g->solid[(size_t)nearest_idx]) {
        return visible_effective_depth_m(g, nearest_idx);
    }
    return OUTER_ABSORB_DEPTH_M;
}

Simulation create_simulation(const GridData *g) {
    Simulation s;
    memset(&s, 0, sizeof(s));
    s.width = g->width;
    s.height = g->height;
    s.padded_height = g->height + 2;
    s.outer_pad = OUTER_ABSORB_CELLS;
    s.solver_width = g->width + 2u * s.outer_pad;
    s.solver_height = g->height + 2u * s.outer_pad;
    s.cells = (size_t)g->width * (size_t)g->height;
    s.padded_cells = (size_t)g->width * (size_t)(g->height + 2);
    s.solver_cells = (size_t)s.solver_width * (size_t)s.solver_height;

    s.height_norm = (float *)tsunami_xcalloc(s.solver_cells, sizeof(float));
    s.height_prev = (float *)tsunami_xcalloc(s.solver_cells, sizeof(float));
    s.velocity = (float *)tsunami_xcalloc(s.solver_cells, sizeof(float));
    s.display_height_m = (float *)tsunami_xcalloc(s.padded_cells, sizeof(float));
    s.display_max_m = (float *)tsunami_xcalloc(s.padded_cells, sizeof(float));
    s.display_arrival_s = (int32_t *)tsunami_xcalloc(s.padded_cells, sizeof(int32_t));
    s.source_height_m = (float *)tsunami_xcalloc(s.cells, sizeof(float));
    s.source_delay_s = (int32_t *)tsunami_xcalloc(s.cells, sizeof(int32_t));
    s.source_applied_m = (float *)tsunami_xcalloc(s.cells, sizeof(float));
    s.station_max_m = (float *)tsunami_xcalloc(g->station_count, sizeof(float));
    s.station_arrival_s = (int32_t *)tsunami_xcalloc(g->station_count, sizeof(int32_t));
    reset_simulation(g, &s);
    return s;
}

void reset_simulation(const GridData *g, Simulation *s) {
    memset(s->height_norm, 0, s->solver_cells * sizeof(float));
    memset(s->height_prev, 0, s->solver_cells * sizeof(float));
    memset(s->velocity, 0, s->solver_cells * sizeof(float));
    memset(s->display_height_m, 0, s->padded_cells * sizeof(float));
    memset(s->display_max_m, 0, s->padded_cells * sizeof(float));
    memset(s->display_arrival_s, 0, s->padded_cells * sizeof(int32_t));
    memset(s->source_height_m, 0, s->cells * sizeof(float));
    memset(s->source_delay_s, 0, s->cells * sizeof(int32_t));
    memset(s->source_applied_m, 0, s->cells * sizeof(float));
    memset(s->station_max_m, 0, g->station_count * sizeof(float));
    memset(s->station_arrival_s, 0, g->station_count * sizeof(int32_t));

    for (uint32_t sy = 0; sy < s->solver_height; ++sy) {
        for (uint32_t sx = 0; sx < s->solver_width; ++sx) {
            size_t p = solver_index_core(s, (int)sx, (int)sy);
            if (solver_solid_at_or_outside(g, s, (int)sx, (int)sy)) {
                s->height_norm[p] = SCRATCH_LAND_VALUE;
            }
        }
    }
}

void free_simulation(Simulation *s) {
    free(s->height_norm);
    free(s->height_prev);
    free(s->velocity);
    free(s->display_height_m);
    free(s->display_max_m);
    free(s->display_arrival_s);
    free(s->source_height_m);
    free(s->source_delay_s);
    free(s->source_applied_m);
    free(s->station_max_m);
    free(s->station_arrival_s);
    memset(s, 0, sizeof(*s));
}

void add_source_grid(const GridData *g, Simulation *s, int grid_x, int grid_y, float height_m, int delay_s) {
    if (grid_x < 0 || grid_y < 0 || grid_x >= (int)g->width || grid_y >= (int)g->height) return;
    size_t idx = (size_t)grid_y * g->width + (uint32_t)grid_x;
    if (g->solid[idx]) return;
    if (s->source_height_m[idx] != 0.0f && s->source_delay_s[idx] != delay_s) {
        s->source_height_m[idx] = height_m;
    } else {
        s->source_height_m[idx] += height_m;
    }
    s->source_delay_s[idx] = delay_s;
    if (s->source_applied_m) s->source_applied_m[idx] = 0.0f;
}

/*
   Return a neighbor's old water height.  For land neighbors, return either the
   current value (hard wall) or a partially absorbing ghost value when coastal
   reflection absorption is enabled.

   Important detail: diagonal neighbors are also blocked when either of the two
   side-adjacent cells is land. Without this "no corner cutting" rule, a wave
   can leak through a one-pixel coastline corner, which looks like it is going
   through land in the GUI.
*/
static float coastal_reflection_keep_for_cell(const GridData *g,
                                                   const Simulation *s,
                                                   int center_sx,
                                                   int center_sy) {
    if (!coastal_reflection_absorption_enabled()) return 1.0f;

    int visible_idx = solver_to_visible_index(g, s, center_sx, center_sy);
    if (visible_idx < 0 || (size_t)visible_idx >= g->width * g->height) {
        return coastal_reflection_keep_per_step();
    }

    /*
       Narrow straits and manually forced-water repair cells are usually real
       water passages that were under-resolved by the color heightmap.  Do not
       apply the coastal partial-reflection/absorption boundary in those cells
       or their immediate gate neighborhood; otherwise every side wall in a
       one-cell passage removes energy and the wave height becomes too low.

       This bypass only affects the coastal partial-reflection boundary.  The
       optional one-cell narrow-sea damping (J / 1 / 2 keys) remains separate.
    */
    if (visible_cell_is_strait_transmission_zone(g, s, center_sx, center_sy)) {
        return 1.0f;
    }

    /*
       v129: remove shallow-water/coastal-depth damping.

       v126 applied additional partial-reflection loss only in water shallower
       than 120 m.  That made ports, bays, and shallow coastal cells lose energy
       before the user could judge whether the scalar source itself was correct.
       Keep the coast-absorb UI/setting for compatibility, but do not apply any
       depth-based shallow-water loss here.  Coastline cells therefore behave as
       hard-wall reflection unless another explicit absorber path handles them.
    */
    return 1.0f;
}

static float coastal_absorbing_ghost_value(const GridData *g,
                                          const Simulation *s,
                                          int center_sx,
                                          int center_sy,
                                          size_t center_p) {
    float center = s->height_prev[center_p];

    /*
       Land neighbors use a ghost value.  Ordinary open coastlines use a
       partially absorbing boundary, but narrow/very shallow passages bypass
       or weaken that absorption so the wave does not become too low while
       passing through an under-resolved strait or bay mouth.
    */
    if (coastal_reflection_absorption_enabled() &&
        solver_to_visible_index(g, s, center_sx, center_sy) >= 0 &&
        visible_cell_is_coastal_water(g, s, center_sx, center_sy)) {
        return center * coastal_reflection_keep_for_cell(g, s, center_sx, center_sy);
    }
    return center;
}

static float neighbor_or_self_blocked(const GridData *g,
                                      const Simulation *s,
                                      int center_sx,
                                      int center_sy,
                                      int dx,
                                      int dy,
                                      size_t center_p) {
    int nx = center_sx + dx;
    int ny = center_sy + dy;

    if (solver_solid_at_or_outside(g, s, nx, ny)) {
        return coastal_absorbing_ghost_value(g, s, center_sx, center_sy, center_p);
    }

    if (dx != 0 && dy != 0) {
        if (solver_solid_at_or_outside(g, s, center_sx + dx, center_sy) ||
            solver_solid_at_or_outside(g, s, center_sx, center_sy + dy)) {
            if (!diagonal_move_crosses_strait_transmission_zone(g, s, center_sx, center_sy, dx, dy)) {
                return coastal_absorbing_ghost_value(g, s, center_sx, center_sy, center_p);
            }
        }
    }

    size_t neighbor_p = solver_index_core(s, nx, ny);
    float v = s->height_prev[neighbor_p];
    if (is_land_value(v)) return coastal_absorbing_ghost_value(g, s, center_sx, center_sy, center_p);
    return v;
}

static int neighbor_open_for_wave_flux(const GridData *g,
                                       const Simulation *s,
                                       int center_sx,
                                       int center_sy,
                                       int dx,
                                       int dy) {
    int nx = center_sx + dx;
    int ny = center_sy + dy;

    if (solver_solid_at_or_outside(g, s, nx, ny)) return 0;

    if (dx != 0 && dy != 0) {
        if (solver_solid_at_or_outside(g, s, center_sx + dx, center_sy) ||
            solver_solid_at_or_outside(g, s, center_sx, center_sy + dy)) {
            if (!diagonal_move_crosses_strait_transmission_zone(g, s, center_sx, center_sy, dx, dy)) {
                return 0;
            }
        }
    }

    {
        size_t neighbor_p = solver_index_core(s, nx, ny);
        if (is_land_value(s->height_prev[neighbor_p])) return 0;
    }
    return 1;
}

static float wave_flux_term_to_neighbor(const GridData *g,
                                        const Simulation *s,
                                        int center_sx,
                                        int center_sy,
                                        int dx,
                                        int dy,
                                        float center_h,
                                        float neighbor_h,
                                        float center_depth_m,
                                        float weight) {
#if WAVE_FACE_DEPTH_FLUX_ENABLED
    if (!neighbor_open_for_wave_flux(g, s, center_sx, center_sy, dx, dy)) {
        return (neighbor_h - center_h) * weight * wave_velocity_gain_from_depth(center_depth_m);
    }
    {
        float neighbor_depth_m = solver_depth_m(g, s, center_sx + dx, center_sy + dy);
        float face_gain = wave_velocity_gain_from_face_depth(center_depth_m, neighbor_depth_m);
        return (neighbor_h - center_h) * weight * face_gain;
    }
#else
    return (neighbor_h - center_h) * weight * wave_velocity_gain_from_depth(center_depth_m);
#endif
}

static float visible_east_gate_absorber_factor(const GridData *g, const Simulation *s, int sx, int sy) {
    if (!g || !s || !g->edge_absorb_gate) return 1.0f;

    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx < 0) return 1.0f;
    if (!g->edge_absorb_gate[(size_t)visible_idx]) return 1.0f;

    int vx = sx - (int)s->outer_pad;
    int distance_from_east_edge = ((int)s->width - 1) - vx;
    if (distance_from_east_edge < 0 || distance_from_east_edge >= EAST_VISIBLE_GATE_CELLS) return 1.0f;

    float t = (float)(EAST_VISIBLE_GATE_CELLS - distance_from_east_edge) / (float)EAST_VISIBLE_GATE_CELLS;
    float factor = 1.0f - ABSORB_MAX_STRENGTH * t * t;
    if (factor < 0.35f) factor = 0.35f;
    if (factor > 1.0f) factor = 1.0f;
    return factor;
}

static float hidden_outer_absorber_factor(const GridData *g, const Simulation *s, int sx, int sy) {
    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;

    int outside_x = 0;
    int outside_y = 0;
    if (vx < 0) outside_x = -vx;
    else if (vx >= (int)s->width) outside_x = vx - (int)s->width + 1;
    if (vy < 0) outside_y = -vy;
    else if (vy >= (int)s->height) outside_y = vy - (int)s->height + 1;

    int outside = outside_x > outside_y ? outside_x : outside_y;
    if (outside <= 0) return 1.0f;
    if (outside > (int)s->outer_pad) outside = (int)s->outer_pad;

    float t = (float)outside / (float)s->outer_pad;
    float factor = 1.0f - ABSORB_MAX_STRENGTH * t * t;
    if (factor < 0.35f) factor = 0.35f;
    if (factor > 1.0f) factor = 1.0f;
    return factor;
}

static int visible_cell_is_coastal_water(const GridData *g, const Simulation *s, int sx, int sy) {
    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx < 0) return 0;
    if ((size_t)visible_idx >= g->width * g->height) return 0;
    if (g->solid[(size_t)visible_idx]) return 0;
    if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx]) return 0;
    if (is_force_water_repair_cell(g, visible_idx)) return 0;

    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;
    const int dx[4] = { 1, -1, 0, 0 };
    const int dy[4] = { 0, 0, 1, -1 };
    for (int k = 0; k < 4; ++k) {
        int nx = vx + dx[k];
        int ny = vy + dy[k];
        if (nx < 0 || ny < 0 || nx >= (int)g->width || ny >= (int)g->height) continue;
        size_t nidx = (size_t)ny * g->width + (uint32_t)nx;
        if (g->solid[nidx]) return 1;
    }
    return 0;
}

static int visible_cell_is_one_cell_narrow_sea(const GridData *g, const Simulation *s, int sx, int sy) {
    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx < 0) return 0;
    if ((size_t)visible_idx >= g->width * g->height) return 0;
    if (g->solid[(size_t)visible_idx]) return 0;
    if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx]) return 0;

    /* Manually forced-water cells are usually the places where the color
       heightmap closed a real narrow sea.  Treat them as one-cell gates. */
    if (is_force_water_repair_cell(g, visible_idx)) return 1;

    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;

    int land_w = (vx <= 0) ? 1 : g->solid[(size_t)vy * g->width + (uint32_t)(vx - 1)] != 0;
    int land_e = (vx >= (int)g->width - 1) ? 1 : g->solid[(size_t)vy * g->width + (uint32_t)(vx + 1)] != 0;
    int land_n = (vy <= 0) ? 1 : g->solid[(size_t)(vy - 1) * g->width + (uint32_t)vx] != 0;
    int land_s = (vy >= (int)g->height - 1) ? 1 : g->solid[(size_t)(vy + 1) * g->width + (uint32_t)vx] != 0;

    /* A one-cell-wide water corridor has land on both opposite sides.
       This catches thin bays/straits without damping ordinary open coastlines. */
    if ((land_n && land_s) || (land_w && land_e)) return 1;
    return 0;
}

static int visible_cell_is_narrow_reflection_bypass_zone(const GridData *g, const Simulation *s, int sx, int sy) {
    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx < 0) return 0;
    if ((size_t)visible_idx >= g->width * g->height) return 0;
    if (g->solid[(size_t)visible_idx]) return 0;
    if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx]) return 0;

    /* Direct one-cell gate / forced-water repair cell. */
    if (visible_cell_is_one_cell_narrow_sea(g, s, sx, sy)) return 1;

    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;

    /* Two or more land sides usually means a corner, bay throat, or narrow
       inner sea at this map resolution.  Do not partially absorb those walls. */
    int land4 = 0;
    const int dx4[4] = { 1, -1, 0, 0 };
    const int dy4[4] = { 0, 0, 1, -1 };
    for (int k = 0; k < 4; ++k) {
        int nx = vx + dx4[k];
        int ny = vy + dy4[k];
        if (nx < 0 || ny < 0 || nx >= (int)g->width || ny >= (int)g->height) {
            land4++;
            continue;
        }
        if (g->solid[(size_t)ny * g->width + (uint32_t)nx]) land4++;
    }
    if (land4 >= 2) return 1;

    /* Cells immediately touching a manually opened gate should also bypass
       partial reflection, otherwise the gate itself is open but its side cells
       still absorb too much energy. */
    for (int oy = -1; oy <= 1; ++oy) {
        for (int ox = -1; ox <= 1; ++ox) {
            if (ox == 0 && oy == 0) continue;
            int nx = vx + ox;
            int ny = vy + oy;
            if (nx < 0 || ny < 0 || nx >= (int)g->width || ny >= (int)g->height) continue;
            int nidx = (int)((size_t)ny * g->width + (uint32_t)nx);
            if (is_force_water_repair_cell(g, nidx)) return 1;
        }
    }

    return 0;
}

static int visible_cell_water_open_for_strait(const GridData *g, int vx, int vy) {
    if (!g || vx < 0 || vy < 0 || vx >= (int)g->width || vy >= (int)g->height) return 0;
    size_t idx = (size_t)vy * g->width + (uint32_t)vx;
    if (g->edge_absorb_gate && g->edge_absorb_gate[idx]) return 1;
    return g->solid[idx] == 0;
}

static int visible_cell_solid_for_strait(const GridData *g, int vx, int vy) {
    if (!g || vx < 0 || vy < 0 || vx >= (int)g->width || vy >= (int)g->height) return 1;
    size_t idx = (size_t)vy * g->width + (uint32_t)vx;
    if (g->edge_absorb_gate && g->edge_absorb_gate[idx]) return 0;
    return g->solid[idx] != 0;
}

static int visible_cell_is_strait_transmission_zone(const GridData *g, const Simulation *s, int sx, int sy) {
    /*
       Broader strait/bay-mouth transmission zone.

       Tokyo Bay entrance and the Seto Inland Sea are not always one-cell gates;
       they are often several under-resolved coastal cells surrounded by nearby
       headlands/islands.  If ordinary coastal absorption is applied on every
       sidewall, both incoming and outgoing waves become too small.  Treat cells
       inside a nearby constriction as transmission cells:
       - bypass coastal partial-absorption,
       - allow diagonal connection through under-resolved corners.
    */
    int visible_idx = solver_to_visible_index(g, s, sx, sy);
    if (visible_idx < 0) return 0;
    if ((size_t)visible_idx >= g->width * g->height) return 0;
    if (g->solid[(size_t)visible_idx]) return 0;
    if (g->edge_absorb_gate && g->edge_absorb_gate[(size_t)visible_idx]) return 0;

    if (visible_cell_is_narrow_reflection_bypass_zone(g, s, sx, sy)) return 1;

    int vx = sx - (int)s->outer_pad;
    int vy = sy - (int)s->outer_pad;

    /*
       Look within two cells.  A strait/bay mouth often has land on both sides
       in one axis but open water through the other axis.  Use radius 2 so
       two- to five-cell-wide entrances are not treated like ordinary coastline.
    */
    int solid_w = 0, solid_e = 0, solid_n = 0, solid_s = 0;
    int open_w = 0, open_e = 0, open_n = 0, open_s = 0;

    for (int d = 1; d <= 2; ++d) {
        if (visible_cell_solid_for_strait(g, vx - d, vy)) solid_w++;
        if (visible_cell_solid_for_strait(g, vx + d, vy)) solid_e++;
        if (visible_cell_solid_for_strait(g, vx, vy - d)) solid_n++;
        if (visible_cell_solid_for_strait(g, vx, vy + d)) solid_s++;

        if (visible_cell_water_open_for_strait(g, vx - d, vy)) open_w++;
        if (visible_cell_water_open_for_strait(g, vx + d, vy)) open_e++;
        if (visible_cell_water_open_for_strait(g, vx, vy - d)) open_n++;
        if (visible_cell_water_open_for_strait(g, vx, vy + d)) open_s++;
    }

    int horizontal_gate = (solid_n > 0 && solid_s > 0 && open_w > 0 && open_e > 0);
    int vertical_gate   = (solid_w > 0 && solid_e > 0 && open_n > 0 && open_s > 0);
    if (horizontal_gate || vertical_gate) return 1;

    /*
       Dense island/coast maze: many land cells nearby but still open water on
       at least two sides.  This catches Seto Inland Sea-style island passages
       without opening ordinary straight coastlines too much.
    */
    int nearby_solid = 0;
    int nearby_water = 0;
    for (int oy = -2; oy <= 2; ++oy) {
        for (int ox = -2; ox <= 2; ++ox) {
            if (ox == 0 && oy == 0) continue;
            int nx = vx + ox;
            int ny = vy + oy;
            if (visible_cell_solid_for_strait(g, nx, ny)) nearby_solid++;
            else if (visible_cell_water_open_for_strait(g, nx, ny)) nearby_water++;
        }
    }

    int open_sides = 0;
    if (open_w > 0) open_sides++;
    if (open_e > 0) open_sides++;
    if (open_n > 0) open_sides++;
    if (open_s > 0) open_sides++;

    if (nearby_solid >= 5 && nearby_water >= 5 && open_sides >= 2) return 1;

    return 0;
}

static int diagonal_move_crosses_strait_transmission_zone(const GridData *g, const Simulation *s,
                                                         int center_sx, int center_sy, int dx, int dy) {
    if (dx == 0 || dy == 0) return 0;

    /*
       The ordinary no-corner-cutting rule is correct for clean coastlines, but
       it can close or over-throttle diagonal flow in pixelated bay mouths and
       island-strait networks.  If the center, target, or side-adjacent cells
       are inside a strait transmission zone, allow the diagonal connection.
    */
    if (visible_cell_is_strait_transmission_zone(g, s, center_sx, center_sy)) return 1;
    if (visible_cell_is_strait_transmission_zone(g, s, center_sx + dx, center_sy + dy)) return 1;
    if (visible_cell_is_strait_transmission_zone(g, s, center_sx + dx, center_sy)) return 1;
    if (visible_cell_is_strait_transmission_zone(g, s, center_sx, center_sy + dy)) return 1;
    return 0;
}

static float narrow_sea_damping_factor(const GridData *g, const Simulation *s, int sx, int sy) {
    if (!narrow_sea_damping_enabled()) return 1.0f;
    if (!visible_cell_is_one_cell_narrow_sea(g, s, sx, sy)) return 1.0f;
    return narrow_sea_keep_per_step();
}

static int coast_cell_returning_toward_zero(float current_h, float predicted_h) {
    if (current_h == 0.0f || predicted_h == 0.0f) return 0;
    if (current_h * predicted_h <= 0.0f) return 0;
    return fabsf(predicted_h) < fabsf(current_h);
}

static void simulate_one_numerical_substep(const GridData *g, Simulation *s, int time_s) {
    const uint32_t SW = s->solver_width;
    const uint32_t SH = s->solver_height;
    const float diag_weight = 0.7071067811865476f;

    memcpy(s->height_prev, s->height_norm, s->solver_cells * sizeof(float));

    for (uint32_t sy = 0; sy < SH; ++sy) {
        for (uint32_t sx = 0; sx < SW; ++sx) {
            size_t p = solver_index_core(s, (int)sx, (int)sy);
            int visible_idx = solver_to_visible_index(g, s, (int)sx, (int)sy);

            if (solver_solid_at_or_outside(g, s, (int)sx, (int)sy) || is_land_value(s->height_prev[p])) {
                s->height_norm[p] = SCRATCH_LAND_VALUE;
                s->velocity[p] = 0.0f;
                if (visible_idx >= 0) {
                    s->display_height_m[(size_t)visible_idx + g->width] = 0.0f;
                    s->display_max_m[(size_t)visible_idx + g->width] = 0.0f;
                    s->display_arrival_s[(size_t)visible_idx + g->width] = 0;
                }
                continue;
            }

            float center = s->height_prev[p];
            float nw = neighbor_or_self_blocked(g, s, (int)sx, (int)sy, -1, -1, p);
            float n  = neighbor_or_self_blocked(g, s, (int)sx, (int)sy,  0, -1, p);
            float ne = neighbor_or_self_blocked(g, s, (int)sx, (int)sy,  1, -1, p);
            float w  = neighbor_or_self_blocked(g, s, (int)sx, (int)sy, -1,  0, p);
            float e  = neighbor_or_self_blocked(g, s, (int)sx, (int)sy,  1,  0, p);
            float sw = neighbor_or_self_blocked(g, s, (int)sx, (int)sy, -1,  1, p);
            float s0 = neighbor_or_self_blocked(g, s, (int)sx, (int)sy,  0,  1, p);
            float se = neighbor_or_self_blocked(g, s, (int)sx, (int)sy,  1,  1, p);

            float depth = solver_depth_m(g, s, (int)sx, (int)sy);
            float wave_flux_accel =
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy, -1, -1, center, nw, depth, diag_weight) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy,  0, -1, center, n,  depth, 1.0f) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy,  1, -1, center, ne, depth, diag_weight) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy, -1,  0, center, w,  depth, 1.0f) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy,  1,  0, center, e,  depth, 1.0f) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy, -1,  1, center, sw, depth, diag_weight) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy,  0,  1, center, s0, depth, 1.0f) +
                wave_flux_term_to_neighbor(g, s, (int)sx, (int)sy,  1,  1, center, se, depth, diag_weight);

            s->velocity[p] += wave_flux_accel;

            /* Only the hidden outer pad is damped.  Do not apply any visible
               east-edge sponge; otherwise the east side differs from the
               south side and the last visible cells can lose wave height. */
            float damping = hidden_outer_absorber_factor(g, s, (int)sx, (int)sy);
            float global_decay = global_height_decay_per_step();
            float total_damping = damping * global_decay;
            s->velocity[p] *= total_damping;
            float next_h = (center + s->velocity[p]) * total_damping;

            /*
               Coastal reflection damping is handled in neighbor_or_self_blocked()
               through a partial-reflection ghost boundary.  Do not damp
               next_h here; damping the falling limb can make water stick to
               complex coastlines and even increase the apparent reflected wave.
            */

            if (visible_idx >= 0 &&
                s->source_height_m[(size_t)visible_idx] != 0.0f &&
                time_s != 0) {
                /*
                   Smooth rupture source injection.

                   Older versions injected the whole cell height only when
                   source_delay_s == time_s.  For giant ruptures whose delay
                   changes along a long trench, that produces visible delay
                   bands, like bamboo nodes.  Instead, each source cell ramps
                   in over g_source_ramp_seconds and only the newly reached
                   fraction is injected at this step.
                */
                size_t src_idx = (size_t)visible_idx;
                float delay = (float)s->source_delay_s[src_idx];
                float ramp = g_source_ramp_seconds;
                if (ramp < 1.0f) ramp = 1.0f;

                float u = ((float)time_s - delay) / ramp;
                float target_frac;
                if (u <= 0.0f) {
                    target_frac = 0.0f;
                } else if (u >= 1.0f) {
                    target_frac = 1.0f;
                } else {
                    target_frac = u * u * (3.0f - 2.0f * u);
                }

                float already = s->source_applied_m ? s->source_applied_m[src_idx] : 0.0f;
                float target_h = s->source_height_m[src_idx] * target_frac;
                float delta_h = target_h - already;

                if (fabsf(delta_h) > 1.0e-7f) {
                    next_h += delta_h / output_shoaling_factor(visible_effective_depth_m(g, visible_idx));
                    if (s->source_applied_m) s->source_applied_m[src_idx] = target_h;
                }
            }

            if (visible_idx >= 0 && visible_cell_is_shimonoseki_choke(g, visible_idx)) {
                float choke_keep = shimonoseki_choke_keep_per_internal_substep();
                s->velocity[p] *= choke_keep;
                next_h *= choke_keep;
            }

            s->height_norm[p] = next_h;
            if (visible_idx >= 0) {
                float visible_h =
                    next_h * output_shoaling_factor(visible_effective_depth_m(g, visible_idx));
                size_t display_idx = (size_t)visible_idx + g->width;
                s->display_height_m[display_idx] = visible_h;

                /*
                   Ocean-cell maximum wave-height layer for MOST-style review.
                   Track only the largest positive water-level rise at each
                   displayed ocean cell.  Negative drawdown/subsidence is not
                   counted as maximum wave height.
                */
                {
                    float positive_h = visible_h > 0.0f ? visible_h : 0.0f;
                    if (positive_h > s->display_max_m[display_idx]) s->display_max_m[display_idx] = positive_h;

                    /*
                       Ocean max-height / propagation review layer: use the actual
                       grid-cell water height.  Station/JESIM gauges still use the
                       separate 1 m shoreline-equivalent conversion, but the map
                       layer should show offshore propagation without converting
                       every ocean cell to a shoreline gauge basis.
                    */
                    if (s->display_arrival_s[display_idx] == 0 && positive_h > 0.0f) {
#if OCEAN_MAX_HEIGHT_ARRIVAL_SHORELINE_EQUIV_ENABLED
                        float arrival_basis_h = positive_h * station_shoreline_equiv_factor(visible_effective_depth_m(g, visible_idx));
#else
                        float arrival_basis_h = positive_h;
#endif

                        if (arrival_basis_h >= OCEAN_CELL_ARRIVAL_THRESHOLD_M) {
                            s->display_arrival_s[display_idx] = time_s;
                        }
                    }
                }
            }
        }
    }
}

void simulate_one_step(const GridData *g, Simulation *s, int time_s) {
    /* v125: true 15-second model tick.  Keep the stable v85 depth coefficient
       and use one internal numerical wave update per 15-second tick.  The loop
       remains so future substep experiments can be re-enabled safely. */
    int substeps = WAVE_INTERNAL_SUBSTEPS;
    if (substeps < 1) substeps = 1;
    for (int sub = 0; sub < substeps; ++sub) {
        simulate_one_numerical_substep(g, s, time_s);
    }
}

void update_station_stats(const GridData *g, Simulation *s, int time_s) {
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        int32_t idx = g->station_sample_index ? g->station_sample_index[i] : g->station_index[i];
        if (idx < 0 || (size_t)idx >= s->cells) continue;
        size_t p = (size_t)idx + s->width;
        /* Keep station_max_m as positive tsunami height for warning ranks,
           and define station arrival for JESIM/news output as the first
           meaningful POSITIVE water-level rise.  Do not let the initial
           source deformation drawdown/subsidence side lock the observation
           arrival time early. */
        float signed_h = station_observed_height_m(g, s, i);
        float positive_h = signed_h > 0.0f ? signed_h : 0.0f;
        if (positive_h > s->station_max_m[i]) s->station_max_m[i] = positive_h;
        if (s->station_arrival_s[i] == 0 && positive_h >= 0.05f) s->station_arrival_s[i] = time_s;
    }
}

void write_station_file(const char *path, const GridData *g, const Simulation *s) {
    FILE *fp = fopen(path, "wb");
    if (!fp) return;
    /* UTF-8 BOM helps Excel open Japanese forecast-zone names correctly. */
    fputs("\xEF\xBB\xBF", fp);
    fprintf(fp, "station_no,grid_index,grid_x,grid_y,sample_grid_index,sample_grid_x,sample_grid_y,forecast_zone_id,forecast_zone_name,max_height_m,arrival_s,color_bin,color_hex\n");
    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        int region_id = g->station_id[i];
        /* station_id=67 is non-forecast/outside. 1..66 are Japan JMA zones; 68..93 are KMA Korea coastal zones. */
        if (region_id < 1 || region_id > 93 || region_id == 67) continue;

        int32_t idx = g->station_index[i];
        if (idx < 0 || (size_t)idx >= s->cells) continue;
        int grid_x = idx % (int)g->width;
        int grid_y = idx / (int)g->width;
        int32_t sample_idx = g->station_sample_index ? g->station_sample_index[i] : idx;
        int sample_x = (sample_idx >= 0) ? (sample_idx % (int)g->width) : -1;
        int sample_y = (sample_idx >= 0) ? (sample_idx / (int)g->width) : -1;
        float max_h = s->station_max_m[i];
        fprintf(fp, "%u,%d,%d,%d,%d,%d,%d,%d,",
                i + 1, idx, grid_x, grid_y, sample_idx, sample_x, sample_y, region_id);
        csv_write_escaped(fp, forecast_zone_name_from_id(region_id));
        fprintf(fp, ",%.6g,%d,%s,%s\n",
                max_h, s->station_arrival_s[i],
                tsunami_height_color_bin(max_h), tsunami_height_color_hex(max_h));
    }
    fclose(fp);
}

void write_forecast_zone_file(const char *path, const GridData *g, const Simulation *s) {
    enum { MAX_ZONE_ID = 93 };
    float zone_max[MAX_ZONE_ID + 1];
    int zone_station_no[MAX_ZONE_ID + 1];
    int zone_grid_index[MAX_ZONE_ID + 1];
    int zone_arrival_s[MAX_ZONE_ID + 1];
    int zone_station_count[MAX_ZONE_ID + 1];

    for (int z = 0; z <= MAX_ZONE_ID; ++z) {
        zone_max[z] = 0.0f;
        zone_station_no[z] = 0;
        zone_grid_index[z] = -1;
        zone_arrival_s[z] = 0;
        zone_station_count[z] = 0;
    }

    for (uint32_t i = 0; i < g->station_count; ++i) {
        if (!station_is_enabled(g, i)) continue;
        int region_id = g->station_id[i];
        if (region_id < 1 || region_id > MAX_ZONE_ID || region_id == 67) continue;
        zone_station_count[region_id]++;
        float h = s->station_max_m[i];
        if (h > zone_max[region_id]) {
            zone_max[region_id] = h;
            zone_station_no[region_id] = (int)i + 1;
            zone_grid_index[region_id] = g->station_index[i];
            zone_arrival_s[region_id] = s->station_arrival_s[i];
        }
    }

    FILE *fp = fopen(path, "wb");
    if (!fp) return;
    fputs("\xEF\xBB\xBF", fp);
    fprintf(fp, "forecast_zone_id,forecast_zone_name,max_height_m,max_height_cm,station_no_of_max,grid_index,grid_x,grid_y,arrival_s,station_count,color_bin,color_hex\n");

    for (int region_id = 1; region_id <= MAX_ZONE_ID; ++region_id) {
        if (region_id == 67) continue;
        int idx = zone_grid_index[region_id];
        int grid_x = (idx >= 0) ? (idx % (int)g->width) : -1;
        int grid_y = (idx >= 0) ? (idx / (int)g->width) : -1;
        float max_h = zone_max[region_id];
        fprintf(fp, "%d,", region_id);
        csv_write_escaped(fp, forecast_zone_name_from_id(region_id));
        fprintf(fp, ",%.6g,%.6g,%d,%d,%d,%d,%d,%d,%s,%s\n",
                max_h, max_h * 100.0f, zone_station_no[region_id], idx, grid_x, grid_y,
                zone_arrival_s[region_id], zone_station_count[region_id],
                tsunami_height_color_bin(max_h), tsunami_height_color_hex(max_h));
    }

    /* station_id=67 is outside/non-forecast and intentionally not exported. */

    fclose(fp);
}

void source_list_append(SourceList *list, int grid_x, int grid_y, float height_m, int delay_s) {
    if (list->count == list->capacity) {
        size_t new_capacity = list->capacity ? list->capacity * 2u : 64u;
        SourceSpec *new_items = (SourceSpec *)realloc(list->items, new_capacity * sizeof(SourceSpec));
        if (!new_items) tsunami_fatal("out of memory");
        list->items = new_items;
        list->capacity = new_capacity;
    }
    list->items[list->count].grid_x = grid_x;
    list->items[list->count].grid_y = grid_y;
    list->items[list->count].height_m = height_m;
    list->items[list->count].slip_m = (height_m > 0.0005f) ? height_m : 0.0f; /* fallback; GUI overwrites uplift/footprint slip when available */
    list->items[list->count].stamp_dip_deg = 0.0f;
    list->items[list->count].moment_area_factor = 1.0f;
    list->items[list->count].fault_id = 0;
    list->items[list->count].fault_patch_id = 0;
    list->items[list->count].stamp_strike_deg = 0.0f;
    list->items[list->count].stamp_depth_km = 0.0f;
    list->items[list->count].stamp_width_km = 0.0f;
    list->items[list->count].stamp_length_km = 0.0f;
    list->items[list->count].delay_s = delay_s;
    list->count++;
}

void source_list_clear(SourceList *list) {
    list->count = 0;
}

void source_list_free(SourceList *list) {
    free(list->items);
    memset(list, 0, sizeof(*list));
}

int parse_source_csv_line(const char *text, int *x, int *y, float *height_m, int *delay_s) {
    char extra = '\0';
    return sscanf(text, " %d , %d , %f , %d %c", x, y, height_m, delay_s, &extra) == 4;
}

int load_sources_csv(SourceList *list, const char *path) {
    FILE *fp = fopen(path, "r");
    if (!fp) return 0;
    char line[512];
    int line_no = 0;
    while (fgets(line, sizeof(line), fp)) {
        line_no++;
        char *p = line;
        while (*p == ' ' || *p == '\t') p++;
        if (*p == '\0' || *p == '\n' || *p == '#') continue;
        if (line_no == 1 && strstr(p, "grid_x") != NULL) continue;
        int x = 0, y = 0, delay_s = 0;
        float h = 0.0f;
        if (parse_source_csv_line(p, &x, &y, &h, &delay_s)) {
            source_list_append(list, x, y, h, delay_s);
        }
    }
    fclose(fp);
    return 1;
}

void apply_sources_to_simulation(const GridData *g, Simulation *s, const SourceList *sources) {
    for (size_t i = 0; i < sources->count; ++i) {
        const SourceSpec *src = &sources->items[i];
        add_source_grid(g, s, src->grid_x, src->grid_y, src->height_m, src->delay_s);
    }
}
