import { INF, MAP_H, MAP_W, SIZE, clamp, createMapFields, fbm, hash2, indexOf, inside, lerp, pickEntities, rand, smoothstep, valueNoise } from "./mapUtils.js";
import {
  aStar,
  extractMaskBorder,
  extractRegionBorderSegments,
  generateRegionalPrefectures,
  makePrefectureMask,
  neighbors8,
} from "./mapGeneratorHelpers.js";

export function buildTerrainTemplate(seed) {
  const deposition = 0.18 + rand(seed, 41) * 0.72;
  const erosion = 0.24 + rand(seed, 42) * 0.68;
  const roughness = 0.34 + rand(seed, 43) * 0.62;
  const coastAxisPick = Math.floor(rand(seed, 10) * 3);
  const coastAngle = coastAxisPick === 0
    ? Math.PI / 2
    : coastAxisPick === 1
      ? 0
      : (rand(seed, 11) > 0.5 ? Math.PI / 4 : -Math.PI / 4) + (rand(seed, 14) - 0.5) * 0.28;
  const ridgeJaggedness = 0.20 + rand(seed, 44) * 0.70;
  const spineCount = 2 + Math.floor(rand(seed, 45) * 2);
  const sideAPlain = 0.035 + rand(seed, 56) * 0.115 + deposition * 0.085;
  const sideBPlain = 0.035 + rand(seed, 57) * 0.115 + deposition * 0.085;
  const spineAngle = coastAngle + Math.PI * (0.28 + rand(seed, 46) * 0.44);
  const spineCurve = (rand(seed, 47) - 0.5) * 0.28;
  const spinePosition = (rand(seed, 48) - 0.5) * 0.56;
  const backboneLongShift = (rand(seed, 68) - 0.5) * 0.34;
  const backboneCenterX = clamp(0.5 + Math.cos(spineAngle + Math.PI / 2) * spinePosition * 0.36 + Math.cos(spineAngle) * backboneLongShift, 0.18, 0.82);
  const backboneCenterY = clamp(0.5 + Math.sin(spineAngle + Math.PI / 2) * spinePosition * 0.36 + Math.sin(spineAngle) * backboneLongShift, 0.18, 0.82);
  const backboneLength = 0.50 + rand(seed, 69) * 0.28;
  const backboneWidth = 0.12 + rand(seed, 83) * 0.12;
  const backboneScratchCount = 34 + Math.floor(rand(seed, 84) * 30);

  return {
    seed,
    spineCount,
    spineAngle,
    spineCurve,
    spinePosition,
    backboneCenterX,
    backboneCenterY,
    backboneLength,
    backboneWidth,
    backboneScratchCount,
    spineStrength: 0.56 + rand(seed, 49) * 0.32,
    spineWidth: 0.034 + rand(seed, 50) * 0.036,
    // v4: 個別の丸い山塊生成を主役にしない。山地は下の folded orogeny field で一括生成する。
    secondaryMountainCount: 0,
    secondaryMountainSize: 0.038 + rand(seed, 52) * 0.060,
    secondaryMountainStrength: 0.40 + rand(seed, 53) * 0.25,
    rangeBreakCount: 4 + Math.floor(rand(seed, 62) * 4),
    rangeBreakWidth: 0.022 + rand(seed, 63) * 0.026,
    rangeBreakStrength: 0.060 + rand(seed, 64) * 0.070,
    plainNoiseSuppression: 0.34 + rand(seed, 65) * 0.22,
    // 高標高がすぐ天井へ張り付いて平頂山化しないよう、
    // ソフトクリップ開始をやや遅らせ、肩を高めに設定する。
    peakSoftStart: 0.905 + rand(seed, 66) * 0.030,
    peakSoftCap: 1.010 + rand(seed, 67) * 0.020,
    orographicStrength: 0.88 + rand(seed, 70) * 0.28,
    orographicCoverage: 0.72 + rand(seed, 71) * 0.18,
    foldDensity: 5.2 + rand(seed, 72) * 2.2,
    foldSharpness: 1.65 + rand(seed, 73) * 0.85,
    fluvialAggression: 1.50 + rand(seed, 74) * 0.48,
    // 粗い格子でも山肌の複雑さが出るよう、
    // 細谷の「見え」は手続きノイズ寄りにし、明示的な小流路本数は抑える。
    drainageDensity: 0.56 + rand(seed, 75) * 0.22,
    gullyIncision: 0.18 + rand(seed, 76) * 0.10,
    dendriticTexture: 0.72 + rand(seed, 77) * 0.24,
    macroNoiseStrength: 0.018 + rand(seed, 85) * 0.014,
    macroNoiseScale: 0.028 + rand(seed, 86) * 0.018,
    globalNoiseStrength: 0.010 + rand(seed, 87) * 0.010,
    headwaterGullyCount: 0,
    alpineMicroRelief: 0.018 + rand(seed, 79) * 0.022,
    basinOutletCount: 10 + Math.floor(rand(seed, 80) * 8),
    meanderStrength: 0.18 + rand(seed, 81) * 0.16,
    terrainNoisePatchiness: 0.42 + rand(seed, 82) * 0.28,
    coastAxis: coastAxisPick === 0 ? "east-west" : coastAxisPick === 1 ? "north-south" : "diagonal",
    coastAngle,
    coastBias: 0.18 + rand(seed, 12) * 0.24,
    coastRoughness: 0.34 + rand(seed, 54) * 0.58,
    coastSides: [
      {
        penetration: 0.18 + rand(seed, 58) * 0.16,
        inletStrength: 0.18 + rand(seed, 59) * 0.56,
        plainWidth: sideAPlain,
      },
      {
        penetration: 0.18 + rand(seed, 60) * 0.16,
        inletStrength: 0.18 + rand(seed, 61) * 0.56,
        plainWidth: sideBPlain,
      },
    ],
    deposition,
    erosion,
    roughness,
    ridgeJaggedness,
    ridgeBranchiness: 0.44 + rand(seed, 55) * 0.66,
    detachedRangeCount: 0,
    alpinePeakCount: 0,
  };
}

function jaggedRidgeContribution(x, y, ridge, seed) {
  const dx = x - ridge.x;
  const dy = y - ridge.y;
  const ca = Math.cos(ridge.angle);
  const sa = Math.sin(ridge.angle);
  const along = dx * ca + dy * sa;
  const perp = -dx * sa + dy * ca;
  const nAlong = along / Math.max(0.001, ridge.length);
  const lengthFade = smoothstep(1 - Math.abs(nAlong));
  if (lengthFade <= 0) return 0;

  // Bend the centerline itself with coherent long/mid waves, then apply ridge falloff.
  const low = (valueNoise(along * 0.85 + ridge.seedOffset, ridge.seedOffset * 0.37, seed + 6100, 28) - 0.5) * 2;
  const mid = (valueNoise(along * 1.7 - ridge.seedOffset, ridge.seedOffset * 0.23, seed + 6200, 13) - 0.5) * 2;
  const sine = Math.sin(along * ridge.kinkFrequency + ridge.kinkPhase);
  const curve = (ridge.curve || 0) * along * along * (along >= 0 ? 1 : -1);
  const axisOffset = low * ridge.axisWobble + mid * ridge.axisWobble * 0.55 + sine * ridge.axisWobble * 0.25 + curve;
  const widthNoise = 0.78 + valueNoise(along * 1.2 + ridge.seedOffset, ridge.seedOffset * 0.19, seed + 6300, 21) * ridge.widthVariation;
  const localWidth = Math.max(0.006, ridge.width * widthNoise);
  const jaggedPerp = perp - axisOffset;
  const serration = 0.76 + valueNoise(x * 1.1 + along * 0.18, y * 1.1 + perp * 0.18, seed + ridge.seedOffset, 7) * 0.48;
  return Math.exp(-(jaggedPerp * jaggedPerp) / (localWidth * localWidth)) * lengthFade * ridge.h * serration;
}

function spineFieldAt(x, y, template, spineIndex) {
  const seed = template.seed || 0;
  const spacing = spineIndex === 0 ? 0 : (spineIndex % 2 ? 0.30 : -0.30);
  const angle = template.spineAngle + (spineIndex - 0.5) * 0.17 + (rand(seed, 700 + spineIndex) - 0.5) * 0.18;
  const ridge = {
    x: 0.5 + Math.cos(angle + Math.PI / 2) * (template.spinePosition + spacing) * 0.45,
    y: 0.5 + Math.sin(angle + Math.PI / 2) * (template.spinePosition + spacing) * 0.45,
    angle,
    width: template.spineWidth * (0.82 + rand(seed, 710 + spineIndex) * 0.38),
    length: 0.78 + rand(seed, 720 + spineIndex) * 0.28,
    h: template.spineStrength * (0.24 + rand(seed, 730 + spineIndex) * 0.12),
    curve: template.spineCurve,
    axisWobble: template.spineWidth * (0.45 + template.ridgeJaggedness * 1.15),
    kinkFrequency: 10 + rand(seed, 740 + spineIndex) * 18,
    kinkPhase: rand(seed, 750 + spineIndex) * Math.PI * 2,
    seedOffset: 7600 + spineIndex * 211,
    widthVariation: 0.18 + template.ridgeJaggedness * 0.34,
  };
  return jaggedRidgeContribution(x, y, ridge, seed);
}

function buildSpineRidges(seed, template) {
  const spines = [];
  const branches = [];
  const centerX = template.backboneCenterX ?? 0.5;
  const centerY = template.backboneCenterY ?? 0.5;
  const baseAngle = template.spineAngle;
  const axisCos = Math.cos(baseAngle);
  const axisSin = Math.sin(baseAngle);
  const crossCos = Math.cos(baseAngle + Math.PI / 2);
  const crossSin = Math.sin(baseAngle + Math.PI / 2);
  const corridorLength = template.backboneLength ?? 0.62;
  const corridorWidth = template.backboneWidth ?? 0.16;

  for (let i = 0; i < template.spineCount; i++) {
    const along = ((i / Math.max(1, template.spineCount - 1)) - 0.5) * corridorLength * 0.48 + (rand(seed, 705 + i) - 0.5) * corridorLength * 0.12;
    const cross = (rand(seed, 706 + i) - 0.5) * corridorWidth * 0.45;
    const angle = baseAngle + (rand(seed, 700 + i) - 0.5) * 0.16;
    const x = clamp(centerX + axisCos * along + crossCos * cross, 0.06, 0.94);
    const y = clamp(centerY + axisSin * along + crossSin * cross, 0.06, 0.94);
    spines.push({
      x, y, angle,
      width: template.spineWidth * (0.52 + rand(seed, 710 + i) * 0.28),
      length: corridorLength * (0.58 + rand(seed, 720 + i) * 0.18),
      h: template.spineStrength * (0.11 + rand(seed, 730 + i) * 0.08),
      curve: template.spineCurve,
      axisWobble: template.spineWidth * (0.42 + template.ridgeJaggedness * 0.95),
      kinkFrequency: 10 + rand(seed, 740 + i) * 16,
      kinkPhase: rand(seed, 750 + i) * Math.PI * 2,
      seedOffset: 7600 + i * 211,
      widthVariation: 0.18 + template.ridgeJaggedness * 0.34,
    });
  }

  const scratchCount = template.backboneScratchCount ?? (36 + Math.floor(template.ridgeBranchiness * 30));
  for (let b = 0; b < scratchCount; b++) {
    const centerBias = ((rand(seed, 810 + b) + rand(seed, 811 + b)) * 0.5 - 0.5) * 2;
    const edgeBias = (rand(seed, 812 + b) - 0.5) * 2;
    const along = centerBias * corridorLength * 0.82;
    const corridorT = clamp(1 - Math.abs(centerBias));
    const lateralSpread = corridorWidth * (0.48 + corridorT * 0.72);
    const cross = edgeBias * lateralSpread;
    const radialT = clamp(1 - Math.abs(edgeBias));
    const density = clamp(corridorT * 0.72 + radialT * 0.28);
    const x = clamp(centerX + axisCos * along + crossCos * cross, 0.04, 0.96);
    const y = clamp(centerY + axisSin * along + crossSin * cross, 0.04, 0.96);
    const branchAngle = baseAngle + (rand(seed, 830 + b) - 0.5) * 0.92 + edgeBias * 0.20;
    branches.push({
      x,
      y,
      angle: branchAngle,
      width: template.spineWidth * (0.16 + rand(seed, 840 + b) * 0.16 + density * 0.08),
      length: 0.07 + rand(seed, 850 + b) * 0.13 + density * 0.06,
      h: template.spineStrength * (0.018 + density * 0.066 + rand(seed, 860 + b) * 0.018),
      curve: template.spineCurve * 0.22,
      axisWobble: template.spineWidth * (0.22 + template.ridgeJaggedness * 0.58),
      kinkFrequency: 15 + rand(seed, 870 + b) * 22,
      kinkPhase: rand(seed, 880 + b) * Math.PI * 2,
      seedOffset: 8800 + b * 37,
      widthVariation: 0.16 + template.ridgeJaggedness * 0.20,
    });
  }
  return { spines, branches };
}


function softUpperClamp(value, start = 0.8, cap = 0.96) {
  if (value <= start) return value;
  if (value <= cap) {
    const t = clamp((value - start) / Math.max(0.001, cap - start));
    // 肩へ向かって緩やかに圧縮するが、cap 未満ではなるべく差を残す。
    return lerp(value, start + (cap - start) * (1 - Math.pow(1 - t, 1.18)), 0.16);
  }
  const overflow = value - cap;
  // cap 超過分も少し残して、山頂が一様な平頂面にならないようにする。
  return cap + overflow * (0.28 / (1 + overflow * 4.2));
}

function elongatedFeatureContribution(x, y, feature, seed) {
  const dx = x - feature.x;
  const dy = y - feature.y;
  const ca = Math.cos(feature.angle);
  const sa = Math.sin(feature.angle);
  const along = dx * ca + dy * sa;
  const perp = -dx * sa + dy * ca;
  const nAlong = along / Math.max(0.001, feature.length);
  if (Math.abs(nAlong) > 1.35) return 0;
  const alongFade = Math.exp(-nAlong * nAlong * 1.7);
  const low = (valueNoise(along * 0.95 + feature.seedOffset, feature.seedOffset * 0.31, seed + 6400, 19) - 0.5) * 2;
  const mid = (valueNoise(along * 1.75 - feature.seedOffset, feature.seedOffset * 0.21, seed + 6500, 9) - 0.5) * 2;
  const axisOffset = low * feature.axisWobble + mid * feature.axisWobble * 0.45;
  const localWidth = Math.max(0.008, feature.width * (0.84 + valueNoise(along * 1.15, feature.seedOffset, seed + 6600, 14) * feature.widthVariation));
  const offsetPerp = perp - axisOffset;
  return Math.exp(-(offsetPerp * offsetPerp) / (localWidth * localWidth)) * alongFade * feature.h;
}

function buildRangeBreaks(seed, template, spines) {
  const rangeBreaks = [];
  for (let i = 0; i < spines.length; i++) {
    const spine = spines[i];
    const count = Math.max(2, template.rangeBreakCount - 1 + Math.floor(rand(seed, 890 + i) * 3));
    for (let b = 0; b < count; b++) {
      const along = (rand(seed, 900 + i * 37 + b) - 0.5) * spine.length * 0.84;
      const lateral = (rand(seed, 910 + i * 37 + b) - 0.5) * spine.width * 0.9;
      rangeBreaks.push({
        x: spine.x + Math.cos(spine.angle) * along + Math.cos(spine.angle + Math.PI / 2) * lateral,
        y: spine.y + Math.sin(spine.angle) * along + Math.sin(spine.angle + Math.PI / 2) * lateral,
        angle: spine.angle + (rand(seed, 920 + i * 37 + b) > 0.5 ? Math.PI / 2 : -Math.PI / 2) + (rand(seed, 930 + i * 37 + b) - 0.5) * 0.42,
        width: template.rangeBreakWidth * (0.75 + rand(seed, 940 + i * 37 + b) * 0.75),
        length: 0.12 + rand(seed, 950 + i * 37 + b) * 0.14,
        h: template.rangeBreakStrength * (0.78 + rand(seed, 960 + i * 37 + b) * 0.55),
        axisWobble: template.rangeBreakWidth * (0.18 + rand(seed, 970 + i * 37 + b) * 0.32),
        widthVariation: 0.14 + rand(seed, 980 + i * 37 + b) * 0.24,
        seedOffset: 9900 + i * 311 + b * 41,
      });
    }
  }
  return rangeBreaks;
}


function buildDetachedRanges(seed, template) {
  const ranges = [];
  const count = template.detachedRangeCount ?? 6;
  for (let i = 0; i < count; i++) {
    const quadrantX = i % 2 === 0 ? 0.24 : 0.76;
    const quadrantY = Math.floor(i / 2) % 2 === 0 ? 0.24 : 0.76;
    const free = rand(seed, 12000 + i) < 0.45;
    const x = free ? 0.12 + rand(seed, 12010 + i) * 0.76 : quadrantX + (rand(seed, 12020 + i) - 0.5) * 0.28;
    const y = free ? 0.12 + rand(seed, 12030 + i) * 0.76 : quadrantY + (rand(seed, 12040 + i) - 0.5) * 0.28;
    const angle = template.spineAngle + (rand(seed, 12050 + i) - 0.5) * Math.PI * 0.95;
    ranges.push({
      x: clamp(x, 0.08, 0.92),
      y: clamp(y, 0.08, 0.92),
      angle,
      width: 0.020 + rand(seed, 12060 + i) * 0.030,
      length: 0.18 + rand(seed, 12070 + i) * 0.28,
      h: 0.075 + rand(seed, 12080 + i) * 0.095,
      curve: (rand(seed, 12090 + i) - 0.5) * 0.10,
      axisWobble: 0.018 + template.ridgeJaggedness * 0.030,
      kinkFrequency: 14 + rand(seed, 12100 + i) * 24,
      kinkPhase: rand(seed, 12110 + i) * Math.PI * 2,
      seedOffset: 12120 + i * 173,
      widthVariation: 0.28 + template.ridgeJaggedness * 0.36,
    });
  }
  return ranges;
}

function buildAlpinePeaks(seed, template, detachedRanges) {
  const peaks = [];
  const count = template.alpinePeakCount ?? 8;
  for (let i = 0; i < count; i++) {
    const attach = detachedRanges.length && rand(seed, 12300 + i) < 0.62;
    const base = attach ? detachedRanges[i % detachedRanges.length] : null;
    const along = base ? (rand(seed, 12310 + i) - 0.5) * base.length * 0.90 : 0;
    const perp = base ? (rand(seed, 12320 + i) - 0.5) * base.width * 4.5 : 0;
    const x = base ? base.x + Math.cos(base.angle) * along + Math.cos(base.angle + Math.PI / 2) * perp : 0.10 + rand(seed, 12330 + i) * 0.80;
    const y = base ? base.y + Math.sin(base.angle) * along + Math.sin(base.angle + Math.PI / 2) * perp : 0.10 + rand(seed, 12340 + i) * 0.80;
    peaks.push({
      x: clamp(x, 0.06, 0.94),
      y: clamp(y, 0.06, 0.94),
      angle: base ? base.angle + (rand(seed, 12350 + i) - 0.5) * 0.9 : rand(seed, 12360 + i) * Math.PI * 2,
      rx: 0.022 + rand(seed, 12370 + i) * 0.035,
      ry: 0.012 + rand(seed, 12380 + i) * 0.024,
      h: 0.070 + rand(seed, 12390 + i) * 0.100,
      seedOffset: 12400 + i * 191,
    });
  }
  return peaks;
}

// v5: 「全域を海底として初期化し、海底から広域隆起で山地を生やす」ための造山場。
// 周期的な褶曲波はワッフル状の縞を作るため廃止し、広い帯状隆起・不均質な断層谷・
// 連続した尾根核を domain-warp 付きで合成する。
function foldedOrogenyAt(px, py, seed, template, coastLower = 0) {
  const baseAngle = template.spineAngle + (rand(seed, 13001) - 0.5) * 0.28;
  const warpX = (fbm(px * 2.0 + 17, py * 2.0 - 31, seed + 13010) - 0.5) * 0.20;
  const warpY = (fbm(px * 2.1 - 43, py * 2.1 + 19, seed + 13020) - 0.5) * 0.20;
  const x = px + warpX;
  const y = py + warpY;

  let beltMass = 0;
  let ridgeCores = 0;
  let structuralValleys = 0;
  const beltCount = 3;

  for (let k = 0; k < beltCount; k++) {
    const angle = baseAngle + (k - 1) * 0.24 + (rand(seed, 13100 + k) - 0.5) * 0.26;
    const ca = Math.cos(angle);
    const sa = Math.sin(angle);
    const along = x * ca + y * sa;
    const cross = -x * sa + y * ca;
    const offset = (rand(seed, 13120 + k) - 0.5) * 0.38 + (k - 1) * 0.10;
    const axisWarp = (fbm(px * 1.65 + k * 11, py * 1.65 - k * 7, seed + 13200 + k) - 0.5) * (0.15 + template.ridgeJaggedness * 0.08);
    const lengthWarp = (valueNoise(px * 1.3 - k * 17, py * 1.3 + k * 13, seed + 13250 + k, 2.4) - 0.5) * 0.34;
    const localCross = cross - offset - axisWarp;
    const localAlong = along + lengthWarp;
    const width = 0.24 + rand(seed, 13300 + k) * 0.14;
    const coreWidth = 0.058 + rand(seed, 13320 + k) * 0.046;
    const lengthFade = smoothstep(1.38 - Math.abs(localAlong - 0.52));
    const broad = Math.exp(-(localCross * localCross) / (width * width)) * lengthFade;
    const core = Math.exp(-(localCross * localCross) / (coreWidth * coreWidth)) * lengthFade;
    const broken = 0.72 + valueNoise(px * 2.2 + k * 21, py * 2.2 - k * 15, seed + 13400 + k, 2.1) * 0.46;
    const rugged = 0.74 + fbm(px * 4.3 + k * 19, py * 4.3 - k * 23, seed + 13480 + k) * 0.54;
    beltMass += broad * (0.285 + k * 0.040) * broken;
    ridgeCores += core * (0.340 + k * 0.030) * rugged;

    // 河川が後で選びやすい弱線。周期的な縞ではなく、ノイズで途切れる断層・構造谷として扱う。
    const valleyNoise = Math.max(0, valueNoise(px * 5.2 + k * 31, py * 5.2 - k * 27, seed + 13540 + k, 3.2) - 0.55);
    const transverse = Math.exp(-((localAlong - (0.26 + rand(seed, 13600 + k) * 0.52)) ** 2) / 0.030);
    structuralValleys += broad * (valleyNoise * 0.18 + transverse * 0.035);
  }

  const broadA = fbm(px * 0.80 + 23, py * 0.80 - 61, seed + 13700);
  const broadB = valueNoise(px * 1.35 - 41, py * 1.35 + 17, seed + 13710, 2.4);
  const continentalPulse = clamp((broadA * 0.58 + broadB * 0.42 - 0.22) / 0.66);
  const edgeDistance = Math.min(px, py, 1 - px, 1 - py);
  const edgeDrowning = 1 - smoothstep(edgeDistance / 0.16);
  const marineLoss = edgeDrowning * (0.22 + template.deposition * 0.10) + coastLower * 0.035;

  const mass = clamp(beltMass * 0.95 + continentalPulse * 0.40 + template.orographicCoverage * 0.12 - structuralValleys * 0.78 - marineLoss * 0.18);
  const ridges = clamp(ridgeCores * 1.18 + mass * 0.20 - structuralValleys * 0.64);
  const uplift = clamp((mass * 0.62 + ridges * 0.34) * template.orographicStrength);
  return { uplift, ridges, valleys: clamp(structuralValleys * 5.0) };
}


// v6: 谷状の細かな侵食テクスチャ。流路網のない場所も「沢に削られた山肌」に見せる。
// 高周波ノイズをそのまま標高に足すのではなく、ノイズ場の細い等値線だけを抽出して
// 傾斜・標高・湿潤度でマスクする。
function dendriticRavineTexture(x, y, seed) {
  const warpX = (fbm(x * 0.72 + 113, y * 0.72 - 71, seed + 15010) - 0.5) * 9.5;
  const warpY = (fbm(x * 0.74 - 59, y * 0.74 + 97, seed + 15020) - 0.5) * 9.5;
  const wx = x + warpX;
  const wy = y + warpY;

  const n1 = valueNoise(wx * 1.15 + 31, wy * 1.15 - 47, seed + 15100, 8.0);
  const n2 = valueNoise(wx * 2.05 - 19, wy * 2.05 + 23, seed + 15110, 4.8);
  const n3 = valueNoise(wx * 3.65 + 71, wy * 3.65 - 11, seed + 15120, 3.0);
  const line1 = Math.pow(clamp(1 - Math.abs(n1 - 0.50) * 6.1), 3.4);
  const line2 = Math.pow(clamp(1 - Math.abs(n2 - 0.47) * 7.6), 3.0);
  const line3 = Math.pow(clamp(1 - Math.abs(n3 - 0.52) * 9.0), 2.6);

  const join = valueNoise(wx * 0.55 + 7, wy * 0.55 - 5, seed + 15140, 11.0);
  return clamp(line1 * 0.48 + line1 * line2 * 0.44 + line2 * line3 * 0.26 + line3 * Math.max(0, join - 0.58) * 0.36);
}

function addPathIncision(field, path, strength, radius = 1) {
  if (!path || path.length < 2) return;
  for (let k = 0; k < path.length; k++) {
    const [px, py] = path[k];
    const downstream = k / Math.max(1, path.length - 1);
    const local = strength * (0.72 + downstream * 0.44);
    for (let dy = -radius; dy <= radius; dy++) {
      for (let dx = -radius; dx <= radius; dx++) {
        const nx = px + dx;
        const ny = py + dy;
        if (!inside(nx, ny)) continue;
        const d = Math.hypot(dx, dy);
        if (d > radius + 0.15) continue;
        const w = d < 0.001 ? 1 : Math.max(0, 1 - d / (radius + 0.35)) * 0.52;
        const i = indexOf(nx, ny);
        field[i] = clamp(field[i] + local * w);
      }
    }
  }
}

// v7: line-only carving is too thin and looks like a texture overlay.  This helper
// converts traced drainage paths into a distance field with a V-shaped cross-section:
// small headwater gullies stay narrow, while longer/merged channels open a slightly
// wider valley floor.  The visual result is produced by terrain change, not by merely
// drawing more blue river lines.
function addValleyDistanceInfluence(incisionField, floorField, path, strength, radius = 2, floorRadius = 0.55) {
  if (!path || path.length < 2) return;
  for (let k = 0; k < path.length; k++) {
    const [px, py] = path[k];
    const downstream = k / Math.max(1, path.length - 1);
    const localRadius = Math.max(1.1, radius * (0.68 + downstream * 0.56));
    const localStrength = strength * (0.62 + downstream * 0.70);
    const r = Math.ceil(localRadius + 1.2);
    for (let dy = -r; dy <= r; dy++) {
      for (let dx = -r; dx <= r; dx++) {
        const nx = px + dx;
        const ny = py + dy;
        if (!inside(nx, ny)) continue;
        const d = Math.hypot(dx, dy);
        if (d > localRadius + 1.0) continue;
        const i = indexOf(nx, ny);
        const vShape = Math.pow(clamp(1 - d / (localRadius + 0.55)), 1.55);
        const floor = smoothstep((floorRadius + 0.25 - d) / Math.max(0.35, floorRadius + 0.25));
        incisionField[i] = clamp(incisionField[i] + localStrength * vShape * 0.48);
        floorField[i] = clamp(floorField[i] + localStrength * floor * 0.42);
      }
    }
  }
}

function applyAlpineMicroRelief(elevation, sea, slope, ridgeField, valleyField, coastalLowland, seaLevel, seed, terrainTemplate, surfaceTextureField = null) {
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const highland = clamp((elevation[i] - (seaLevel + 0.12)) / 0.36);
      const alpineMask = clamp((elevation[i] - (seaLevel + 0.20)) / 0.28) * clamp(slope[i] * 1.15 + ridgeField[i] * 0.44 - valleyField[i] * 0.26) * (1 - coastalLowland[i] * 0.75);
      const ruggedMask = clamp(highland * (0.28 + slope[i] * 0.95 + ridgeField[i] * 0.38 - valleyField[i] * 0.16));
      const patch = clamp(0.35 + (fbm(x * 0.055 + 80, y * 0.055 - 34, seed + 18018) - 0.5) * 1.45 + terrainTemplate.terrainNoisePatchiness * 0.22);
      const activeMask = Math.max(alpineMask, ruggedMask * 0.56) * (0.52 + patch * 0.48);
      if (activeMask <= 0.02) continue;
      const warpX = x + (fbm(x * 0.18 + 24, y * 0.18 - 17, seed + 18021) - 0.5) * 5.5;
      const warpY = y + (fbm(x * 0.18 - 37, y * 0.18 + 13, seed + 18022) - 0.5) * 5.5;
      const coarse = (valueNoise(warpX * 0.72, warpY * 0.72, seed + 18023, 4.8) - 0.5) * 2;
      const medium = (valueNoise(warpX * 1.18 - 11, warpY * 1.18 + 19, seed + 18024, 6.7) - 0.5) * 2;
      const fine = (valueNoise(warpX * 1.95 + 17, warpY * 1.95 - 9, seed + 18025, 9.4) - 0.5) * 2;
      const ridged = 1 - Math.abs((valueNoise(warpX * 1.36 - 7, warpY * 1.36 + 21, seed + 18026, 6.2) - 0.5) * 2);
      const cellular = (fbm(warpX * 0.52 + 9, warpY * 0.52 - 4, seed + 18027) - 0.5) * 2;
      const perturb = (coarse * 0.42 + medium * 0.30 + fine * 0.18 + (ridged - 0.5) * 0.88 + cellular * 0.22) * terrainTemplate.alpineMicroRelief * activeMask;
      elevation[i] = clamp(elevation[i] + perturb, seaLevel + 0.006, 0.998);
      if (surfaceTextureField) {
        const textureSignal = clamp(Math.abs(perturb) * 44 + ridged * activeMask * 0.40 + Math.abs(cellular) * activeMask * 0.24 + patch * activeMask * 0.16);
        surfaceTextureField[i] = clamp(surfaceTextureField[i] + textureSignal);
      }
      if (perturb < 0) valleyField[i] = clamp(valleyField[i] + (-perturb) * activeMask * 4.5);
    }
  }
}


function breakHighPlateaus(elevation, sea, slope, ridgeField, valleyField, flowAccum, coastalLowland, seaLevel, seed, terrainTemplate, visibleRavineField = null, surfaceTextureField = null) {
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const high = clamp((elevation[i] - (seaLevel + 0.18)) / 0.28);
      const flat = clamp((0.26 - slope[i]) * 5.2);
      const mountain = clamp(ridgeField[i] * 0.78 + high * 0.62 - valleyField[i] * 0.16) * (1 - coastalLowland[i] * 0.82);
      const active = high * flat * mountain;
      if (active <= 0.02) continue;

      const warpX = x + (fbm(x * 0.11 + 13, y * 0.11 - 17, seed + 18101) - 0.5) * 4.6;
      const warpY = y + (fbm(x * 0.11 - 29, y * 0.11 + 7, seed + 18102) - 0.5) * 4.6;
      const dend = dendriticRavineTexture(warpX * 1.05, warpY * 1.05, seed + 18103);
      const broad = Math.abs((fbm(warpX * 0.16 + 4, warpY * 0.16 - 9, seed + 18104) - 0.5) * 2);
      const summit = Math.max(0, valueNoise(warpX * 1.85 + 11, warpY * 1.85 - 23, seed + 18105, 4.7) - 0.57);
      const drainage = Math.pow(flowAccum[i], 0.50);
      const carve = active * (0.0038 + dend * 0.0062 + broad * 0.0026 + drainage * 0.0032) * (0.92 + terrainTemplate.dendriticTexture * 0.30);
      const bump = active * summit * (0.0028 + terrainTemplate.roughness * 0.0022);
      elevation[i] = clamp(elevation[i] - carve + bump, seaLevel + 0.006, 0.998);
      valleyField[i] = clamp(valleyField[i] + carve * 18.0);
      ridgeField[i] = clamp(ridgeField[i] + bump * 8.0);
      if (visibleRavineField) visibleRavineField[i] = clamp(visibleRavineField[i] + carve * 8.6);
      if (surfaceTextureField) surfaceTextureField[i] = clamp(surfaceTextureField[i] + active * (0.16 + dend * 0.22 + broad * 0.10));
    }
  }
}

function carveOutletChannel(elevation, sea, lake, river, valleyField, basinField, flowAccum, seaLevel, startIndex, targetIndex, seed, bonusSeed = 0) {
  const start = { x: startIndex % MAP_W, y: Math.floor(startIndex / MAP_W) };
  const goal = { x: targetIndex % MAP_W, y: Math.floor(targetIndex / MAP_W) };
  const path = aStar(start, goal, (x, y, cx, cy) => {
    const i = indexOf(x, y);
    const ci = indexOf(cx, cy);
    if (sea[i]) return 0.12;
    const uphill = Math.max(0, elevation[i] - elevation[ci]);
    return Math.max(0.16, 0.55 + uphill * 80 + Math.max(0, elevation[i] - seaLevel) * 0.10 - flowAccum[i] * 0.65 - valleyField[i] * 0.38 - basinField[i] * 0.12 + (hash2(x, y, seed + bonusSeed) - 0.5) * 0.05);
  });
  if (path.length < 2) return [];
  const startElev = elevation[startIndex];
  const targetElev = sea[targetIndex] ? seaLevel - 0.002 : Math.min(startElev - 0.010, elevation[targetIndex]);
  for (let k = 0; k < path.length; k++) {
    const [x, y] = path[k];
    const i = indexOf(x, y);
    const t = k / Math.max(1, path.length - 1);
    const base = lerp(startElev - 0.001, targetElev + 0.004, t);
    const wiggle = (hash2(x, y, seed + 22000 + bonusSeed) - 0.5) * 0.0025;
    const floorLimit = seaLevel + 0.020 + basinField[i] * 0.014;
    elevation[i] = Math.min(elevation[i], Math.max(floorLimit, base + wiggle));
    valleyField[i] = clamp(valleyField[i] + 0.20 + (1 - t) * 0.12);
    basinField[i] = Math.max(0, basinField[i] - 0.10);
    river[i] = Math.max(river[i], 0.18 + flowAccum[i] * 0.42 + t * 0.10);
    if (lake[i] && k < path.length - 1) lake[i] = 0;
    if (!sea[i]) {
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (sea[ni]) continue;
        const d = Math.hypot(nx - x, ny - y);
        const widen = Math.max(0, 0.010 - d * 0.003);
        if (widen > 0) {
          elevation[ni] = Math.min(elevation[ni], Math.max(seaLevel + 0.020, elevation[i] + 0.008 + d * 0.004));
          valleyField[ni] = clamp(valleyField[ni] + widen * 12);
        }
      }
    }
  }
  return path;
}

export function generateTerrainAndRivers(seed) {
  let prefectureMask;
  let prefectureBorder;

  const {
    elevation,
    moisture,
    slope,
    sea,
    ocean,
    lake,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    basinField,
    coastalLowland,
    flowAccum,
    erosionField,
    depositionField,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    naturalBarrierScore,
    flowTo,
    portSuitability,
    crossingSuitability,
    passSuitability,
  } = createMapFields();

  const terrainTemplate = buildTerrainTemplate(seed);
  const coastAngle = terrainTemplate.coastAngle;
  const coastX = Math.cos(coastAngle);
  const coastY = Math.sin(coastAngle);
  const coastThreshold = terrainTemplate.coastBias;
  const coastStrength = 0.10 + (1 - terrainTemplate.deposition) * 0.12 + rand(seed, 13) * 0.09;
  const { spines, branches } = buildSpineRidges(seed, terrainTemplate);
  const detachedRanges = buildDetachedRanges(seed, terrainTemplate);
  const alpinePeaks = buildAlpinePeaks(seed, terrainTemplate, detachedRanges);
  const rangeBreaks = buildRangeBreaks(seed, terrainTemplate, spines);

  function coastPressureAt(x, y, wx = x, wy = y) {
    const nx = x / (MAP_W - 1) - 0.5;
    const ny = y / (MAP_H - 1) - 0.5;
    const axis = nx * coastX + ny * coastY;
    const waveA = (fbm(wx * 0.72 + 31, wy * 0.72 - 17, seed + 2222) - 0.5) * (0.05 + terrainTemplate.coastRoughness * terrainTemplate.coastSides[0].inletStrength * 0.18) +
      (valueNoise(wx + 19, wy - 23, seed + 2233, 18) - 0.5) * (0.03 + terrainTemplate.coastSides[0].inletStrength * 0.10);
    const waveB = (fbm(wx * 0.68 - 41, wy * 0.68 + 29, seed + 3222) - 0.5) * (0.05 + terrainTemplate.coastRoughness * terrainTemplate.coastSides[1].inletStrength * 0.18) +
      (valueNoise(wx - 13, wy + 37, seed + 3233, 16) - 0.5) * (0.03 + terrainTemplate.coastSides[1].inletStrength * 0.10);
    const sideA = smoothstep((axis + waveA - (0.50 - terrainTemplate.coastSides[0].penetration)) / Math.max(0.08, terrainTemplate.coastSides[0].plainWidth * 2.4));
    const sideB = smoothstep((-axis + waveB - (0.50 - terrainTemplate.coastSides[1].penetration)) / Math.max(0.08, terrainTemplate.coastSides[1].plainWidth * 2.4));
    return { sideA, sideB, pressure: Math.max(sideA, sideB), signedAxis: axis };
  }

  const seaLevel = 0.275;

  const mountainBlobs = Array.from({ length: terrainTemplate.secondaryMountainCount }, (_, i) => {
    const spine = spines[i % spines.length];
    const nearSpine = rand(seed, 98 + i) < 0.72;
    const edgeBias = rand(seed, 99 + i) < 0.28;
    const along = (rand(seed, 100 + i) - 0.5) * spine.length * 0.95;
    const side = rand(seed, 101 + i) > 0.5 ? 1 : -1;
    const offset = (0.055 + rand(seed, 102 + i) * 0.22) * side;
    let x = nearSpine ? spine.x + Math.cos(spine.angle) * along + Math.cos(spine.angle + Math.PI / 2) * offset : rand(seed, 103 + i);
    let y = nearSpine ? spine.y + Math.sin(spine.angle) * along + Math.sin(spine.angle + Math.PI / 2) * offset : rand(seed, 104 + i);
    if (edgeBias) {
      const edgeSide = Math.floor(rand(seed, 105 + i) * 4);
      if (edgeSide === 0) x = Math.min(x, 0.08 + rand(seed, 106 + i) * 0.10);
      if (edgeSide === 1) x = Math.max(x, 0.92 - rand(seed, 107 + i) * 0.10);
      if (edgeSide === 2) y = Math.min(y, 0.08 + rand(seed, 108 + i) * 0.10);
      if (edgeSide === 3) y = Math.max(y, 0.92 - rand(seed, 109 + i) * 0.10);
    }
    const coastSide = (x - 0.5) * coastX + (y - 0.5) * coastY;
    const mountainSide = coastSide >= 0 ? 1 : -1;
    if (rand(seed, 110 + i) < 0.46 && Math.abs(coastSide) > 0.28 - coastThreshold * 0.35) {
      x -= coastX * mountainSide * (0.05 + rand(seed, 111 + i) * 0.11);
      y -= coastY * mountainSide * (0.05 + rand(seed, 112 + i) * 0.11);
    }
    const angle = nearSpine ? spine.angle + (rand(seed, 302 + i) - 0.5) * 0.75 : rand(seed, 303 + i) * Math.PI * 2;
    const baseRadius = terrainTemplate.secondaryMountainSize * Math.min(MAP_W, MAP_H);
    return {
      x: clamp(x) * MAP_W,
      y: clamp(y) * MAP_H,
      angle,
      rx: baseRadius * (0.95 + rand(seed, 300 + i) * 1.10),
      ry: baseRadius * (0.34 + rand(seed, 301 + i) * 0.46),
      h: terrainTemplate.secondaryMountainStrength * (0.11 + rand(seed, 400 + i) * 0.23),
    };
  });

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const nx = x / (MAP_W - 1) - 0.5;
      const ny = y / (MAP_H - 1) - 0.5;
      const i = indexOf(x, y);

      const warpX = (fbm(x * 0.62 + 180, y * 0.62 - 90, seed + 3101) - 0.5) * 13;
      const warpY = (fbm(x * 0.62 - 70, y * 0.62 + 210, seed + 3201) - 0.5) * 13;
      const wx = x + warpX;
      const wy = y + warpY;

      let mountains = 0;
      for (const blob of mountainBlobs) {
        const dx = wx - blob.x;
        const dy = wy - blob.y;
        const ca = Math.cos(blob.angle);
        const sa = Math.sin(blob.angle);
        const along = (dx * ca + dy * sa) / Math.max(1, blob.rx);
        const perp = (-dx * sa + dy * ca) / Math.max(1, blob.ry);
        const d2 = along * along + perp * perp;
        const rugged = 0.82 + valueNoise(wx * 0.18 + blob.x, wy * 0.18 - blob.y, seed + 12600, 8) * 0.42;
        mountains += Math.exp(-d2 * 2.55) * blob.h * rugged;
      }

      const px = wx / (MAP_W - 1);
      const py = wy / (MAP_H - 1);
      let spineRidges = 0;
      for (let si = 0; si < spines.length; si++) spineRidges += jaggedRidgeContribution(px, py, spines[si], seed);
      let branchRidges = 0;
      for (const ridge of branches) branchRidges += jaggedRidgeContribution(px, py, ridge, seed);
      let detachedRidges = 0;
      for (const ridge of detachedRanges) detachedRidges += jaggedRidgeContribution(px, py, ridge, seed);
      let alpineMassifs = 0;
      for (const peak of alpinePeaks) {
        const dx = px - peak.x;
        const dy = py - peak.y;
        const ca = Math.cos(peak.angle);
        const sa = Math.sin(peak.angle);
        const along = (dx * ca + dy * sa) / Math.max(0.002, peak.rx);
        const perp = (-dx * sa + dy * ca) / Math.max(0.002, peak.ry);
        const d2 = along * along + perp * perp;
        const crag = 0.78 + valueNoise(px * 38 + peak.seedOffset, py * 38 - peak.seedOffset, seed + 12700, 5) * 0.52;
        alpineMassifs += Math.exp(-d2 * 1.85) * peak.h * crag;
      }
      let rangeBreakField = 0;
      for (const feature of rangeBreaks) rangeBreakField += elongatedFeatureContribution(px, py, feature, seed);
      const ridges = Math.max(0, spineRidges + branchRidges + detachedRidges * 0.95 + alpineMassifs * 0.70 - rangeBreakField * 0.90);

      const coast = coastPressureAt(x, y, wx, wy);
      const coastLower = coast.pressure;
      const folded = foldedOrogenyAt(px, py, seed, terrainTemplate, coastLower);
      const orogenicUplift = folded.uplift;
      const orogenicRidges = folded.ridges;
      const orogenicValleys = folded.valleys;
      // v5 terrain lifecycle:
      // 1) 全セルを海面下の海底として置く。
      // 2) 造山帯・広域隆起・尾根核を海底から持ち上げる。
      // 3) この後の流路計算で谷・扇状地・沖積平野を作る。
      const terrainLarge = fbm(wx * 0.30 + 40, wy * 0.30 - 60, seed + 710);
      const terrainRegional = fbm(wx * 0.72 + 80, wy * 0.72 - 20, seed + 777);
      const terrainLocal = fbm(wx * 1.55 + 17, wy * 1.55 - 31, seed + 1777);
      const terrainFine = valueNoise(wx * 2.15 + 11, wy * 2.15 - 19, seed + 2444, 5.5);
      const fineDissection = (Math.abs(terrainLocal - 0.5) * 0.050 + Math.abs(terrainFine - 0.5) * 0.024) * (0.62 + terrainTemplate.roughness * 0.58);
      const edgeDistance = Math.min(x, y, MAP_W - 1 - x, MAP_H - 1 - y) / Math.min(MAP_W, MAP_H);
      const deepEdge = 1 - smoothstep(edgeDistance / 0.15);
      const backboneDX = px - terrainTemplate.backboneCenterX;
      const backboneDY = py - terrainTemplate.backboneCenterY;
      const backboneAlong = (backboneDX * Math.cos(terrainTemplate.spineAngle) + backboneDY * Math.sin(terrainTemplate.spineAngle)) / Math.max(0.001, terrainTemplate.backboneLength);
      const backboneCross = (-backboneDX * Math.sin(terrainTemplate.spineAngle) + backboneDY * Math.cos(terrainTemplate.spineAngle)) / Math.max(0.001, terrainTemplate.backboneWidth);
      const backboneCore = clamp(1 - Math.sqrt(backboneAlong * backboneAlong + backboneCross * backboneCross));
      const basin = 0.040 * (terrainRegional - 0.5) + 0.030 * (terrainLarge - 0.5) + (backboneCore - 0.5) * 0.010;
      const protoHighland = clamp(orogenicUplift * 1.20 + orogenicRidges * 0.72 + spineRidges * 0.98 + branchRidges * 0.76 + detachedRidges * 0.66 + alpineMassifs * 0.72 + mountains * 0.24 - rangeBreakField * 1.35 - orogenicValleys * 0.42 + backboneCore * 0.10);
      // 海域は画面端の一律沈降ではなく、海岸圧・低地性・非山地性から開く。
      const marineOpening = clamp(
        coastLower * (0.95 + terrainTemplate.deposition * 0.35) +
        Math.max(0, 0.42 - protoHighland) * 0.22 +
        Math.max(0, -basin) * 0.20 -
        backboneCore * 0.10
      );
      const protoLowland = clamp((1 - protoHighland) * 0.40 + marineOpening * 0.28 + Math.max(0, -basin) * 0.34 + orogenicValleys * 0.34);
      const plainNoiseSuppression = protoLowland * terrainTemplate.plainNoiseSuppression;
      const subduedTerrainLocal = lerp(terrainLocal, 0.5, plainNoiseSuppression * 0.70);
      const subduedTerrainFine = lerp(terrainFine, 0.5, plainNoiseSuppression * 0.82);
      const subduedDissection = fineDissection * (1 - plainNoiseSuppression * 0.86);
      const seafloor = seaLevel -
        (0.188 + (1 - terrainLarge) * 0.045 + (1 - terrainRegional) * 0.034 + coastLower * (0.070 + terrainTemplate.deposition * 0.035) + marineOpening * 0.045 + terrainTemplate.deposition * 0.016);
      const platformEmergence = (0.050 + terrainTemplate.orographicCoverage * 0.032) * (1 - marineOpening * 0.62);
      const broadEmergence = Math.pow(protoHighland, 1.02) * (0.305 + terrainTemplate.orographicCoverage * 0.130);
      const ridgeEmergence =
        orogenicUplift * 0.245 +
        orogenicRidges * 0.225 +
        spineRidges * 0.520 +
        branchRidges * 0.430 +
        detachedRidges * 0.310 +
        alpineMassifs * 0.360 +
        mountains * 0.080;
      const shelfDepression = Math.max(0, -basin) * 0.032 + marineOpening * (0.070 + terrainTemplate.deposition * 0.030);
      const rawElevation =
        seafloor +
        platformEmergence +
        broadEmergence +
        ridgeEmergence +
        (terrainLarge - 0.5) * 0.055 +
        (terrainRegional - 0.5) * 0.040 +
        (subduedTerrainLocal - 0.5) * 0.036 +
        (subduedTerrainFine - 0.5) * 0.014 +
        subduedDissection +
        basin -
        shelfDepression -
        rangeBreakField * (0.155 + terrainTemplate.erosion * 0.080) -
        orogenicValleys * (0.035 + terrainTemplate.erosion * 0.045);

      elevation[i] = clamp(softUpperClamp(rawElevation, terrainTemplate.peakSoftStart, terrainTemplate.peakSoftCap));
      arcSpineField[i] = clamp(orogenicRidges * 1.05 + orogenicUplift * 0.58 + spineRidges * 1.65 + detachedRidges * 0.98 + alpineMassifs * 0.88);
      branchRidgeField[i] = clamp(branchRidges * 1.70 + orogenicValleys * 0.18);
      ridgeField[i] = clamp(arcSpineField[i] * 0.78 + branchRidgeField[i] * 0.42 + orogenicRidges * 0.48 + orogenicUplift * 0.25 + subduedDissection * 0.85 - rangeBreakField * 0.72 - orogenicValleys * 0.36);
      basinField[i] = clamp(Math.max(0, -basin) * 2.6 + rangeBreakField * 1.25 + orogenicValleys * 1.00 + Math.max(0, seaLevel + 0.085 - elevation[i]) * (1.10 + terrainTemplate.deposition * 0.65));
      coastalLowland[i] = 0;
      moisture[i] = clamp(0.44 * fbm(wx + 400, wy - 200, seed + 333) + 0.18 * valueNoise(wx, wy, seed + 343, 11) + 0.22 * (1 - Math.abs(ny * 1.7)) + 0.12 * deepEdge + 0.18 * (1 - clamp((elevation[i] - seaLevel) / 0.35)) - Math.max(0, elevation[i] - 0.62) * 0.22);
    }
  }

  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      // v5: 海岸線は後から切るのではなく、海底からの隆起量が海面を超えた場所だけを陸にする。
      if (elevation[i] < seaLevel) sea[i] = 1;
      if (sea[i]) elevation[i] = Math.min(elevation[i], seaLevel - 0.014 + hash2(x, y, seed + 2311) * 0.010);
    }
  }

  // Edge-connected water is ocean. Isolated water is only kept when it reads as
  // a small mountain/valley lake or lagoon; oversized round basins become wet lowland.
  const waterSeen = new Uint8Array(SIZE);
  const oceanQueue = [];
  for (let x = 0; x < MAP_W; x++) {
    for (const y of [0, MAP_H - 1]) {
      const i = indexOf(x, y);
      if (sea[i] && !waterSeen[i]) {
        waterSeen[i] = 1;
        ocean[i] = 1;
        oceanQueue.push(i);
      }
    }
  }
  for (let y = 0; y < MAP_H; y++) {
    for (const x of [0, MAP_W - 1]) {
      const i = indexOf(x, y);
      if (sea[i] && !waterSeen[i]) {
        waterSeen[i] = 1;
        ocean[i] = 1;
        oceanQueue.push(i);
      }
    }
  }
  for (let q = 0; q < oceanQueue.length; q++) {
    const cur = oceanQueue[q];
    const [x, y] = [cur % MAP_W, Math.floor(cur / MAP_W)];
    for (const [nx, ny] of neighbors8(x, y)) {
      const ni = indexOf(nx, ny);
      if (!sea[ni] || waterSeen[ni]) continue;
      waterSeen[ni] = 1;
      ocean[ni] = 1;
      oceanQueue.push(ni);
    }
  }
  for (let i = 0; i < SIZE; i++) {
    if (!sea[i] || waterSeen[i]) continue;
    const queue = [i];
    const component = [i];
    waterSeen[i] = 1;
    let sx = 0, sy = 0, perimeter = 0, ridgeSum = 0, valleySum = 0, coastTouch = 0;
    for (let q = 0; q < queue.length; q++) {
      const cur = queue[q];
      const x = cur % MAP_W;
      const y = Math.floor(cur / MAP_W);
      sx += x;
      sy += y;
      ridgeSum += ridgeField[cur];
      valleySum += valleyField[cur];
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (!sea[ni]) {
          perimeter++;
          if (coastalLowland[ni] > 0.12 || coastPressureAt(nx, ny).pressure > 0.42) coastTouch++;
          continue;
        }
        if (waterSeen[ni]) continue;
        waterSeen[ni] = 1;
        queue.push(ni);
        component.push(ni);
      }
    }
    const area = component.length;
    const cx = sx / area;
    const cy = sy / area;
    let radiusSum = 0;
    for (const ci of component) {
      const x = ci % MAP_W;
      const y = Math.floor(ci / MAP_W);
      radiusSum += Math.hypot(x - cx, y - cy);
    }
    const meanRadius = radiusSum / Math.max(1, area);
    const circularity = perimeter > 0 ? (4 * Math.PI * area) / (perimeter * perimeter) : 1;
    const mountainLake = area <= 38 && ridgeSum / area > 0.28;
    const valleyLake = area <= 70 && valleySum / area > 0.24 && circularity < 0.58;
    const lagoon = area <= 110 && coastTouch / Math.max(1, perimeter) > 0.18 && circularity < 0.70;
    const rareSpecial = area <= 145 && circularity < 0.52 && hash2(Math.round(cx), Math.round(cy), seed + 2401) > 0.88;
    const keepLake = mountainLake || valleyLake || lagoon || rareSpecial;
    for (const ci of component) {
      if (keepLake) {
        lake[ci] = 1;
        continue;
      }
      sea[ci] = 0;
      elevation[ci] = Math.max(seaLevel + 0.012, seaLevel + Math.min(0.055, meanRadius * 0.004) + hash2(ci, area, seed + 2402) * 0.012);
      basinField[ci] = clamp(basinField[ci] + 0.42);
      valleyField[ci] = clamp(valleyField[ci] + 0.18);
      depositionalLowland[ci] = clamp(depositionalLowland[ci] + 0.28);
      depositionField[ci] = clamp(depositionField[ci] + 0.035);
    }
  }

  // Seed ごとの過剰な平坦化・過剰な高原化を抑える救済正規化。
  // 海面を超えた陸だけを対象に、90/97 パーセンタイルを「山がちな島弧」の範囲へ寄せる。
  const emergedElevations = [];
  for (let i = 0; i < SIZE; i++) if (!sea[i]) emergedElevations.push(elevation[i]);
  emergedElevations.sort((a, b) => a - b);
  if (emergedElevations.length > 100) {
    const q = (p) => emergedElevations[Math.max(0, Math.min(emergedElevations.length - 1, Math.floor((emergedElevations.length - 1) * p)))];
    const q65 = q(0.65);
    const q90 = q(0.90);
    const q97 = q(0.97);
    const targetQ90 = 0.680 + terrainTemplate.roughness * 0.070;
    const targetQ97 = 0.840 + terrainTemplate.roughness * 0.095;
    const scale97 = clamp((targetQ97 - seaLevel) / Math.max(0.045, q97 - seaLevel), 0.52, 1.50);
    const shift90 = clamp(targetQ90 - q90, -0.30, 0.22);
    for (let i = 0; i < SIZE; i++) {
      if (sea[i]) continue;
      const highShoulder = smoothstep((elevation[i] - q65) / Math.max(0.045, q97 - q65));
      const ridgeBoost = clamp(ridgeField[i] * 0.42 + arcSpineField[i] * 0.28 - valleyField[i] * 0.18 - coastalLowland[i] * 0.15);
      let adjusted = seaLevel + (elevation[i] - seaLevel) * lerp(1, scale97, highShoulder * 0.92 + ridgeBoost * 0.22);
      adjusted += shift90 * highShoulder * (0.48 + ridgeBoost * 0.30);
      elevation[i] = clamp(softUpperClamp(adjusted, 0.895, targetQ97 + 0.070), seaLevel + 0.006, 0.992);
    }
  }

  // 海岸は「切断面」ではなく沈降・波食・堆積で丸める。
  // 強い尾根が海へ落ちる場所は残し、低い場所だけを浜堤・海岸低地へ寄せる。
  for (let y = 0; y < MAP_H; y++) {
    for (let x = 0; x < MAP_W; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let nearestSea = INF;
      let nearestOcean = INF;
      for (let dy = -8; dy <= 8; dy++) {
        for (let dx = -8; dx <= 8; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny) || !sea[indexOf(nx, ny)]) continue;
          nearestSea = Math.min(nearestSea, Math.hypot(dx, dy));
          if (ocean[indexOf(nx, ny)]) nearestOcean = Math.min(nearestOcean, Math.hypot(dx, dy));
        }
      }
      if (nearestSea <= 8) {
        const marineInfluence = smoothstep((8 - nearestSea) / 8);
        const ridgeResistance = smoothstep((ridgeField[i] - 0.24) / 0.42);
        const coastalShelf = seaLevel + 0.015 + nearestSea * (0.015 + terrainTemplate.deposition * 0.008) + Math.max(0, fbm(x * 1.1, y * 1.1, seed + 2350) - 0.5) * (0.010 + terrainTemplate.coastRoughness * 0.012);
        const lowlandBlend = marineInfluence * (1 - ridgeResistance) * (0.42 + terrainTemplate.deposition * 0.30);
        if (elevation[i] > coastalShelf) elevation[i] = lerp(elevation[i], coastalShelf, lowlandBlend);
        if (nearestOcean <= 8) {
          const plainReach = clamp(5.2 + terrainTemplate.deposition * 7.0, 5.5, 10.5);
          coastalLowland[i] = clamp((1 - nearestOcean / plainReach) * (0.70 + terrainTemplate.deposition * 0.55) * (1 - ridgeField[i] * 0.62));
        }
      }
    }
  }

  // Explicit alpine punctuation.  The base ridge system defines broad relief,
  // while these narrow, detached high points make several visually legible
  // mountain groups instead of one round central mass.
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i] || coastalLowland[i] > 0.42) continue;
      const px = x / (MAP_W - 1);
      const py = y / (MAP_H - 1);
      let peakSignal = 0;
      for (const peak of alpinePeaks) {
        const dx = px - peak.x;
        const dy = py - peak.y;
        const ca = Math.cos(peak.angle);
        const sa = Math.sin(peak.angle);
        const along = (dx * ca + dy * sa) / Math.max(0.002, peak.rx);
        const perp = (-dx * sa + dy * ca) / Math.max(0.002, peak.ry);
        const d2 = along * along + perp * perp;
        peakSignal += Math.exp(-d2 * 2.20) * peak.h;
      }
      if (peakSignal <= 0.026) continue;
      const crag = Math.max(0, valueNoise(x * 2.4 + 73, y * 2.4 - 91, seed + 12880, 3.5) - 0.36);
      const target = clamp(0.64 + peakSignal * 2.45 + crag * 0.092, seaLevel + 0.006, 0.992);
      elevation[i] = Math.max(elevation[i], target);
      ridgeField[i] = clamp(ridgeField[i] + peakSignal * 4.6 + crag * 0.28);
      arcSpineField[i] = clamp(arcSpineField[i] + peakSignal * 3.2);
      basinField[i] = Math.max(0, basinField[i] - peakSignal * 1.2);
      depositionalLowland[i] = Math.max(0, depositionalLowland[i] - peakSignal * 1.5);
    }
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[indexOf(x, y)] = clamp(Math.sqrt(gx * gx + gy * gy) * 10.5);
    }
  }

  const landOrder = [];
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let low = i;
      let best = elevation[i] + 0.012 * hash2(x, y, seed + 2468);
      let localMean = 0;
      let localMax = elevation[i];
      let localMin = elevation[i];
      let nCount = 0;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const ev = elevation[ni];
        localMean += ev;
        localMax = Math.max(localMax, ev);
        localMin = Math.min(localMin, ev);
        nCount++;
        const directed = ev + 0.008 * hash2(nx, ny, seed + 2469);
        if (directed < best || sea[ni]) {
          best = directed;
          low = ni;
        }
      }
      if (low !== i) flowTo[i] = low;
      localMean /= Math.max(1, nCount);
      const hollow = Math.max(0, localMean - elevation[i]);
      const relief = localMax - localMin;
      valleyField[i] = clamp(hollow * 8.4 + Math.max(0, 0.42 - elevation[i]) * 0.32 + moisture[i] * 0.08 - ridgeField[i] * 0.18);
      basinField[i] = clamp(basinField[i] + hollow * 2.4 + (relief < 0.055 && elevation[i] < 0.55 ? 0.18 : 0));
      flowAccum[i] = 0.82 + moisture[i] * 0.88 + valleyField[i] * 0.78 + Math.max(0, elevation[i] - seaLevel) * 0.14;
      landOrder.push(i);
    }
  }
  landOrder.sort((a, b) => elevation[b] - elevation[a]);
  for (const i of landOrder) {
    const to = flowTo[i];
    if (to >= 0 && to !== i) flowAccum[to] += flowAccum[i] * 0.91;
  }
  let maxFlowAccum = 0;
  for (let i = 0; i < SIZE; i++) if (!sea[i]) maxFlowAccum = Math.max(maxFlowAccum, flowAccum[i]);
  if (maxFlowAccum > 0) {
    for (let i = 0; i < SIZE; i++) flowAccum[i] = clamp(flowAccum[i] / maxFlowAccum);
  }
  for (let i = 0; i < SIZE; i++) {
    if (!sea[i]) valleyField[i] = clamp(valleyField[i] * 0.62 + Math.pow(flowAccum[i], 0.48) * 0.62);
  }

  function recomputeDrainageFields({ reinforceValleys = false } = {}) {
    flowTo.fill(-1);
    flowAccum.fill(0);
    const order = [];
    for (let y = 1; y < MAP_H - 1; y++) {
      for (let x = 1; x < MAP_W - 1; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        let low = i;
        let best = elevation[i] + 0.010 * hash2(x, y, seed + 18440);
        let localMean = 0;
        let localMax = elevation[i];
        let localMin = elevation[i];
        let nCount = 0;
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          const ev = elevation[ni];
          localMean += ev;
          localMax = Math.max(localMax, ev);
          localMin = Math.min(localMin, ev);
          nCount++;
          const drainageAttraction = valleyField[ni] * 0.020 + erosionField[ni] * 0.030 + gullyIncisionField?.[ni] * 0.020;
          const directed = ev + 0.007 * hash2(nx, ny, seed + 18441) - drainageAttraction;
          if (directed < best || sea[ni]) {
            best = directed;
            low = ni;
          }
        }
        if (low !== i) flowTo[i] = low;
        localMean /= Math.max(1, nCount);
        const hollow = Math.max(0, localMean - elevation[i]);
        const relief = localMax - localMin;
        if (reinforceValleys) {
          valleyField[i] = clamp(valleyField[i] * 0.72 + hollow * 7.8 + Math.max(0, 0.055 - relief) * Math.max(0, elevation[i] - seaLevel - 0.10) * 1.2);
          basinField[i] = clamp(basinField[i] + hollow * 1.8 + (relief < 0.040 && elevation[i] < 0.58 ? 0.10 : 0));
        }
        flowAccum[i] = 0.72 + moisture[i] * 0.72 + valleyField[i] * 0.72 + Math.max(0, elevation[i] - seaLevel) * 0.10;
        order.push(i);
      }
    }
    order.sort((a, b) => elevation[b] - elevation[a]);
    for (const i of order) {
      const to = flowTo[i];
      if (to >= 0 && to !== i) flowAccum[to] += flowAccum[i] * 0.93;
    }
    let maxAccum = 0;
    for (let i = 0; i < SIZE; i++) if (!sea[i]) maxAccum = Math.max(maxAccum, flowAccum[i]);
    if (maxAccum > 0) {
      for (let i = 0; i < SIZE; i++) if (!sea[i]) flowAccum[i] = clamp(flowAccum[i] / maxAccum);
    }
  }

  // v6/v7: 可視河川だけでなく、山地の無数の沢・ガリーを先に掘る。
  // これを描画用の川には使わず、地形侵食だけに使うことで「日本的な谷密度」を出す。
  function traceErosionGully(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    let lastDx = 0;
    let lastDy = 0;
    const path = [];
    const seen = new Set();
    for (let step = 0; step < 92; step++) {
      const i = indexOf(x, y);
      if (seen.has(i) || sea[i]) break;
      seen.add(i);
      path.push([x, y]);
      if (elevation[i] < seaLevel + 0.055 && path.length > 8) break;
      let best = null;
      let bestValue = INF;
      const preferred = flowTo[i];
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (seen.has(ni)) continue;
        const dx = nx - x;
        const dy = ny - y;
        const drop = elevation[i] - elevation[ni];
        const uphill = Math.max(0, -drop);
        if (!sea[ni] && uphill > 0.035 + flowAccum[i] * 0.10) continue;
        const sameDirection = lastDx || lastDy ? (dx * lastDx + dy * lastDy) / Math.max(0.001, Math.hypot(dx, dy) * Math.hypot(lastDx, lastDy)) : 0;
        const bend = Math.abs(dx * lastDy - dy * lastDx);
        const preferredBonus = ni === preferred ? 0.80 : 0;
        const value =
          elevation[ni] * 0.96 +
          uphill * 30.0 -
          Math.max(0, drop) * 2.9 -
          valleyField[ni] * 1.15 -
          Math.pow(flowAccum[ni], 0.55) * 1.05 -
          moisture[ni] * 0.22 -
          preferredBonus -
          Math.max(0, sameDirection) * 0.045 +
          bend * 0.030 +
          (hash2(nx, ny, seed + bonusSeed + step * 31) - 0.5) * 0.055;
        if (value < bestValue) {
          bestValue = value;
          best = [nx, ny, dx, dy];
        }
      }
      if (!best) break;
      x = best[0];
      y = best[1];
      lastDx = best[2];
      lastDy = best[3];
    }
    return path;
  }

  const gullyIncisionField = new Float32Array(SIZE);
  const visibleRavineField = new Float32Array(SIZE);
  const surfaceTextureField = new Float32Array(SIZE);
  const gullyCandidates = [];
  for (let y = 3; y < MAP_H - 3; y++) {
    for (let x = 3; x < MAP_W - 3; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const highland = clamp((elevation[i] - (seaLevel + 0.075)) / 0.42);
      const highFlat = highland * clamp(1 - slope[i] * 2.7) * clamp((elevation[i] - (seaLevel + 0.170)) / 0.34) * (1 - coastalLowland[i] * 0.70);
      const relief = clamp(slope[i] * 1.10 + ridgeField[i] * 0.34 + valleyField[i] * 0.22 + highFlat * 0.86);
      const drainage = clamp(Math.pow(flowAccum[i], 0.36) * 0.74 + moisture[i] * 0.26 + dendriticRavineTexture(x, y, seed) * 0.18 + highFlat * 0.58);
      const score = highland * relief * drainage + highFlat * 0.44 + hash2(x, y, seed + 15200) * 0.10;
      if (score > 0.13) gullyCandidates.push({ x, y, score });
    }
  }
  // 明示的な細流路大量生成は粗い格子では効きにくく、過剰掘削の原因にもなる。
  // 細かな山肌表現は renderer 側の手続きノイズへ移し、ここでは実際の小河川は掘らない。
  const gullySources = [];
  const erosionGullyPaths = [];
  const valleyNetworkIncision = new Float32Array(SIZE);
  const valleyFloorField = new Float32Array(SIZE);
  for (const source of gullySources) {
    const path = traceErosionGully(source.x, source.y, 15250 + source.x * 17 + source.y * 23);
    if (path.length < 5) continue;
    erosionGullyPaths.push(path);
    const startIndex = indexOf(source.x, source.y);
    const highFlat = clamp(1 - slope[startIndex] * 2.4) * clamp((elevation[startIndex] - (seaLevel + 0.16)) / 0.35);
    const strength = terrainTemplate.gullyIncision * clamp(0.016 + slope[startIndex] * 0.026 + Math.pow(flowAccum[startIndex], 0.42) * 0.020 + valleyField[startIndex] * 0.010 + highFlat * 0.018);
    addPathIncision(gullyIncisionField, path, strength, path.length > 22 ? 2 : 1);
    const valleyRadius = path.length > 54 ? 3.2 : path.length > 28 ? 2.4 : 1.7;
    addValleyDistanceInfluence(valleyNetworkIncision, valleyFloorField, path, strength * (0.55 + highFlat * 0.55), valleyRadius, path.length > 42 ? 0.80 : 0.48);
    addPathIncision(visibleRavineField, path, strength * (8.4 + highFlat * 3.4), path.length > 34 ? 2 : 1);
    addValleyDistanceInfluence(visibleRavineField, valleyFloorField, path, strength * (2.2 + highFlat * 1.8), Math.max(1.5, valleyRadius * 0.72), 0.30);
  }


  // v10 の明示的な支谷トレースは削除。粗い解像度では可視変化が薄く、
  // 処理だけが増えるため、支谷表現は後段の手続きノイズに集約する。

  // v11: red-box style dissected relief should not be a local accident.
  // Spread that branchy mountain texture across the whole mountainous massif,
  // while still letting real drainage control the strongest valleys.
  const massifDissectionField = new Float32Array(SIZE);
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const highland = clamp((elevation[i] - (seaLevel + 0.055)) / 0.46);
      const mountainMask = highland * clamp(slope[i] * 1.35 + ridgeField[i] * 0.42 + arcSpineField[i] * 0.24 + branchRidgeField[i] * 0.18 - basinField[i] * 0.22) * (1 - coastalLowland[i] * 0.70);
      if (mountainMask <= 0.04) continue;
      const highFlat = clamp((elevation[i] - (seaLevel + 0.16)) / 0.34) * clamp(1 - slope[i] * 2.7) * mountainMask;
      const warpX = x + (fbm(x * 0.075 + 17, y * 0.075 - 29, seed + 16201) - 0.5) * 7.8;
      const warpY = y + (fbm(x * 0.075 - 53, y * 0.075 + 11, seed + 16202) - 0.5) * 7.8;
      const ravA = dendriticRavineTexture(warpX * 0.70, warpY * 0.70, seed + 16203);
      const ravB = dendriticRavineTexture(warpX * 1.02 + 37, warpY * 1.02 - 19, seed + 16204);
      const coarse = Math.abs((fbm(warpX * 0.085 + 21, warpY * 0.085 - 8, seed + 16205) - 0.5) * 2);
      const patch = clamp(0.64 + (fbm(x * 0.028 + 80, y * 0.028 - 41, seed + 16206) - 0.5) * 0.78);
      const branchiness = clamp(ravA * 0.64 + ravB * 0.42 + coarse * 0.18);
      const dissection = mountainMask * patch * (branchiness * (1.00 + highFlat * 0.72) + highFlat * 0.26);
      if (dissection <= 0.03) continue;
      massifDissectionField[i] = dissection;
      visibleRavineField[i] = clamp(visibleRavineField[i] + dissection * (0.34 + highFlat * 0.18));
      surfaceTextureField[i] = clamp(surfaceTextureField[i] + dissection * (0.40 + highFlat * 0.32) + coarse * mountainMask * 0.08);
      valleyField[i] = clamp(valleyField[i] + dissection * 0.22);
    }
  }


  // Coarse grids cannot resolve endless tiny tributaries.  Here the “fine
  // ravines” are not just a paint/shading overlay: they are folded into the DEM
  // itself as a small alternating cut/crest signal, then slope is recomputed.
  // This keeps the detail persistent for rivers, labels and all later terrain uses.
  const proceduralRavineField = new Float32Array(SIZE);
  const proceduralReliefField = new Float32Array(SIZE);
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const highland = clamp((elevation[i] - (seaLevel + 0.07)) / 0.42);
      const mountainMask = highland * clamp(slope[i] * 0.84 + ridgeField[i] * 0.48 + arcSpineField[i] * 0.22 - basinField[i] * 0.20) * (1 - coastalLowland[i] * 0.78);
      if (mountainMask <= 0.025) continue;
      const warpX = x + (fbm(x * 0.10 + 91, y * 0.10 - 47, seed + 16601) - 0.5) * 7.8;
      const warpY = y + (fbm(x * 0.10 - 33, y * 0.10 + 28, seed + 16602) - 0.5) * 7.8;
      const ravA = dendriticRavineTexture(warpX * 0.74, warpY * 0.74, seed + 16603);
      const ravB = dendriticRavineTexture(warpX * 1.18 + 23, warpY * 1.18 - 41, seed + 16604);
      const ridged = 1 - Math.abs((valueNoise(warpX * 0.92 + 17, warpY * 0.92 - 9, seed + 16605, 6.2) - 0.5) * 2);
      const patch = clamp(0.66 + (fbm(x * 0.038 + 80, y * 0.038 - 51, seed + 16606) - 0.5) * 0.88);
      const micro = clamp((ravA * 0.70 + ravB * 0.46 + ridged * 0.22) * patch * mountainMask);
      if (micro <= 0.025) continue;

      const crestNoise = Math.max(0, ridged - 0.45) * mountainMask * patch;
      const cut = micro * (0.024 + highland * 0.024 + terrainTemplate.roughness * 0.012 + terrainTemplate.macroNoiseStrength * 0.20);
      const crest = crestNoise * (0.008 + highland * 0.010 + terrainTemplate.globalNoiseStrength * 0.18);
      const relief = crest - cut;
      proceduralRavineField[i] = micro;
      proceduralReliefField[i] = relief;
      elevation[i] = clamp(elevation[i] + relief, seaLevel + 0.006, 0.998);
      visibleRavineField[i] = clamp(visibleRavineField[i] + micro * 0.66);
      surfaceTextureField[i] = clamp(surfaceTextureField[i] + micro * 0.78 + Math.max(0, crestNoise) * 0.20);
      valleyField[i] = clamp(valleyField[i] + micro * 0.20);
      erosionField[i] = clamp(erosionField[i] + cut * 0.40);
    }
  }
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const upland = clamp((elevation[i] - seaLevel) / 0.56);
      const macroA = (fbm(x * terrainTemplate.macroNoiseScale + 301, y * terrainTemplate.macroNoiseScale - 119, seed + 16701) - 0.5) * 2;
      const macroB = (valueNoise(x * (terrainTemplate.macroNoiseScale * 1.9) - 87, y * (terrainTemplate.macroNoiseScale * 1.9) + 63, seed + 16702, 10.5) - 0.5) * 2;
      const globalA = (fbm(x * 0.016 + 57, y * 0.016 - 21, seed + 16703) - 0.5) * 2;
      const globalB = (valueNoise(x * 0.030 + 19, y * 0.030 - 44, seed + 16704, 14.0) - 0.5) * 2;
      const macroMask = clamp(0.30 + upland * 0.78 - coastalLowland[i] * 0.44);
      const macroRelief = (macroA * terrainTemplate.macroNoiseStrength + macroB * terrainTemplate.macroNoiseStrength * 0.62 + globalA * terrainTemplate.globalNoiseStrength + globalB * terrainTemplate.globalNoiseStrength * 0.58) * macroMask;
      elevation[i] = clamp(elevation[i] + macroRelief, seaLevel + 0.006, 0.998);
      visibleRavineField[i] = clamp(visibleRavineField[i] + Math.abs(macroRelief) * 6.4 * upland);
      surfaceTextureField[i] = clamp(surfaceTextureField[i] + Math.abs(macroRelief) * 9.0 * macroMask);
      valleyField[i] = clamp(valleyField[i] + Math.max(0, -macroRelief) * 2.4);
      ridgeField[i] = clamp(ridgeField[i] + Math.max(0, macroRelief) * 1.9);
    }
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[i] = clamp(Math.sqrt(gx * gx + gy * gy) * 10.5);
    }
  }

  // First-order fluvial shaping: cut valley floors on steep/high-flow cells and
  // deposit gently in coastal lowlands and basin floors. This gives visible
  // river valleys without destroying the macro terrain structure.
  const shapedElevation = new Float32Array(elevation);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const flow = Math.pow(flowAccum[i], 0.58);
      const incisionNoise = 0.82 + hash2(x, y, seed + 8120) * 0.36;
      const firstOrderPower = smoothstep((flowAccum[i] - 0.010) / 0.095);
      const highlandMask = clamp((elevation[i] - (seaLevel + 0.060)) / 0.44);
      const ravineTexture = dendriticRavineTexture(x, y, seed);
      const highFlat = highlandMask * clamp(1 - slope[i] * 2.9) * clamp((elevation[i] - (seaLevel + 0.170)) / 0.34) * (1 - coastalLowland[i] * 0.72);
      const basinProtection = clamp(basinField[i] * (1 - slope[i] * 2.2) * (1 - ridgeField[i] * 0.65));
      const textureMask = clamp(highlandMask * (slope[i] * 1.40 + ridgeField[i] * 0.32 + valleyField[i] * 0.26 + highFlat * 0.82) * (0.40 + moisture[i] * 0.78));
      const syntheticRavine = proceduralRavineField[i];
      const textureCut = clamp(terrainTemplate.dendriticTexture * ravineTexture * textureMask * (0.0016 + terrainTemplate.erosion * 0.0022 + slope[i] * 0.006 + highFlat * 0.003));
      const gullyCut = clamp(gullyIncisionField[i] * (0.34 + slope[i] * 0.76 + highlandMask * 0.24 + highFlat * 0.26));
      const networkCut = clamp(valleyNetworkIncision[i] * (0.18 + highlandMask * 0.22 + slope[i] * 0.24));
      const massifDissection = massifDissectionField[i];
      const visualRavine = clamp(visibleRavineField[i] * (0.16 + highlandMask * 0.32 + slope[i] * 0.18) + massifDissection * (0.10 + highFlat * 0.08) + syntheticRavine * 0.26);
      const plateauCut = clamp(highFlat * (ravineTexture * 0.002 + valleyNetworkIncision[i] * 0.06 + Math.pow(flowAccum[i], 0.52) * 0.0018 + massifDissection * 0.004 + syntheticRavine * 0.003) + visualRavine * 0.004);
      const sourceProtection = clamp((0.060 - flowAccum[i]) / 0.060) * clamp((elevation[i] - (seaLevel + 0.08)) / 0.40);
      const fluvialCore = terrainTemplate.fluvialAggression * firstOrderPower * flow * (0.018 + terrainTemplate.erosion * 0.028 + slope[i] * (0.050 + terrainTemplate.erosion * 0.050) + ridgeField[i] * (0.008 + terrainTemplate.erosion * 0.014)) * incisionNoise;
      let steepValley = clamp((fluvialCore + textureCut + gullyCut + networkCut + plateauCut) * (1 - basinProtection * 0.70));
      let lateralCut = clamp((terrainTemplate.fluvialAggression * firstOrderPower * Math.pow(flowAccum[i], 0.76) * valleyField[i] * (0.014 + terrainTemplate.erosion * 0.024) + valleyFloorField[i] * (0.050 + highlandMask * 0.08) + gullyIncisionField[i] * 0.06) * (1 - basinProtection * 0.82));
      if (sourceProtection > 0) {
        const protect = 1 - sourceProtection * 0.82;
        steepValley *= protect;
        lateralCut *= 1 - sourceProtection * 0.74;
      }
      const lowSettling = clamp(flow * (coastalLowland[i] * (0.018 + terrainTemplate.deposition * 0.040) + basinField[i] * (0.018 + terrainTemplate.deposition * 0.040) + (elevation[i] < 0.40 ? 0.006 + terrainTemplate.deposition * 0.018 : 0)) * (1 - slope[i] * 0.82) * (1 - ridgeField[i] * 0.45) + basinProtection * (0.006 + terrainTemplate.deposition * 0.012));
      erosionField[i] = steepValley + lateralCut + networkCut * 0.20;
      depositionField[i] = lowSettling;
      depositionalLowland[i] = clamp(lowSettling * 6.5 + basinField[i] * terrainTemplate.deposition * 0.28 + coastalLowland[i] * terrainTemplate.deposition * 0.34 + valleyFloorField[i] * 0.18);
      valleyField[i] = clamp(valleyField[i] + ravineTexture * textureMask * 0.22 + gullyIncisionField[i] * 1.8 + valleyNetworkIncision[i] * 1.5 + visibleRavineField[i] * 1.30 + massifDissection * 0.80 + syntheticRavine * 0.42 + valleyFloorField[i] * 1.0 + steepValley * 1.4);
      surfaceTextureField[i] = clamp(surfaceTextureField[i] + ravineTexture * textureMask * 0.30 + visibleRavineField[i] * 0.66 + massifDissection * 0.72 + syntheticRavine * 0.74 + gullyIncisionField[i] * 1.2);
      const basinFloorGuard = seaLevel + 0.010 + basinProtection * 0.030 + coastalLowland[i] * 0.010;
      const headwaterGuard = elevation[i] - (0.010 + slope[i] * 0.014 + highlandMask * 0.010 + Math.pow(flowAccum[i], 0.60) * 0.028);
      const carved = elevation[i] - steepValley - lateralCut + lowSettling * (0.72 + basinProtection * 0.35);
      shapedElevation[i] = clamp(Math.max(basinFloorGuard, sourceProtection > 0 ? headwaterGuard : seaLevel + 0.006, carved), seaLevel + 0.006, 1);
    }
  }
  elevation.set(shapedElevation);

  // Final orographic pass: ensure true alpine/high-mountain cells remain after
  // river incision and lowland smoothing.  Uplift is confined to ridge cores and
  // fades out in valley floors so drainage still reads correctly.
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const dissectionLock = clamp(valleyNetworkIncision[i] * 2.6 + gullyIncisionField[i] * 2.4 + erosionField[i] * 2.0);
      const ridgeCore = clamp(arcSpineField[i] * 0.74 + branchRidgeField[i] * 0.58 + ridgeField[i] * 0.42 - valleyField[i] * 0.38 - flowAccum[i] * 0.28 - dissectionLock * 0.18);
      const highBase = clamp((elevation[i] - 0.55) / 0.25);
      const alpine = clamp(ridgeCore * 0.88 + highBase * 0.18 - coastalLowland[i] * 0.45 - depositionalLowland[i] * 0.36 - dissectionLock * 0.20);
      if (alpine <= 0.08) continue;
      const summitTexture = Math.max(0, valueNoise(x * 2.7 + 31, y * 2.7 - 41, seed + 9771, 3.0) - 0.38);
      const uplift = Math.pow(alpine, 1.55) * (0.032 + terrainTemplate.roughness * 0.040 + summitTexture * 0.032);
      const summitCap = 0.970 + Math.min(0.045, ridgeCore * 0.042) + summitTexture * 0.020;
      elevation[i] = clamp(elevation[i] + uplift, seaLevel + 0.006, summitCap);
      ridgeField[i] = clamp(ridgeField[i] + uplift * 1.15);
      erosionField[i] = Math.max(0, erosionField[i] - uplift * 0.25);
    }
  }

  applyAlpineMicroRelief(elevation, sea, slope, ridgeField, valleyField, coastalLowland, seaLevel, seed, terrainTemplate, surfaceTextureField);
  breakHighPlateaus(elevation, sea, slope, ridgeField, valleyField, flowAccum, coastalLowland, seaLevel, seed, terrainTemplate, visibleRavineField, surfaceTextureField);

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[i] = clamp(Math.sqrt(gx * gx + gy * gy) * 10.5);
      visibleRavineField[i] = clamp(visibleRavineField[i] + Math.pow(flowAccum[i], 0.56) * clamp(slope[i] * 1.5 + valleyField[i] * 0.28) * 0.24);
      valleyField[i] = clamp(valleyField[i] + visibleRavineField[i] * 0.42 + erosionField[i] * 2.1 + depositionField[i] * 0.8 - ridgeField[i] * 0.06);
      basinField[i] = clamp(basinField[i] + depositionField[i] * 1.6);
    }
  }

  // v7: after actual incision, recompute flow direction/accumulation so the
  // visible river network follows the carved valleys rather than the pre-erosion surface.
  recomputeDrainageFields({ reinforceValleys: true });

  // v8: breach the most prominent enclosed basins so盆地 often acquire an outlet
  // instead of remaining as unexplained closed depressions.
  const sinkCandidates = [];
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      if (flowTo[i] >= 0 && flowTo[i] !== i) continue;
      const basinScore = basinField[i] * 1.15 + valleyField[i] * 0.55 + flowAccum[i] * 0.38 + (lake[i] ? 0.16 : 0) - coastalLowland[i] * 0.30;
      if (basinScore > 0.24) sinkCandidates.push({ i, x, y, score: basinScore });
    }
  }
  sinkCandidates.sort((a, b) => b.score - a.score);
  const usedOutletStarts = new Set();
  for (const candidate of sinkCandidates.slice(0, terrainTemplate.basinOutletCount)) {
    const startIndex = candidate.i;
    const startKey = `${candidate.x},${candidate.y}`;
    if (usedOutletStarts.has(startKey)) continue;
    let targetIndex = -1;
    let bestScore = INF;
    for (let j = 0; j < SIZE; j++) {
      if (j === startIndex) continue;
      if (sea[j]) {
        const d = Math.hypot((j % MAP_W) - candidate.x, Math.floor(j / MAP_W) - candidate.y);
        const score = d * 0.48 - 7.5;
        if (score < bestScore) { bestScore = score; targetIndex = j; }
        continue;
      }
      const elevDelta = elevation[j] - elevation[startIndex];
      if (elevDelta > 0.050) continue;
      if (flowAccum[j] < Math.max(0.08, flowAccum[startIndex] + 0.010) && basinField[j] < 0.18) continue;
      const dx = (j % MAP_W) - candidate.x;
      const dy = Math.floor(j / MAP_W) - candidate.y;
      const d = Math.hypot(dx, dy);
      const score = d * 0.34 + Math.max(0, elevDelta) * 120 - flowAccum[j] * 18 - valleyField[j] * 4 - (lake[j] ? 1.5 : 0);
      if (score < bestScore) { bestScore = score; targetIndex = j; }
    }
    if (targetIndex < 0) continue;
    const outletPath = carveOutletChannel(elevation, sea, lake, river, valleyField, basinField, flowAccum, seaLevel, startIndex, targetIndex, seed, 22800 + candidate.x * 17 + candidate.y * 31);
    if (outletPath.length > 2) {
      for (const [px, py] of outletPath.slice(0, 6)) usedOutletStarts.add(`${px},${py}`);
    }
  }
  recomputeDrainageFields({ reinforceValleys: true });

  const sourceCandidates = [];
  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gullyHint = gullyIncisionField[i] * 1.8 + dendriticRavineTexture(x, y, seed) * 0.12;
      const score = elevation[i] * 0.22 + moisture[i] * 0.25 + ridgeField[i] * 0.030 + arcSpineField[i] * 0.030 + branchRidgeField[i] * 0.018 + flowAccum[i] * 1.12 + valleyField[i] * 0.58 + gullyHint * 0.42 + basinField[i] * 0.12 + coastalLowland[i] * 0.07 + hash2(x, y, seed + 9000) * 0.05;
      if (elevation[i] > 0.30 && elevation[i] < 0.94 && moisture[i] > 0.16 && (flowAccum[i] > 0.003 || valleyField[i] > 0.035 || gullyIncisionField[i] > 0.006 || slope[i] > 0.18) && ridgeField[i] < 0.98) sourceCandidates.push({ x, y, score });
    }
  }

  const majorRiverCandidates = sourceCandidates
    .filter((p) => {
      const i = indexOf(p.x, p.y);
      return elevation[i] > 0.24 && elevation[i] < 0.80 && moisture[i] > 0.18 && (flowAccum[i] > 0.070 || valleyField[i] > 0.18);
    })
    .map((p) => {
      const i = indexOf(p.x, p.y);
      const inland = Math.min(p.x, p.y, MAP_W - 1 - p.x, MAP_H - 1 - p.y) / Math.min(MAP_W, MAP_H);
      return {
        ...p,
        score: p.score + flowAccum[i] * 1.55 + valleyField[i] * 0.60 + inland * 0.55 - Math.abs(elevation[i] - 0.50) * 0.16,
      };
    });
  const majorRiverDesired = rand(seed, 9120) < 0.30 ? 0 : (rand(seed, 9121) < 0.24 ? 2 : 1);
  const majorSources = pickEntities(majorRiverCandidates, {
    max: majorRiverDesired,
    minDistance: 12,
    threshold: 0.54,
    seed: seed + 9122,
    jitter: 0.01,
  });

  const sources = pickEntities(sourceCandidates, {
    max: 48 + Math.floor(rand(seed, 910) * 30),
    minDistance: 4,
    threshold: 0.22 + rand(seed, 911) * 0.05,
    seed,
  });

  function nearestWaterGoal(from) {
    let bestSea = null;
    let bestScore = INF;
    for (let y = 0; y < MAP_H; y++) {
      for (let x = 0; x < MAP_W; x++) {
        const i = indexOf(x, y);
        if (!sea[i]) continue;
        const d = Math.hypot(x - from.x, y - from.y);
        const score = d - coastalLowland[indexOf(Math.max(0, Math.min(MAP_W - 1, from.x)), Math.max(0, Math.min(MAP_H - 1, from.y)))] * 2;
        if (score < bestScore) {
          bestScore = score;
          bestSea = { x, y };
        }
      }
    }
    return bestSea;
  }

  function riverRouteCost(x, y, cx, cy) {
    const i = indexOf(x, y);
    const ci = indexOf(cx, cy);
    if (sea[i]) return 0.18;
    const uphill = Math.max(0, elevation[i] - elevation[ci]);
    const downhill = Math.max(0, elevation[ci] - elevation[i]);
    if (!sea[i] && uphill > 0.070 && flowAccum[i] < flowAccum[ci] + 0.010) return INF;
    return Math.max(
      0.18,
      1 +
        uphill * 86 +
        slope[i] * 0.38 +
        elevation[i] * 0.42 -
        downhill * 2.1 -
        valleyField[i] * 1.24 -
        flowAccum[i] * 1.18 -
        moisture[i] * 0.22 -
        coastalLowland[i] * 0.36
    );
  }

  function forceRiverToWater(path) {
    if (!path.length) return path;
    const [ex, ey] = path[path.length - 1];
    if (sea[indexOf(ex, ey)]) return path;
    const goal = nearestWaterGoal({ x: ex, y: ey });
    if (!goal) return path;
    const startElevation = elevation[indexOf(ex, ey)];
    const tail = aStar({ x: ex, y: ey }, goal, (x, y, cx, cy) => {
      const i = indexOf(x, y);
      const ci = indexOf(cx, cy);
      if (!sea[i] && elevation[i] > Math.max(startElevation + 0.160, elevation[ci] + 0.090)) return INF;
      return riverRouteCost(x, y, cx, cy);
    });
    if (tail.length <= 2) return path;
    return path.concat(tail.slice(1));
  }

  function confluenceAnglePenalty(nx, ny, dx, dy, lengthSoFar) {
    if (lengthSoFar < 7 || river[indexOf(nx, ny)] < 0.24) return 0;
    let best = 0.16;
    const inLen = Math.hypot(dx, dy) || 1;
    for (const [rx, ry] of neighbors8(nx, ny)) {
      if (river[indexOf(rx, ry)] < 0.22) continue;
      const rdx = rx - nx;
      const rdy = ry - ny;
      const cos = clamp((dx * rdx + dy * rdy) / Math.max(0.001, inLen * Math.hypot(rdx, rdy)), -1, 1);
      const angle = Math.acos(cos);
      const shallow = angle < 0.45 ? 0.28 : 0;
      best = Math.min(best, Math.abs(angle - Math.PI * 0.62) * 0.045 + shallow);
    }
    return best;
  }

  function traceRiverPath(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    let lastDx = 0;
    let lastDy = 0;
    const path = [];
    const seen = new Set();
    let accum = 0;

    for (let step = 0; step < 600; step++) {
      const i = indexOf(x, y);
      if (seen.has(i)) break;
      seen.add(i);
      path.push([x, y]);
      river[i] += 0.64 + path.length / 128 + flowAccum[i] * 0.92;
      accum += river[i] + flowAccum[i];
      if (sea[i]) break;

      let best = null;
      let bestValue = INF;
      const currentElevation = elevation[i];
      const preferred = flowTo[i];

      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const dx = nx - x;
        const dy = ny - y;
        const drop = currentElevation - elevation[ni];
        const uphill = Math.max(0, -drop);
        if (!sea[ni] && uphill > 0.040 && flowAccum[ni] < flowAccum[i] + 0.020) continue;
        let surrounding = 0;
        let surroundingCount = 0;
        for (const [vx, vy] of neighbors8(nx, ny)) {
          surrounding += elevation[indexOf(vx, vy)];
          surroundingCount++;
        }
        const valley = Math.max(0, surrounding / Math.max(1, surroundingCount) - elevation[ni]);
        const sameDirection = lastDx || lastDy ? (dx * lastDx + dy * lastDy) / Math.max(0.001, Math.hypot(dx, dy) * Math.hypot(lastDx, lastDy)) : 0;
        const lowlandMeander = clamp((1 - slope[i] * 2.4) * (0.35 + valleyField[i] * 0.60 + Math.pow(flowAccum[i], 0.35) * 0.42));
        const straightPenalty = Math.max(0, sameDirection) * (0.12 + terrainTemplate.meanderStrength * (0.42 + lowlandMeander * 0.35));
        const turnPenalty = sameDirection < -0.62 ? 0.10 : 0;
        const sideSwing = dx * lastDy - dy * lastDx;
        const bendMag = Math.abs(sideSwing);
        const meanderWave = Math.sin((path.length + bonusSeed * 0.011) * (0.44 + terrainTemplate.meanderStrength * 0.28) + hash2(startX, startY, seed + bonusSeed) * Math.PI * 2);
        const targetBend = Math.sign(meanderWave);
        const meanderBias = targetBend !== 0 ? Math.max(0, sideSwing * targetBend) * (0.055 + terrainTemplate.meanderStrength * 0.110 + lowlandMeander * 0.070) : 0;
        const antiStraight = bendMag * (0.024 + terrainTemplate.meanderStrength * 0.070 + lowlandMeander * 0.035);
        const flowBonus = ni === preferred ? (0.36 + flowAccum[ni] * 0.18) : 0;
        const junctionPenalty = confluenceAnglePenalty(nx, ny, dx, dy, path.length);
        const noise = (hash2(nx, ny, seed + bonusSeed + step * 11) - 0.5) * 0.04;
        const value =
          elevation[ni] * 1.45 +
          uphill * 88 -
          Math.max(0, drop) * 2.05 -
          valley * 1.05 -
          valleyField[ni] * (2.00 + lowlandMeander * 0.28) -
          flowAccum[ni] * (1.22 + lowlandMeander * 0.12) -
          moisture[ni] * 0.18 -
          coastalLowland[ni] * (0.38 + lowlandMeander * 0.22) -
          (river[ni] > 0 ? 0.34 : 0) -
          flowBonus +
          slope[ni] * 0.04 +
          straightPenalty +
          turnPenalty +
          junctionPenalty * 1.10 -
          meanderBias -
          antiStraight -
          noise -
          (sea[ni] ? 0.6 : 0);

        if (value < bestValue) {
          bestValue = value;
          best = [nx, ny, dx, dy];
        }
      }
      if (!best) break;
      x = best[0];
      y = best[1];
      lastDx = best[2];
      lastDy = best[3];
    }

    const forced = forceRiverToWater(path);
    if (forced.length > path.length) {
      for (const [rx, ry] of forced.slice(path.length)) {
        const ri = indexOf(rx, ry);
        river[ri] += 0.50 + flowAccum[ri] * 0.68;
        accum += river[ri] + flowAccum[ri];
      }
    }
    return { path: forced, accum };
  }

  function traceSmallStreamPath(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    let lastDx = 0;
    let lastDy = 0;
    const path = [];
    const seen = new Set();
    for (let step = 0; step < 210; step++) {
      const i = indexOf(x, y);
      if (seen.has(i)) break;
      seen.add(i);
      path.push([x, y]);
      river[i] += 0.026 + flowAccum[i] * 0.045;
      if ((river[i] > 0.62 && path.length > 9) || sea[i]) break;
      let best = null;
      let bestValue = INF;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        const dx = nx - x;
        const dy = ny - y;
        const drop = elevation[i] - elevation[ni];
        const sameDirection = lastDx || lastDy ? (dx * lastDx + dy * lastDy) / Math.max(0.001, Math.hypot(dx, dy) * Math.hypot(lastDx, lastDy)) : 0;
        const swing = dx * lastDy - dy * lastDx;
        const lowlandMeander = clamp((1 - slope[i] * 2.0) * (0.28 + valleyField[i] * 0.80 + Math.pow(flowAccum[i], 0.35) * 0.38));
        const meanderWave = Math.sin((step + bonusSeed * 0.009) * (0.52 + terrainTemplate.meanderStrength * 0.25) + hash2(startX, startY, seed + 23123) * Math.PI * 2);
        const targetSwing = Math.sign(meanderWave);
        const lateralBonus = targetSwing !== 0 ? Math.max(0, swing * targetSwing) * (0.038 + terrainTemplate.meanderStrength * 0.070 + lowlandMeander * 0.060) : 0;
        const value = elevation[ni] * 1.10 + Math.max(0, -drop) * 17.5 - Math.max(0, drop) * 1.52 - valleyField[ni] * (1.42 + lowlandMeander * 0.18) - flowAccum[ni] * (0.74 + lowlandMeander * 0.06) - moisture[ni] * 0.16 + Math.max(0, sameDirection) * (0.062 + lowlandMeander * 0.03) - Math.abs(swing) * 0.016 - lateralBonus + (hash2(nx, ny, seed + bonusSeed + step * 13) - 0.5) * 0.090;
        if (value < bestValue) { bestValue = value; best = [nx, ny, dx, dy]; }
      }
      if (!best) break;
      x = best[0];
      y = best[1];
      lastDx = best[2];
      lastDy = best[3];
    }
    return path;
  }

  const riverPaths = [];
  const majorRiverPathSet = new Set();
  const riverScores = [];
  for (const source of majorSources) {
    const { path, accum } = traceRiverPath(source.x, source.y, 2000 + source.x * 13 + source.y * 19);
    if (path.length > 14) {
      riverPaths.push(path);
      majorRiverPathSet.add(path);
      riverScores.push(path.length * 1.35 + accum * 0.32);
    }
  }
  for (const source of sources) {
    const { path, accum } = traceRiverPath(source.x, source.y, 0);
    if (path.length > 6) {
      riverPaths.push(path);
      riverScores.push(path.length + accum * 0.18);
    }
  }

  const preliminaryMainRiverCells = new Set(riverPaths.slice().sort((a, b) => b.length - a.length).slice(0, 5).flatMap((path) => path.map(([x, y]) => `${x},${y}`)));
  const tributarySources = pickEntities(sourceCandidates
    .filter((p) => !preliminaryMainRiverCells.has(`${p.x},${p.y}`))
    .map((p) => ({ ...p, score: p.score + flowAccum[indexOf(p.x, p.y)] * 0.75 + valleyField[indexOf(p.x, p.y)] * 0.24 })), {
      max: 42 + Math.floor(rand(seed, 915) * 28),
      minDistance: 3,
      threshold: 0.16,
      seed: seed + 916,
      jitter: 0.02,
    });
  for (const source of tributarySources) {
    const { path, accum } = traceRiverPath(source.x, source.y, 4000 + source.x * 7 + source.y * 11);
    if (path.length > 8) {
      riverPaths.push(path);
      riverScores.push(path.length * 0.92 + accum * 0.17);
    }
  }

  const streamPaths = [];
  const streamSources = pickEntities(sourceCandidates
    .map((p) => ({ ...p, score: valleyField[indexOf(p.x, p.y)] * 0.46 + flowAccum[indexOf(p.x, p.y)] * 0.36 + moisture[indexOf(p.x, p.y)] * 0.18 + hash2(p.x, p.y, seed + 918) * 0.05 }))
    .filter((p) => p.score > 0.095), {
      max: 82 + Math.floor(rand(seed, 919) * 52),
      minDistance: 1.5,
      threshold: 0.062,
      seed: seed + 919,
      jitter: 0.015,
    });
  for (const source of streamSources) {
    const path = traceSmallStreamPath(source.x, source.y, 7000 + source.x * 5 + source.y * 17);
    if (path.length > 4) streamPaths.push(path);
  }

  if (riverPaths.length === 0 && sourceCandidates.length > 0) {
    const fallback = sourceCandidates.slice().sort((a, b) => b.score - a.score)[0];
    let bestSea = null;
    let bestSeaDist = INF;
    for (let y = 0; y < MAP_H; y++) {
      for (let x = 0; x < MAP_W; x++) {
        if (!sea[indexOf(x, y)]) continue;
        const d = Math.hypot(x - fallback.x, y - fallback.y);
        if (d < bestSeaDist) {
          bestSeaDist = d;
          bestSea = { x, y };
        }
      }
    }
    if (bestSea) {
      const fallbackPath = aStar(fallback, bestSea, (x, y, cx, cy) => {
        const i = indexOf(x, y);
        const ci = indexOf(cx, cy);
        if (sea[i]) return 0.25;
        const uphill = Math.max(0, elevation[i] - elevation[ci]) * 24;
        const downhill = Math.max(0, elevation[ci] - elevation[i]) * 1.8;
        return Math.max(0.24, 1 + uphill + slope[i] * 0.7 + elevation[i] * 0.8 - downhill - Math.min(0.55, river[i] * 0.1));
      });
      if (fallbackPath.length > 6) {
        let accum = 0;
        for (const [x, y] of fallbackPath) {
          const i = indexOf(x, y);
          river[i] += 0.42;
          accum += river[i];
        }
        riverPaths.push(fallbackPath);
        riverScores.push(fallbackPath.length + accum * 0.18);
      }
    }
  }

  function sanitizeDownhillRiverPath(path, tolerance = 0.040) {
    if (!path || path.length < 2) return path || [];
    const out = [path[0]];
    for (let k = 1; k < path.length; k++) {
      const [px, py] = out[out.length - 1];
      const [x, y] = path[k];
      const pi = indexOf(px, py);
      const i = indexOf(x, y);
      if (!sea[i] && elevation[i] > elevation[pi] + tolerance) break;
      out.push(path[k]);
      if (sea[i]) break;
    }
    return out.length >= 2 ? out : [];
  }
  function trimMountainHeadwaters(path) {
    if (!path || path.length < 4) return path || [];
    let start = 0;
    while (start < path.length - 3) {
      const [x, y] = path[start];
      const i = indexOf(x, y);
      if (sea[i]) break;
      if (elevation[i] <= 0.79 && (valleyField[i] >= 0.13 || flowAccum[i] >= 0.030)) break;
      start++;
    }
    return path.slice(start);
  }
  for (let r = 0; r < riverPaths.length; r++) riverPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(riverPaths[r]), 0.055);
  for (let r = riverPaths.length - 1; r >= 0; r--) if (riverPaths[r].length < 2) riverPaths.splice(r, 1);
  for (let r = 0; r < streamPaths.length; r++) streamPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(streamPaths[r]), 0.040);
  for (let r = streamPaths.length - 1; r >= 0; r--) if (streamPaths[r].length < 2) streamPaths.splice(r, 1);
  river.fill(0);
  for (const path of riverPaths) {
    const major = majorRiverPathSet.has(path);
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] += major ? (0.72 + k / 118 + flowAccum[i] * 0.98) : (0.46 + k / 170 + flowAccum[i] * 0.72);
    }
  }
  for (const path of streamPaths) {
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] += 0.020 + flowAccum[i] * 0.032;
    }
  }

  const expandedRiver = new Float32Array(river);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (river[i] <= 0) continue;
      for (const [nx, ny] of neighbors8(x, y)) {
        expandedRiver[indexOf(nx, ny)] = Math.max(expandedRiver[indexOf(nx, ny)], river[i] * 0.26);
      }
    }
  }
  river.set(expandedRiver);

  // Second fluvial pass uses the actual traced river network. Main channels cut
  // visible V-shaped valleys; lower reaches accumulate alluvial deposits.
  const fluvialElevation = new Float32Array(elevation);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i] || river[i] <= 0.02) continue;
      const r = clamp(river[i] / 2.6);
      const smallPower = smoothstep((r - 0.025) / 0.16);
      const mediumPower = smoothstep((r - 0.20) / 0.34);
      const largePower = smoothstep((r - 0.45) / 0.42);
      const actionPower = clamp(mediumPower * 0.20 + largePower * 0.44);
      const headwaterProtect = clamp((0.54 - r) / 0.54) * clamp((elevation[i] - (seaLevel + 0.08)) / 0.42);
      const highlandProtect = clamp((elevation[i] - 0.54) / 0.34) * clamp((0.58 - flowAccum[i]) / 0.58);
      const localWallProtect = clamp((slope[i] - 0.22) * 1.8) * clamp((elevation[i] - 0.54) / 0.34);
      const cutLimiter = 1 - clamp(headwaterProtect * 0.96 + highlandProtect * 0.76 + localWallProtect * 0.55);
      const channelCutRaw = terrainTemplate.fluvialAggression * actionPower * Math.pow(r, 0.82) * (0.0045 + terrainTemplate.erosion * 0.0055 + slope[i] * (0.006 + terrainTemplate.erosion * 0.008) + ridgeField[i] * (0.0015 + terrainTemplate.erosion * 0.0025));
      const valleyWidenRaw = terrainTemplate.fluvialAggression * (mediumPower * 0.08 + largePower * 0.18) * Math.pow(r, 0.92) * (0.0015 + terrainTemplate.erosion * 0.0030 + Math.max(0, elevation[i] - seaLevel) * (0.0018 + terrainTemplate.erosion * 0.0030) + valleyField[i] * (0.0015 + terrainTemplate.erosion * 0.0030));
      const maxRiverCut = 0.0018 + mediumPower * 0.0032 + largePower * 0.0065 + Math.pow(flowAccum[i], 0.65) * 0.0045;
      const channelCut = Math.min(maxRiverCut, clamp(channelCutRaw * cutLimiter));
      const valleyWiden = Math.min(maxRiverCut * 0.65, clamp(valleyWidenRaw * cutLimiter));
      const alluvium = clamp((mediumPower * 0.32 + largePower * 0.70) * Math.pow(r, 0.86) * (coastalLowland[i] * (0.010 + terrainTemplate.deposition * 0.030) + basinField[i] * (0.007 + terrainTemplate.deposition * 0.020) + (slope[i] < 0.10 ? 0.004 + terrainTemplate.deposition * 0.012 : 0)) * (1 - ridgeField[i] * 0.45));
      erosionField[i] = clamp(erosionField[i] + channelCut + valleyWiden);
      depositionField[i] = clamp(depositionField[i] + alluvium);
      depositionalLowland[i] = clamp(depositionalLowland[i] + alluvium * 5.5);
      fluvialElevation[i] = clamp(elevation[i] - channelCut - valleyWiden + alluvium, seaLevel + 0.005, 1);
      valleyField[i] = clamp(valleyField[i] + r * 0.28 + channelCut * 3.0);
      basinField[i] = clamp(basinField[i] + alluvium * 3.2);
    }
  }
  // No lateral terrain carving from traced river lines. Side-valley complexity is
  // already present in the DEM through proceduralReliefField/proceduralRavineField.

  // Template-driven deposition is limited to plausible low-energy places:
  // river mouths, basin floors, coastal plains, and slope breaks below ridges.
  const depositionElevation = new Float32Array(fluvialElevation);
  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let nearSea = 0;
      let localRiver = river[i];
      let highSide = 0;
      let lowSide = 1;
      for (let dy = -4; dy <= 4; dy++) {
        for (let dx = -4; dx <= 4; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          const ni = indexOf(nx, ny);
          const d = Math.hypot(dx, dy);
          if (d > 4.25) continue;
          if (sea[ni]) nearSea = Math.max(nearSea, 1 - d / 4.25);
          localRiver = Math.max(localRiver, river[ni] / (1 + d * 0.5));
          highSide = Math.max(highSide, fluvialElevation[ni]);
          lowSide = Math.min(lowSide, fluvialElevation[ni]);
        }
      }
      const reliefDrop = clamp((highSide - lowSide - 0.075) * 4.5);
      const lowlandPotential = clamp(
        basinField[i] * 0.44 +
        coastalLowland[i] * 0.52 +
        Math.pow(flowAccum[i], 0.56) * 0.32 +
        plain[i] * 0.18 +
        localRiver * 0.16 -
        ridgeField[i] * 0.48 -
        slope[i] * 0.52 -
        Math.max(0, fluvialElevation[i] - 0.55) * 1.35
      );
      const delta = clamp(nearSea * localRiver * coastalLowland[i] * (0.32 + terrainTemplate.deposition * 1.25) * (1 - ridgeField[i] * 0.55));
      const fan = clamp(reliefDrop * localRiver * valleyField[i] * (0.20 + terrainTemplate.deposition * 0.95) * (1 - coastalLowland[i] * 0.45));
      const lowland = clamp(lowlandPotential * terrainTemplate.deposition + delta * 0.72 + fan * 0.42);
      if (lowland <= 0.01) continue;
      deltaField[i] = clamp(deltaField[i] + delta);
      alluvialFanField[i] = clamp(alluvialFanField[i] + fan);
      depositionalLowland[i] = clamp(depositionalLowland[i] + lowland);
      depositionField[i] = clamp(depositionField[i] + lowland * 0.050);
      erosionField[i] = Math.max(0, erosionField[i] - lowland * 0.018);
      const floor = seaLevel + 0.008 + basinField[i] * 0.012 + coastalLowland[i] * 0.010;
      depositionElevation[i] = clamp(lerp(fluvialElevation[i], Math.max(floor, fluvialElevation[i] - 0.032), lowland * 0.55), seaLevel + 0.005, 1);
    }
  }
  fluvialElevation.set(depositionElevation);

  // Restore rugged summit relief after strong river incision. This prevents highlands
  // from becoming unnaturally flat or visually concave while keeping valleys cut.
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const high = clamp((fluvialElevation[i] - 0.62) / 0.26);
      const summit = high * clamp(ridgeField[i] * 1.4 - flowAccum[i] * 0.8);
      const rugged = (valueNoise(x * 2.1 + 19, y * 2.1 - 23, seed + 9661, 3.2) - 0.5) * 0.035;
      const uplift = summit * (0.018 + Math.max(0, rugged));
      if (uplift > 0) {
        fluvialElevation[i] = clamp(fluvialElevation[i] + uplift, seaLevel + 0.005, 0.985);
        erosionField[i] = Math.max(0, erosionField[i] - uplift * 0.6);
      }
    }
  }

  // Guard against one-cell over-incision: a river cell next to a 0.8-0.9 ridge
  // must not collapse to near sea level just because a routed channel crossed it.
  const guardedFluvialElevation = new Float32Array(fluvialElevation);
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i] || river[i] <= 0.10) continue;
      let highNeighbor = fluvialElevation[i];
      let meanNeighbor = 0;
      let nCount = 0;
      for (const [nx, ny] of neighbors8(x, y)) {
        const ni = indexOf(nx, ny);
        if (sea[ni]) continue;
        highNeighbor = Math.max(highNeighbor, fluvialElevation[ni]);
        meanNeighbor += fluvialElevation[ni];
        nCount++;
      }
      meanNeighbor /= Math.max(1, nCount);
      const maxAllowedDrop = 0.18 + clamp(river[i] / 2.6) * 0.06 + coastalLowland[i] * 0.08 + basinField[i] * 0.04;
      const floorFromWall = highNeighbor - maxAllowedDrop;
      const floorFromMean = meanNeighbor - 0.13;
      if (highNeighbor > 0.64 && highNeighbor - fluvialElevation[i] > 0.26) {
        guardedFluvialElevation[i] = Math.max(fluvialElevation[i], Math.min(highNeighbor - 0.08, Math.max(floorFromWall, floorFromMean)));
        erosionField[i] = Math.max(0, erosionField[i] * 0.55);
      }
    }
  }
  fluvialElevation.set(guardedFluvialElevation);

  elevation.set(fluvialElevation);

  // Broad alluvial/coastal/basin plains. The plain score alone is not enough;
  // the elevation surface must also be locally calm, otherwise every lowland
  // still reads as rugged terrain. Smooth only low, wet depositional cells and
  // leave ridges/headwaters untouched.
  for (let pass = 0; pass < 4 + Math.round(terrainTemplate.deposition * 3); pass++) {
    const nextElevation = new Float32Array(elevation);
    for (let y = 2; y < MAP_H - 2; y++) {
      for (let x = 2; x < MAP_W - 2; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const lowland = clamp(
          coastalLowland[i] * 0.72 +
          basinField[i] * 0.54 +
          depositionalLowland[i] * 0.52 +
          deltaField[i] * 0.34 +
          alluvialFanField[i] * 0.22 +
          valleyField[i] * 0.18 +
          Math.pow(flowAccum[i], 0.58) * 0.20 -
          gullyIncisionField[i] * 3.0 -
          ridgeField[i] * 0.62 -
          Math.max(0, elevation[i] - 0.54) * 1.65 -
          slope[i] * 0.74
        );
        if (lowland <= 0.12) continue;
        let sum = 0;
        let weight = 0;
        let localMin = 1;
        let localMax = 0;
        for (let dy = -2; dy <= 2; dy++) {
          for (let dx = -2; dx <= 2; dx++) {
            const nx = x + dx;
            const ny = y + dy;
            const ni = indexOf(nx, ny);
            if (sea[ni]) continue;
            const d = Math.hypot(dx, dy);
            if (d > 2.25) continue;
            localMin = Math.min(localMin, elevation[ni]);
            localMax = Math.max(localMax, elevation[ni]);
            const compatible = clamp(1 - Math.abs(elevation[ni] - elevation[i]) / 0.11);
            const w = compatible / (1 + d);
            sum += elevation[ni] * w;
            weight += w;
          }
        }
        if (weight <= 0) continue;
        const localMean = sum / weight;
        const localRelief = localMax - localMin;
        const flatBias = clamp(1 - localRelief / 0.10);
        const terrace = Math.round(localMean * 42) / 42;
        const target = lerp(localMean, terrace, 0.18 + flatBias * 0.24);
        const flattenStrength = lowland * (0.32 + terrainTemplate.deposition * 0.24 + flatBias * 0.22);
        nextElevation[i] = clamp(lerp(elevation[i], target, flattenStrength), seaLevel + 0.006, 1);
        if (lowland > 0.55) {
          depositionField[i] = clamp(depositionField[i] + lowland * (0.010 + terrainTemplate.deposition * 0.018));
          erosionField[i] = Math.max(0, erosionField[i] - lowland * 0.012);
        }
      }
    }
    elevation.set(nextElevation);
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[i] = clamp(Math.sqrt(gx * gx + gy * gy) * 11.2);
    }
  }

  // Re-trim visible river paths after fluvial reshaping changes local elevation.
  for (let r = 0; r < riverPaths.length; r++) riverPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(riverPaths[r]), 0.052);
  for (let r = riverPaths.length - 1; r >= 0; r--) if (riverPaths[r].length < 2) riverPaths.splice(r, 1);
  for (let r = 0; r < streamPaths.length; r++) streamPaths[r] = sanitizeDownhillRiverPath(trimMountainHeadwaters(streamPaths[r]), 0.022);
  for (let r = streamPaths.length - 1; r >= 0; r--) if (streamPaths[r].length < 2) streamPaths.splice(r, 1);

  function pathKey(path) {
    return path.map(([x, y]) => `${x},${y}`).join("|");
  }

  function buildPathCellSet(paths) {
    const set = new Set();
    for (const path of paths) for (const [x, y] of path) set.add(`${x},${y}`);
    return set;
  }

  function riverPathStats(path) {
    let maxRiver = 0;
    let sumRiver = 0;
    let maxFlow = 0;
    let sumFlow = 0;
    let populatedCorridor = 0;
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      maxRiver = Math.max(maxRiver, river[i]);
      sumRiver += river[i];
      maxFlow = Math.max(maxFlow, flowAccum[i]);
      sumFlow += flowAccum[i];
      populatedCorridor += plain[i] * 0.18 + valleyField[i] * 0.28 + coastalLowland[i] * 0.10 + basinField[i] * 0.08;
    }
    const [lx, ly] = path[path.length - 1];
    const li = indexOf(lx, ly);
    const outletToWater = Boolean(sea[li] || lake[li]);
    const lowerReach = path.slice(Math.max(0, path.length - Math.min(path.length, 8)));
    const lowerReachStrength = lowerReach.reduce((sum, [x, y]) => sum + river[indexOf(x, y)], 0) / Math.max(1, lowerReach.length);
    const meanRiver = sumRiver / Math.max(1, path.length);
    const meanFlow = sumFlow / Math.max(1, path.length);
    const corridorMean = populatedCorridor / Math.max(1, path.length);
    const score =
      path.length * 0.92 +
      maxRiver * 8.4 +
      meanRiver * 4.4 +
      maxFlow * 8.2 +
      meanFlow * 2.8 +
      lowerReachStrength * 3.2 +
      corridorMean * 5.2 +
      (outletToWater ? 5.0 : 0);
    return { length: path.length, maxRiver, meanRiver, maxFlow, meanFlow, lowerReachStrength, corridorMean, outletToWater, score };
  }

  let rankedRivers = riverPaths
    .map((path, i) => ({ path, score: riverScores[i] || 0, stats: riverPathStats(path), key: pathKey(path) }))
    .filter((item) => item.path.length >= 5)
    .sort((a, b) => (b.stats.score + b.score * 0.25) - (a.stats.score + a.score * 0.25));

  let mainRivers = rankedRivers
    .filter((item) => item.stats.length >= 8)
    .slice(0, Math.min(8, rankedRivers.length))
    .map((item) => item.path);

  if (mainRivers.length === 0 && riverPaths.length > 0) mainRivers.push(riverPaths[0]);
  if (mainRivers.length === 0) {
    let start = null;
    let startScore = -INF;
    for (let y = 4; y < MAP_H - 4; y++) {
      for (let x = 4; x < MAP_W - 4; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        const score = elevation[i] * 0.55 + moisture[i] * 0.35 - slope[i] * 0.15;
        if (score > startScore) {
          startScore = score;
          start = { x, y };
        }
      }
    }
    if (start) {
      let goal = null;
      let goalDist = INF;
      for (let y = 0; y < MAP_H; y++) {
        for (let x = 0; x < MAP_W; x++) {
          if (!sea[indexOf(x, y)]) continue;
          const d = Math.hypot(x - start.x, y - start.y);
          if (d < goalDist) {
            goalDist = d;
            goal = { x, y };
          }
        }
      }
      if (goal) {
        const fallbackPath = aStar(start, goal, (x, y, cx, cy) => {
          const i = indexOf(x, y);
          const ci = indexOf(cx, cy);
          if (sea[i]) return 0.2;
          const uphillBias = Math.max(0, elevation[i] - elevation[ci]) * 22;
          const downhillBias = Math.max(0, elevation[ci] - elevation[i]) * 1.7;
          return Math.max(0.25, 1 + uphillBias + slope[i] * 0.65 + elevation[i] * 0.8 - downhillBias);
        });
        if (fallbackPath.length > 4) {
          riverPaths.push(fallbackPath);
          mainRivers.push(fallbackPath);
          for (const [x, y] of fallbackPath) river[indexOf(x, y)] += 0.55;
        }
      }
    }
  }

  // v4: 急峻な地形では自然流下トレースが短く切れる seed があるため、
  // 高地から海へ抜ける中〜大規模河川の骨格を数本だけ補完する。
  if (mainRivers.length < 4 && sourceCandidates.length > 0) {
    const usedKeys = new Set(mainRivers.map((path) => pathKey(path)));
    const starts = sourceCandidates.slice()
      .sort((a, b) => (b.score + elevation[indexOf(b.x, b.y)] * 0.8 + valleyField[indexOf(b.x, b.y)] * 0.6) - (a.score + elevation[indexOf(a.x, a.y)] * 0.8 + valleyField[indexOf(a.x, a.y)] * 0.6));
    for (const start of starts) {
      if (mainRivers.length >= 4) break;
      const tooClose = mainRivers.some((path) => path.some(([px, py], k) => k % 8 === 0 && Math.hypot(px - start.x, py - start.y) < 10));
      if (tooClose) continue;
      const goal = nearestWaterGoal(start);
      if (!goal) continue;
      const path = aStar(start, goal, (x, y, cx, cy) => {
        const i = indexOf(x, y);
        const ci = indexOf(cx, cy);
        if (sea[i]) return 0.18;
        const uphill = Math.max(0, elevation[i] - elevation[ci]);
        const downhill = Math.max(0, elevation[ci] - elevation[i]);
        return Math.max(0.22, 1 + uphill * 42 + slope[i] * 0.42 + elevation[i] * 0.32 - downhill * 2.4 - valleyField[i] * 1.65 - flowAccum[i] * 1.20 - moisture[i] * 0.18 - coastalLowland[i] * 0.38);
      });
      if (path.length < 9) continue;
      const key = pathKey(path);
      if (usedKeys.has(key)) continue;
      usedKeys.add(key);
      mainRivers.push(path);
      riverPaths.push(path);
      riverScores.push(path.length * 1.05);
      for (let k = 0; k < path.length; k++) {
        const [rx, ry] = path[k];
        river[indexOf(rx, ry)] = Math.max(river[indexOf(rx, ry)], 0.62 + k / 180 + flowAccum[indexOf(rx, ry)] * 0.72);
      }
    }
  }

  // v5: 最終的に主河川が内陸で途切れる場合は、海または湖まで河口部を補完する。
  // これは地形生成後の河川作用を明示的に効かせ、山地から海への侵食軸を保証するため。
  for (let r = 0; r < mainRivers.length; r++) {
    const path = mainRivers[r];
    if (!path || path.length < 2) continue;
    const connected = path.some(([x, y], k) => k > path.length * 0.45 && neighbors8(x, y).some(([nx, ny]) => sea[indexOf(nx, ny)] || lake[indexOf(nx, ny)]));
    if (connected) continue;
    const forced = forceRiverToWater(path);
    if (forced.length > path.length) {
      mainRivers[r] = forced;
      riverPaths.push(forced);
      riverScores.push(forced.length * 1.12);
      for (let k = path.length; k < forced.length; k++) {
        const [rx, ry] = forced[k];
        const ri = indexOf(rx, ry);
        river[ri] = Math.max(river[ri], 0.70 + k / 180 + flowAccum[ri] * 0.70);
        valleyField[ri] = clamp(valleyField[ri] + 0.18);
      }
    }
  }

  const mainRiverCells = buildPathCellSet(mainRivers);
  const mainRiverKeys = new Set(mainRivers.map((path) => pathKey(path)));

  function pathTouchesMain(path) {
    for (const [x, y] of path) {
      if (mainRiverCells.has(`${x},${y}`)) return true;
      for (const [nx, ny] of neighbors8(x, y)) {
        if (mainRiverCells.has(`${nx},${ny}`)) return true;
      }
    }
    return false;
  }

  rankedRivers = riverPaths
    .map((path, i) => ({ path, score: riverScores[i] || 0, stats: riverPathStats(path), key: pathKey(path) }))
    .filter((item) => item.path.length >= 5)
    .sort((a, b) => (b.stats.score + b.score * 0.25) - (a.stats.score + a.score * 0.25));

  const tributaryRivers = [];
  const hiddenRiverPaths = [];
  for (const item of rankedRivers) {
    if (mainRiverKeys.has(item.key)) continue;
    const joinsMain = pathTouchesMain(item.path);
    const visibleMedium =
      item.stats.score >= 18 &&
      item.stats.length >= 7 &&
      (joinsMain || item.stats.outletToWater || item.stats.maxRiver >= 0.95 || item.stats.lowerReachStrength >= 0.70);
    if (visibleMedium) tributaryRivers.push(item.path);
    else hiddenRiverPaths.push(item.path);
  }

  for (const path of mainRivers) {
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] = Math.max(river[i], 0.92 + k / 150 + flowAccum[i] * 0.96);
    }
  }
  for (const path of tributaryRivers) {
    for (let k = 0; k < path.length; k++) {
      const [x, y] = path[k];
      const i = indexOf(x, y);
      river[i] = Math.max(river[i], 0.58 + k / 195 + flowAccum[i] * 0.62);
    }
  }

  function traceFlowLinkedMinorStream(startX, startY, bonusSeed = 0) {
    let x = startX;
    let y = startY;
    const path = [];
    const seen = new Set();
    for (let step = 0; step < 120; step++) {
      const i = indexOf(x, y);
      if (sea[i] || seen.has(i)) break;
      seen.add(i);
      path.push([x, y]);
      if (path.length > 7 && river[i] > 0.42) break;
      let next = flowTo[i];
      if (next < 0 || next === i || sea[next]) break;
      let best = next;
      let bestScore = elevation[next] * 1.05 - flowAccum[next] * 0.85 - valleyField[next] * 1.20 - moisture[next] * 0.10;
      const cx = x;
      const cy = y;
      // Micro-streams can braid into the closest descent when flowTo falls into a tiny sink.
      for (const [nx, ny] of neighbors8(cx, cy)) {
        const ni = indexOf(nx, ny);
        if (sea[ni]) continue;
        const uphill = Math.max(0, elevation[ni] - elevation[i]);
        if (uphill > 0.024 && flowAccum[ni] < flowAccum[i] + 0.006) continue;
        const score = elevation[ni] * 1.05 + uphill * 16 - flowAccum[ni] * 0.82 - valleyField[ni] * 1.22 - moisture[ni] * 0.10 + (hash2(nx, ny, seed + bonusSeed + step * 19) - 0.5) * 0.035;
        if (score < bestScore) {
          bestScore = score;
          best = ni;
        }
      }
      if (best < 0 || best === i) break;
      x = best % MAP_W;
      y = Math.floor(best / MAP_W);
    }
    return path;
  }

  const minorCandidateCells = [];
  for (let y = 3; y < MAP_H - 3; y += 1) {
    for (let x = 3; x < MAP_W - 3; x += 1) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      if (elevation[i] < 0.30 || elevation[i] > 0.96) continue;
      const drainage = valleyField[i] * 0.52 + Math.pow(flowAccum[i], 0.48) * 0.38 + moisture[i] * 0.18 + slope[i] * 0.08 - ridgeField[i] * 0.10;
      const stochastic = hash2(x, y, seed + 9340);
      if (drainage > 0.085 && stochastic > 0.10) {
        minorCandidateCells.push({ x, y, score: drainage + stochastic * 0.055 });
      }
    }
  }
  const minorSources = pickEntities(minorCandidateCells, {
    max: 90 + Math.floor(rand(seed, 9341) * 60),
    minDistance: 2,
    threshold: 0.070,
    seed: seed + 9342,
    jitter: 0.02,
  });
  const derivedSmallStreams = [];
  const occupiedMinorStarts = new Set();
  for (const source of minorSources) {
    const startKey = `${source.x},${source.y}`;
    if (occupiedMinorStarts.has(startKey)) continue;
    const path = traceFlowLinkedMinorStream(source.x, source.y, 11000 + source.x * 13 + source.y * 17);
    if (path.length >= 3) {
      derivedSmallStreams.push(path);
      for (const [x, y] of path.slice(0, 4)) occupiedMinorStarts.add(`${x},${y}`);
      for (let k = 0; k < path.length; k++) {
        const [x, y] = path[k];
        const i = indexOf(x, y);
        river[i] = Math.max(river[i], 0.045 + Math.min(0.16, flowAccum[i] * 0.10) + Math.min(0.055, k / 1900));
      }
    }
  }

  const smallStreams = streamPaths.filter((path) => path.length >= 4)
    .concat(hiddenRiverPaths.filter((path) => path.length >= 5))
    .concat(derivedSmallStreams);

  for (const path of smallStreams) {
    if (!path || path.length < 3) continue;
    const strength = path.length > 24 ? 0.020 : 0.013;
    addPathIncision(visibleRavineField, path, strength, path.length > 18 ? 2 : 1);
  }

  // v9: keep enclosed盆地 as habitable basins rather than over-incised pits.
  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const basinFloor = clamp(basinField[i] * (1 - slope[i] * 3.0) * (1 - coastalLowland[i] * 0.7) * (1 - Math.min(1, river[i] * 0.7)));
      if (basinFloor <= 0.16) continue;
      elevation[i] = clamp(elevation[i] + basinFloor * 0.020, seaLevel + 0.008, 1);
      depositionalLowland[i] = clamp(depositionalLowland[i] + basinFloor * 0.12);
      plain[i] = clamp(plain[i] + basinFloor * 0.10);
    }
  }
  // Final one-cell canyon guard.  Apply it to all land cells, not only river
  // cells, because a traced channel or earlier basin operation can leave a 0.3
  // cell directly beside a 0.9 ridge.  The guard preserves broad basins/coasts but
  // prevents single-cell cliff trenches.
  for (let pass = 0; pass < 4; pass++) {
    const guarded = new Float32Array(elevation);
    for (let y = 1; y < MAP_H - 1; y++) {
      for (let x = 1; x < MAP_W - 1; x++) {
        const i = indexOf(x, y);
        if (sea[i]) continue;
        let highNeighbor = elevation[i];
        let meanNeighbor = 0;
        let nCount = 0;
        let seaTouch = 0;
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          if (sea[ni]) { seaTouch++; continue; }
          highNeighbor = Math.max(highNeighbor, elevation[ni]);
          meanNeighbor += elevation[ni];
          nCount++;
        }
        meanNeighbor /= Math.max(1, nCount);
        if (highNeighbor > 0.64 && highNeighbor - elevation[i] > 0.22) {
          const coastalAllowance = coastalLowland[i] * 0.10 + (seaTouch ? 0.08 : 0);
          const basinAllowance = basinField[i] * 0.055;
          const riverAllowance = clamp(river[i] / 2.8) * 0.035;
          const allowedDrop = 0.19 + coastalAllowance + basinAllowance + riverAllowance;
          guarded[i] = Math.max(elevation[i], Math.max(highNeighbor - allowedDrop, meanNeighbor - 0.085));
          valleyField[i] = clamp(valleyField[i] * 0.86);
          erosionField[i] = Math.max(0, erosionField[i] * 0.40);
        }
      }
    }
    elevation.set(guarded);
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)];
      const gy = elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)];
      slope[i] = clamp(Math.sqrt(gx * gx + gy * gy) * 10.5);
    }
  }

  prefectureMask = makePrefectureMask(seed, sea, elevation, slope, river);
  prefectureBorder = extractMaskBorder(prefectureMask, sea);
  const regionalPrefectures = generateRegionalPrefectures(seed, sea, elevation, slope, river, ridgeField, flowAccum, prefectureMask);
  const prefectureRegionId = regionalPrefectures.regionId;
  const regionalDebug = regionalPrefectures.debug;
  const regionalPrefectureBorders = extractRegionBorderSegments(prefectureRegionId, sea);

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const low = 1 - clamp((elevation[i] - 0.28) / 0.4);
      const flat = 1 - slope[i];
      const valleyPlain = valleyField[i] * 0.44 + basinField[i] * 0.36 + coastalLowland[i] * 0.55 + depositionalLowland[i] * 0.34 + deltaField[i] * 0.28 + alluvialFanField[i] * 0.20;
      plain[i] = clamp(low * 0.44 + flat * 0.58 + valleyPlain - ridgeField[i] * 0.28 - (elevation[i] > 0.62 ? 0.48 : 0));

      let nearRiver = 0;
      for (let dy = -4; dy <= 4; dy++) {
        for (let dx = -4; dx <= 4; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          nearRiver = Math.max(nearRiver, river[indexOf(nx, ny)] / (1 + Math.hypot(dx, dy)));
        }
      }

      const fan = clamp(Math.max(alluvialFanField[i], valleyField[i] * (1 - coastalLowland[i]) * (elevation[i] > 0.34 && elevation[i] < 0.58 ? 0.9 : 0.35)) * (1 - slope[i] * 0.55));
      floodplain[i] = clamp(nearRiver * plain[i] * 0.92 + coastalLowland[i] * nearRiver * 0.22 + deltaField[i] * 0.18);
      agriculture[i] = clamp(plain[i] * 0.58 + fan * 0.30 + basinField[i] * 0.2 + depositionalLowland[i] * 0.24 + deltaField[i] * 0.18 + moisture[i] * 0.14 + clamp(nearRiver) * 0.32 - slope[i] * 0.34 - ridgeField[i] * 0.18 - floodplain[i] * 0.06);
    }
  }

  for (let y = 2; y < MAP_H - 2; y++) {
    for (let x = 2; x < MAP_W - 2; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      let seaNear = 0;
      let riverNear = 0;
      let sheltered = 0;

      for (let dy = -5; dy <= 5; dy++) {
        for (let dx = -5; dx <= 5; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (!inside(nx, ny)) continue;
          const d = Math.hypot(dx, dy);
          if (sea[indexOf(nx, ny)]) seaNear += 1 / (1 + d);
          riverNear = Math.max(riverNear, river[indexOf(nx, ny)] / (1 + d));
        }
      }

      for (let dy = -2; dy <= 2; dy++) {
        for (let dx = -2; dx <= 2; dx++) {
          const nx = x + dx;
          const ny = y + dy;
          if (inside(nx, ny) && !sea[indexOf(nx, ny)]) sheltered += 1;
        }
      }

      const isDelta = (riverNear > 0.22 && coastalLowland[i] > 0.18) || deltaField[i] > 0.16;
      const bayShelter = sheltered * 0.012 + seaNear * 0.055 + coastalLowland[i] * 0.16;
      portSuitability[i] = clamp(bayShelter + riverNear * 0.24 + (isDelta ? 0.22 : 0) + deltaField[i] * 0.18 + depositionalLowland[i] * 0.08 + plain[i] * 0.08 - slope[i] * 0.48 - ridgeField[i] * 0.16);
    }
  }

  for (let y = 3; y < MAP_H - 3; y++) {
    for (let x = 3; x < MAP_W - 3; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const r = river[i];
      if (r < 0.2 || r > 1.85) continue;
      let bankPlain = 0;
      for (const [nx, ny] of neighbors8(x, y)) bankPlain += plain[indexOf(nx, ny)];
      crossingSuitability[i] = clamp(r * 0.34 + (bankPlain / 8) * 0.54 + valleyField[i] * 0.18 - slope[i] * 0.55 - floodplain[i] * 0.06);
    }
  }

  for (let y = 4; y < MAP_H - 4; y++) {
    for (let x = 4; x < MAP_W - 4; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const e = elevation[i];
      if (e < 0.43 || e > 0.82) continue;
      const ewHigh = (elevation[indexOf(x - 3, y)] + elevation[indexOf(x + 3, y)]) / 2;
      const nsHigh = (elevation[indexOf(x, y - 3)] + elevation[indexOf(x, y + 3)]) / 2;
      const diagLow = Math.min(
        elevation[indexOf(x - 3, y - 3)],
        elevation[indexOf(x + 3, y + 3)],
        elevation[indexOf(x - 3, y + 3)],
        elevation[indexOf(x + 3, y - 3)]
      );
      passSuitability[i] = clamp((Math.max(ewHigh, nsHigh) - e) * 2.2 + (e - diagLow) * 0.55 + valleyField[i] * 0.28 - ridgeField[i] * 0.18 - slope[i] * 0.2);
    }
  }

  for (let y = 1; y < MAP_H - 1; y++) {
    for (let x = 1; x < MAP_W - 1; x++) {
      const i = indexOf(x, y);
      if (sea[i]) continue;
      const gx = Math.abs(elevation[indexOf(x + 1, y)] - elevation[indexOf(x - 1, y)]);
      const gy = Math.abs(elevation[indexOf(x, y + 1)] - elevation[indexOf(x, y - 1)]);
      const slopeBreak = clamp((gx + gy) * 3.2 + Math.max(0, slope[i] - 0.28) * 0.72);
      const majorRiver = clamp(Math.max(0, river[i] - 0.34) * 1.45 + Math.max(0, flowAccum[i] - 0.42) * 0.58);
      const basinRim = clamp(basinField[i] * Math.max(0, slope[i] - 0.16) * 1.25 + ridgeField[i] * basinField[i] * 0.32);
      naturalBarrierScore[i] = clamp(
        arcSpineField[i] * 0.80 +
        branchRidgeField[i] * 0.62 +
        ridgeField[i] * 0.54 +
        majorRiver * 0.62 +
        slopeBreak * 0.34 +
        basinRim * 0.36 -
        valleyField[i] * 0.30 -
        depositionalLowland[i] * 0.42 -
        coastalLowland[i] * 0.20 -
        plain[i] * 0.18
      );
    }
  }

  function countWaterComponents(mask, minArea = 1) {
    const seen = new Uint8Array(SIZE);
    let count = 0;
    for (let i = 0; i < SIZE; i++) {
      if (!mask[i] || seen[i]) continue;
      const queue = [i];
      seen[i] = 1;
      let area = 0;
      for (let q = 0; q < queue.length; q++) {
        const cur = queue[q];
        area++;
        const x = cur % MAP_W;
        const y = Math.floor(cur / MAP_W);
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          if (!mask[ni] || seen[ni]) continue;
          seen[ni] = 1;
          queue.push(ni);
        }
      }
      if (area >= minArea) count++;
    }
    return count;
  }

  function countSmallLandIslands(maxArea = 8) {
    const seen = new Uint8Array(SIZE);
    let count = 0;
    for (let i = 0; i < SIZE; i++) {
      if (sea[i] || seen[i]) continue;
      const queue = [i];
      seen[i] = 1;
      let area = 0;
      let touchesEdge = false;
      for (let q = 0; q < queue.length; q++) {
        const cur = queue[q];
        area++;
        const x = cur % MAP_W;
        const y = Math.floor(cur / MAP_W);
        if (x === 0 || y === 0 || x === MAP_W - 1 || y === MAP_H - 1) touchesEdge = true;
        for (const [nx, ny] of neighbors8(x, y)) {
          const ni = indexOf(nx, ny);
          if (sea[ni] || seen[ni]) continue;
          seen[ni] = 1;
          queue.push(ni);
        }
      }
      if (!touchesEdge && area <= maxArea) count++;
    }
    return count;
  }

  const spineValues = [...arcSpineField].filter((_, i) => !sea[i]).sort((a, b) => b - a);
  const strongSpineSample = Math.max(1, Math.floor(spineValues.length * 0.05));
  const primarySpineStrength = spineValues.slice(0, strongSpineSample).reduce((sum, value) => sum + value, 0) / strongSpineSample;
  const riverConnectivityRate = mainRivers.length
    ? mainRivers.filter((path) => path.some(([x, y], k) => k > path.length * 0.45 && neighbors8(x, y).some(([nx, ny]) => sea[indexOf(nx, ny)] || lake[indexOf(nx, ny)]))).length / mainRivers.length
    : 0;
  const depositionLowlandArea = [...depositionalLowland].filter((value, i) => !sea[i] && value > 0.24).length;
  const terrainDebug = {
    primarySpineStrength,
    riverConnectivityRate,
    smallIslandCount: countSmallLandIslands(8),
    largeInlandLakeCount: countWaterComponents(Float32Array.from(lake, (value) => value ? 1 : 0), 120),
    depositionLowlandArea,
    smallStreamCount: smallStreams.length,
    erosionGullyCount: erosionGullyPaths.length,
    branchRavineCount: 0,
  };

  return {
    terrainTemplate,
    seaLevel,
    elevation,
    moisture,
    slope,
    sea,
    ocean,
    lake,
    river,
    floodplain,
    plain,
    agriculture,
    ridgeField,
    valleyField,
    visibleRavineField,
    surfaceTextureField,
    basinField,
    coastalLowland,
    flowAccum,
    erosionField,
    depositionField,
    arcSpineField,
    branchRidgeField,
    depositionalLowland,
    alluvialFanField,
    deltaField,
    naturalBarrierScore,
    portSuitability,
    crossingSuitability,
    passSuitability,
    prefectureMask,
    prefectureBorder,
    prefectureRegionId,
    regionalDebug,
    terrainDebug,
    regionalPrefectureBorders,
    riverPaths,
    mainRivers,
    tributaryRivers,
    smallStreams,
  };
}