diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md
index e67b984ffd..f41ed8fec4 100644
--- a/docs/src/inner_layer.md
+++ b/docs/src/inner_layer.md
@@ -12,6 +12,7 @@ Two layer models are available, selected by `[Tearing] inner_model`:
 
 For each surface the `Dispersion` search scans the complex frequency plane, extracts the roots of the dispersion relation and classifies them, and `[Tearing.search]` sets the scan.
 With `[Tearing.CriticalResonantField]` the stage also computes, for the `SLAYER` layer, the resonant field at which a static error field overcomes the layer's shielding: the error-field penetration threshold.
+`[Tearing.IslandFlattening]` flattens the equilibrium pressure across islands of given width and reports the stored energy and pressure they remove; it needs no layer solve and runs alone when it is the only part of `[Tearing]`.
 
 The same GGJ layer serves the `[Match]` stage of the [Workflow](workflow.md), which does the complementary job: instead of finding the free modes it imposes the layer response on the driven force-free states and writes the matched solution to `SingularSurfaces/Match/`.
 
@@ -48,6 +49,24 @@ Each surface uses its own stability parameter ``\alpha = S^{-1/3}(\Delta_\mathrm
 Penetration is the first resolved saddle-node bifurcation on the path from ``Q_0`` toward the electron fluid at ``Q_e``, above which the shielded branch no longer exists; small numerical extrema and the unshielded zone next to ``Q_e`` are excluded, and a surface without a resolved bifurcation reports NaN.
 `br_crit` is normalized by the toroidal field, and `br_crit_tesla = br_crit·|B0|` is the quantity to compare with `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, the total resonant field of the ideal response; the comparison is approximate, since the theory's ``b_r`` is a cylindrical field at the layer.
 
+### Island flattening
+
+Inside a magnetic island the pressure is flattened along the field, so ``dp/d\psi = 0`` there; outside it the gradient is unchanged.
+For a set of islands of full width ``w`` (ψ_N) centred on their rational surfaces, the flattened intervals are merged where islands overlap, so an overlap is counted once, and the pressure is anchored at its value on the last profile point and integrated inward, as in Fortran GPEC's saturated-island profile.
+For merged intervals ``[a_k, b_k]`` the integral is exact:
+
+```math
+p_\mathrm{mod}(\psi) = p(\psi) + \sum_{k:\, b_k > \psi} \left[p(b_k) - p(\max(\psi, a_k))\right],
+\qquad
+\Delta W = \tfrac{3}{2}\int (p - p_\mathrm{mod})\, dV .
+```
+
+Every point inward of an interval drops by ``p(a_k) - p(b_k)``, the pedestal-top pressure change when the interval is at the edge.
+Given a bootstrap profile (the `jbs_B` keyword of `island_flattening`), the result also carries the part the islands remove: to leading order ``j_\mathrm{bs} \propto p'``, so it vanishes on the flattened intervals and is unchanged elsewhere, and `jbs_fraction_lost` is the fraction of its volume integral on them.
+For ``\langle j_\mathrm{bs}\cdot B\rangle`` that is a fraction of ``\int\langle j_\mathrm{bs}\cdot B\rangle\,dV``; for the toroidal bootstrap current pass ``j_\phi\langle R^{-1}\rangle``, since ``dA = \langle R^{-1}\rangle\,dV/2\pi`` on axisymmetric surfaces.
+The pressure is anchored at the last point of the equilibrium profile grid, so for an island reaching the edge the result depends on where that grid ends.
+The flattening is static: the equilibrium is not re-solved with the modified pressure, and island widths are inputs rather than evolved.
+
 ## Running it
 
 From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` and `vac_flag = true` (or `"Galerkin"` with `rpec_flag = true`), which provides the Δ′ matrix, and with a kinetic file on the equilibrium:
@@ -78,6 +97,15 @@ D_c_type = "none"                 # critical-Δ offset from parallel transport (
 P_perp_model = "D_perp"           # how the perpendicular Prandtl number is built (below)
 ```
 
+The flattening takes full island widths in ψ_N keyed by surface; every q = m/n surface of the equilibrium gets its key's width:
+
+```toml
+[Tearing.IslandFlattening]
+widths = { "2/1" = 0.05, "3/2" = 0.03 }
+```
+
+In a script, `island_flattening(ffs.equil; widths=Dict("2/1" => 0.05))` gives the same result.
+
 In a script, `solve(TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl()), SLAYERModel(); config=TearingConfig(...))` runs the same stage, and `solve_inner(GGJModel(), p, γ)` evaluates one GGJ layer at a growth rate ``\gamma``.
 GGJ's η and ρ default from the kinetic data (η by the resistivity closure SLAYER uses, ρ = μ·m_p·n_e); pass `eta` and `rho` to `TearingProblem` to override them.
 
@@ -103,6 +131,9 @@ The DIII-D-like kinetic file of the examples carries analytic χ profiles of the
 | `Tearing/Diagnostics/ValidRoots`, `Poles`, `FilteredRoots` | the root classification of the scan |
 | `Tearing/Scan/Surface_<k>/` | the scanned Δ(Q) with `store_scan = true` |
 | `Tearing/CriticalResonantField/br_crit`, `br_crit_tesla`, `Q_peak`, `Q0` | the penetration threshold as ``b_r/B_\phi`` and in tesla, the Q of the balance maximum and the natural rotation |
+| `Tearing/IslandFlattening/W`, `Delta_W`, `Delta_W_over_W`, `Delta_p0` | the stored energy, the part the islands remove and its fraction, and the on-axis pressure drop |
+| `Tearing/IslandFlattening/Islands/`, `Intervals/` | each island (`rational_m`, `rational_n`, `rational_psi`, its clipped interval) and the merged flattened intervals with the pressure drop `Delta_p` across each |
+| `Tearing/IslandFlattening/Profiles/psi`, `p`, `p_mod` | the equilibrium and flattened pressure (Pa), with `jbs_B` and `jbs_B_mod` when a bootstrap profile was given |
 
 ## Numerical notes
 
@@ -175,3 +206,9 @@ Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion]
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField]
 ```
+
+### IslandFlattening
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.IslandFlattening]
+```
diff --git a/docs/src/workflow.md b/docs/src/workflow.md
index 38b4cdd8ba..fe9c911324 100644
--- a/docs/src/workflow.md
+++ b/docs/src/workflow.md
@@ -161,7 +161,7 @@ The keys each section accepts are the fields of its control struct, listed in th
 | `[ForcingTerms]`, `[[ForcingTerms.coil_set]]` | the external field: a file or coil sets | [Forcing Terms](forcing_terms.md) |
 | `[PerturbedEquilibrium]` | the plasma response stage and its singular-coupling options | [Perturbed Equilibrium](perturbed_equilibrium.md) |
 | `[KineticForces]` | the NTV torque stage and the kinetic solve of `ForceFreeStates` | [Kinetic Forces](kinetic_forces.md) |
-| `[Tearing]`, `[Tearing.search]`, `[Tearing.CriticalResonantField]`, `[SLAYER]` | the tearing stage, its root search, the critical field and the SLAYER layer model | [Tearing](inner_layer.md) |
+| `[Tearing]`, `[Tearing.search]`, `[Tearing.CriticalResonantField]`, `[Tearing.IslandFlattening]`, `[SLAYER]` | the tearing stage, its root search, the critical field, the island pressure flattening and the SLAYER layer model | [Tearing](inner_layer.md) |
 | `[ErrorFields]`, `[ErrorFields.MonteCarlo]`, `[ErrorFields.scenario]`, `[ErrorFields.Risk]`, `[ErrorFields.NTV]` | the error-field tolerance stage | [Error Fields](error_fields.md) |
 
 ## Running it
@@ -216,7 +216,7 @@ All results are written to one HDF5 file (`[Run] output_file`, default `gpec.h5`
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/` |
 | `KineticForces/` | `<method>/` torque and energy profiles; `PerSpecies/<species>/<method>/` for multi-ion runs |
 | `ErrorFields/` | `CoilSensitivities/`, `MonteCarlo/`, `Risk/`, `NTV/` |
-| `Tearing/` | `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/` |
+| `Tearing/` | `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`, `IslandFlattening/` |
 | `SurfaceGeometries/` | Plasma and wall surface point clouds for visualization |
 
 ### Replaying a run
diff --git a/regression-harness/cases/diiid_n1_island_flattening.toml b/regression-harness/cases/diiid_n1_island_flattening.toml
new file mode 100644
index 0000000000..31ac6fdc61
--- /dev/null
+++ b/regression-harness/cases/diiid_n1_island_flattening.toml
@@ -0,0 +1,68 @@
+# Regression case: DIII-D-like n=1 with static pressure flattening across prescribed islands.
+# Same deck as diiid_n1 plus [Tearing.IslandFlattening] on the 3/2, 2/1 and 3/1 surfaces, so the
+# merged intervals, the stored energy they remove and the flattened pressure are tracked.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_n1_island_flattening"
+description = "DIII-D-like equilibrium, n=1, pressure flattened across 3/2, 2/1 and 3/1 islands"
+example_dir = "examples/DIIID-like_ideal_example"
+
+[overrides]
+"Tearing.IslandFlattening.widths" = { "3/2" = 0.04, "2/1" = 0.06, "3/1" = 0.03 }
+
+[quantities.W]
+h5path = "Tearing/IslandFlattening/W"
+type = "real_scalar"
+extract = "value"
+label = "stored energy W"
+noise_threshold = 1e-6
+order = 10
+
+[quantities.delta_W]
+h5path = "Tearing/IslandFlattening/Delta_W"
+type = "real_scalar"
+extract = "value"
+label = "stored energy removed"
+noise_threshold = 1e-6
+order = 11
+
+[quantities.delta_W_over_W]
+h5path = "Tearing/IslandFlattening/Delta_W_over_W"
+type = "real_scalar"
+extract = "value"
+label = "fraction of W removed"
+noise_threshold = 1e-12
+order = 12
+
+[quantities.delta_p0]
+h5path = "Tearing/IslandFlattening/Delta_p0"
+type = "real_scalar"
+extract = "value"
+label = "on-axis pressure drop"
+noise_threshold = 1e-6
+order = 13
+
+[quantities.island_psi]
+h5path = "Tearing/IslandFlattening/Islands/rational_psi"
+type = "real_vector"
+extract = "all_real"
+label = "island rational surfaces"
+noise_threshold = 1e-12
+order = 20
+
+[quantities.interval_delta_p]
+h5path = "Tearing/IslandFlattening/Intervals/Delta_p"
+type = "real_vector"
+extract = "all_real"
+label = "pressure drop per merged interval"
+noise_threshold = 1e-6
+order = 21
+
+[quantities.p_mod]
+h5path = "Tearing/IslandFlattening/Profiles/p_mod"
+type = "real_vector"
+extract = "checksum"
+label = "flattened pressure (checksum)"
+noise_threshold = 0
+order = 30
diff --git a/src/Deck.jl b/src/Deck.jl
index ea0f4b17ce..2f99ec7951 100644
--- a/src/Deck.jl
+++ b/src/Deck.jl
@@ -360,7 +360,9 @@ end
 """
 function tearing_from_deck(inputs::AbstractDict)
     haskey(inputs, "Tearing") || return nothing
-    table = inputs["Tearing"]
+    # A [Tearing] holding only [Tearing.IslandFlattening] runs the flattening without a layer solve.
+    table = _without(inputs["Tearing"], "IslandFlattening")
+    isempty(table) && haskey(inputs["Tearing"], "IslandFlattening") && return nothing
     get(table, "omega_E_kHz", Dict{String,Any}()) isa AbstractDict ||
         throw(ArgumentError("[Tearing] omega_E_kHz must be a table keyed by m/n, e.g. omega_E_kHz = { \"2/1\" = 3.0 }"))
     fields = _struct_fields(TearingProblem, (:ffs, :critical_resonant_field))
@@ -377,6 +379,21 @@ function tearing_from_deck(inputs::AbstractDict)
     return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config)
 end
 
+"""
+    island_flattening_from_deck(inputs) -> Union{Nothing,Dict{String,Float64}}
+
+The full island widths (ψ_N) keyed by `"m/n"` of `[Tearing.IslandFlattening] widths`; `nothing` without that section.
+"""
+function island_flattening_from_deck(inputs::AbstractDict)
+    table = get(get(inputs, "Tearing", Dict{String,Any}()), "IslandFlattening", nothing)
+    table === nothing && return nothing
+    kwargs = section_keywords(table, "Tearing.IslandFlattening", Dict{Symbol,Any}(:widths => AbstractDict))
+    widths = get(kwargs, :widths, nothing)
+    (widths isa AbstractDict && !isempty(widths) && all(v -> v isa Real, values(widths))) ||
+        throw(ArgumentError("[Tearing.IslandFlattening] widths must be a table of full widths in ψ_N keyed by m/n, e.g. widths = { \"2/1\" = 0.05 }"))
+    return Dict{String,Float64}(String(k) => Float64(v) for (k, v) in widths)
+end
+
 """
     perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple}
 
@@ -471,6 +488,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
             error("[Match] needs Δ′: use integrator = \"Galerkin\" with [Galerkin] rpec_flag = true, or \"Riccati\" with vac_flag = true")
     end
     tearing = tearing_from_deck(inputs)
+    island_flattening = island_flattening_from_deck(inputs)
     perturbed_equilibrium = perturbed_equilibrium_from_deck(inputs)
     forcing = forcing_terms_control(inputs, path)
     if haskey(inputs, "ErrorFields")
@@ -483,7 +501,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it")
     haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) &&
         error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it")
-    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing,
+    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, island_flattening,
         kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing,
         error_fields=error_fields_from_deck(inputs, path))
 end
@@ -578,6 +596,10 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an
         merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing),
             model, config))
     end
+    if deck.island_flattening !== nothing
+        tearing = get!(inputs, "Tearing", Dict{String,Any}())
+        tearing["IslandFlattening"] = Dict{String,Any}("widths" => Dict{String,Any}(deck.island_flattening))
+    end
     deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces))
     pe = deck.perturbed_equilibrium
     pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive)))
@@ -626,6 +648,8 @@ deck_sections(m::ForceFreeStates.MatchResult) =
     merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ"))))
 deck_sections(kf::KineticForcesResult) = deck_sections(kf.control)
 deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config)
+deck_sections(r::IslandFlatteningResult) =
+    Dict{String,Any}("Tearing" => Dict{String,Any}("IslandFlattening" => Dict{String,Any}("widths" => Dict{String,Any}(r.widths))))
 deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES)
 
 function deck_sections(ffs::ForceFreeStatesResult)
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index bda2b7c8d8..f451c77bd4 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -74,7 +74,7 @@ include("Rerun.jl")
 using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult
 using .ForceFreeStates: MatchProblem
 using .InnerLayer: GGJModel, SLAYERModel
-using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl
+using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl, IslandFlatteningResult, island_flattening
 using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult
 using .Equilibrium: PlasmaEquilibrium, KineticProfileData, IonSpecies, attach_kinetic_profiles!
 using .ForcingTerms: RMPField
@@ -140,7 +140,7 @@ equilibrium and writes nothing. The raw forcing file contents are captured in
 `Input/RawInputs/ForcingTerms/`, and `preloaded_coil_sets` lets the rerun path inject the coil
 geometry stored in `Input/RawInputs/Coils/`.
 
-Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each
+Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, island_flattening)`, each
 `nothing` when its section is absent. `equil` is the equilibrium the stages used (the solve's, on
 the two-pass grid), `ffs` the result the later stages read (the matched one when the deck has
 `[Match]`), and `match` its `MatchResult`.
@@ -158,7 +158,7 @@ function main_from_inputs(
     runtimes = Vector{Pair{String,Float64}}()
     deck = parse_deck(inputs, path)
     (; output_file, verbose) = deck.run
-    local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = nothing
+    local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = island_flat = nothing
 
     equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do
         Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
@@ -224,12 +224,20 @@ function main_from_inputs(
         save_to_hdf5(output_file, tearing)
     end
 
+    # After the tearing save, which rewrites Tearing/.
+    if deck.island_flattening !== nothing
+        island_flat = _timed_stage(runtimes, "island_flattening", "Island Flattening") do
+            island_flattening(equil; widths=deck.island_flattening)
+        end
+        save_to_hdf5(output_file, island_flat)
+    end
+
     total_dt = time() - total_start
     push!(runtimes, "total" => total_dt)
     ffs === nothing || _write_runtimes!(output_file, runtimes)
     @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
 
-    return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)
+    return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, island_flattening=island_flat)
 end
 
 """
@@ -306,6 +314,7 @@ export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fiel
 export KineticProfileData, IonSpecies, attach_kinetic_profiles!
 export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse
 export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel
+export IslandFlatteningResult, island_flattening
 export LocalStabilityProblem, Ballooning, LocalStabilityResult
 export KineticForcesControl, KineticForcesResult, ErrorFieldsResult
 
diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl
index dbf1ed5b9f..6418cc7702 100644
--- a/src/Tearing/HDF5Output.jl
+++ b/src/Tearing/HDF5Output.jl
@@ -1,9 +1,9 @@
 import ..Utilities: save_to_hdf5
-import ..Utilities.HDF5Annotations: replace_group!, write_ragged!
+import ..Utilities.HDF5Annotations: replace_group!, ensure_group, write_ragged!
 
 # HDF5Output.jl
 #
-# Writer for a `TearingResult`. Layout under the caller's parent group (inputs live only under Input/):
+# Writers for a `TearingResult` and an `IslandFlatteningResult`. Layout under the caller's parent group (inputs live only under Input/):
 #
 #   Tearing/
 #   ├── PerSurface/            -- layer parameters (SLAYER or GGJ set), Δ' matrix, omega_E, q_shift
@@ -12,7 +12,8 @@ import ..Utilities.HDF5Annotations: replace_group!, write_ragged!
 #   ├── Diagnostics/           -- ValidRoots, Poles, FilteredRoots (+ *Winding), flat-plus-offsets ragged;
 #   │                             uncoupled entries are in each surface's E×B frame
 #   ├── Scan/                  -- optional: full Q/Δ scan data
-#   └── CriticalResonantField/ -- optional: torque-balance critical b_r/B_φ per surface
+#   ├── CriticalResonantField/ -- optional: torque-balance critical b_r/B_φ per surface
+#   └── IslandFlattening/      -- its own result, written after the TearingResult in a deck run
 
 using HDF5
 
@@ -366,3 +367,61 @@ function _write_single_scan!(g, data::AMRResult)
     g["truncated"] = Int(data.truncated)
     return nothing
 end
+
+# ---------- island flattening ----------
+const ISLAND_FLATTENING_H5_ANNOTATIONS = [
+    "Islands/rational_m" => (; long_name="poloidal mode number m of each island", dims=("island",)),
+    "Islands/rational_n" => (; long_name="toroidal mode number n of each island", dims=("island",)),
+    "Islands/rational_psi" => (; long_name="rational surface ψ_N of each island", dims=("island",)),
+    "Islands/psi_in" => (; long_name="inner end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)),
+    "Islands/psi_out" => (; long_name="outer end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)),
+    "Islands/clipped" => (; long_name="flag: the island reaches past an end of the profile grid", dims=("island",)),
+    "Intervals/psi_in" => (; long_name="inner end of each merged flattened interval (ψ_N)", dims=("interval",)),
+    "Intervals/psi_out" => (; long_name="outer end of each merged flattened interval (ψ_N)", dims=("interval",)),
+    "Intervals/Delta_p" => (; long_name="pressure drop p(ψ_in) − p(ψ_out) across each merged interval", units="Pa", dims=("interval",)),
+    "Profiles/psi" => (; long_name="normalized poloidal flux ψ_N (profile grid with the interval ends added)", dims=("psi",)),
+    "Profiles/p" => (; long_name="equilibrium pressure", units="Pa", dims=("psi",)),
+    "Profiles/p_mod" => (; long_name="pressure flattened across the islands, anchored at the boundary", units="Pa", dims=("psi",)),
+    "Profiles/jbs_B" => (; long_name="given bootstrap profile (⟨j_bs·B⟩, or a current density for a current fraction)", dims=("psi",)),
+    "Profiles/jbs_B_mod" => (; long_name="given bootstrap profile with the flattened intervals removed (j_bs ∝ p′ to leading order)", dims=("psi",)),
+    "W" => (; long_name="stored energy (3/2)∫p dV", units="J"),
+    "Delta_W" => (; long_name="stored energy removed by the flattening (3/2)∫(p − p_mod) dV", units="J"),
+    "Delta_W_over_W" => (; long_name="fraction of the stored energy removed"),
+    "Delta_p0" => (; long_name="on-axis pressure drop", units="Pa"),
+    "jbs_fraction_lost" => (; long_name="fraction of the volume integral of the given bootstrap profile removed (NaN without one)")
+]
+
+"""
+    save_to_hdf5(parent, result::IslandFlatteningResult)
+
+Write `result` into `Tearing/IslandFlattening/` of `parent`, replacing any earlier one and
+keeping the rest of `Tearing/`.
+"""
+function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, r::IslandFlatteningResult)
+    g = replace_group!(ensure_group(parent, "Tearing"), "IslandFlattening")
+    isl = r.islands
+    g["Islands/rational_m"] = [i.m for i in isl]
+    g["Islands/rational_n"] = [i.n for i in isl]
+    g["Islands/rational_psi"] = [i.psi for i in isl]
+    g["Islands/psi_in"] = [i.psi_in for i in isl]
+    g["Islands/psi_out"] = [i.psi_out for i in isl]
+    g["Islands/clipped"] = Int[i.clipped for i in isl]
+    g["Intervals/psi_in"] = first.(r.intervals)
+    g["Intervals/psi_out"] = last.(r.intervals)
+    g["Intervals/Delta_p"] = r.delta_p
+    g["Profiles/psi"] = r.psi
+    g["Profiles/p"] = r.p
+    g["Profiles/p_mod"] = r.p_mod
+    if !isempty(r.jbs_B)
+        g["Profiles/jbs_B"] = r.jbs_B
+        g["Profiles/jbs_B_mod"] = r.jbs_B_mod
+    end
+    g["W"] = r.W
+    g["Delta_W"] = r.delta_W
+    g["Delta_W_over_W"] = r.delta_W / r.W
+    g["Delta_p0"] = r.delta_p0
+    g["jbs_fraction_lost"] = r.jbs_fraction_lost
+    Utilities.HDF5Annotations.annotate!(g, ISLAND_FLATTENING_H5_ANNOTATIONS)
+    Utilities.record_inputs!(parent, r)
+    return g
+end
diff --git a/src/Tearing/IslandFlattening/Flattening.jl b/src/Tearing/IslandFlattening/Flattening.jl
new file mode 100644
index 0000000000..9e1ca90ec5
--- /dev/null
+++ b/src/Tearing/IslandFlattening/Flattening.jl
@@ -0,0 +1,169 @@
+# Flattening.jl
+#
+# The pressure gradient vanishes inside each island and is unchanged outside it; the
+# pressure keeps its boundary value at the last profile point and is integrated inward,
+# as in the Fortran GPEC saturated-island profile. Islands that overlap flatten their
+# union once. With the flattened intervals [a_k, b_k] merged, the integral is exact:
+#
+#   p_mod(ψ) = p(ψ) + Σ_{k: b_k > ψ} [p(b_k) − p(max(ψ, a_k))],
+#
+# so every point inward of an interval drops by p(a_k) − p(b_k), and points inside it sit
+# at p(b_k). The stored energy removed is ΔW = (3/2) ∫ (p − p_mod) dV.
+
+"""
+    FlattenedIsland
+
+One island of the set: mode numbers `m`, `n`, rational surface `psi` (ψ_N), full width `w`
+(ψ_N) and the flattened interval `psi_in ≤ ψ ≤ psi_out`, clipped to the profile grid.
+`clipped` is true when the island reaches past either end of the grid.
+"""
+struct FlattenedIsland
+    m::Int
+    n::Int
+    psi::Float64
+    w::Float64
+    psi_in::Float64
+    psi_out::Float64
+    clipped::Bool
+end
+
+"""
+    IslandFlatteningResult
+
+Static flattening of the equilibrium pressure across `islands`.
+
+  - `widths`: the requested full widths keyed by `"m/n"` (ψ_N); every rational surface of each
+    `m/n` gets the width
+  - `islands`: one `FlattenedIsland` per rational surface found
+  - `intervals`: the flattened intervals, merged where islands overlap, inward to outward
+  - `delta_p`: pressure drop across each merged interval, p(a_k) − p(b_k) (Pa); for an edge
+    interval it is the pedestal-top pressure change
+  - `psi`, `p`, `p_mod`: the profile grid with the interval ends added, the equilibrium and
+    flattened pressures on it (Pa)
+  - `W`, `delta_W`: stored energy (3/2)∫p dV and the part the islands remove (J)
+  - `delta_p0`: on-axis pressure drop, the sum of `delta_p` (Pa)
+  - `jbs_B`, `jbs_B_mod`, `jbs_fraction_lost`: the given bootstrap profile on `psi`, the same
+    with the flattened intervals removed, and the fraction of its volume integral on them;
+    empty and NaN when no profile was given. To leading order j_bs ∝ p′, so it vanishes where
+    p′ does and is unchanged elsewhere. For ⟨j_bs·B⟩ the fraction is of ∫⟨j_bs·B⟩dV; for the
+    toroidal bootstrap current pass j_φ⟨1/R⟩, since dA = ⟨1/R⟩dV/2π.
+"""
+struct IslandFlatteningResult
+    widths::Dict{String,Float64}
+    islands::Vector{FlattenedIsland}
+    intervals::Vector{Tuple{Float64,Float64}}
+    delta_p::Vector{Float64}
+    psi::Vector{Float64}
+    p::Vector{Float64}
+    p_mod::Vector{Float64}
+    W::Float64
+    delta_W::Float64
+    delta_p0::Float64
+    jbs_B::Vector{Float64}
+    jbs_B_mod::Vector{Float64}
+    jbs_fraction_lost::Float64
+end
+
+"""
+    merged_intervals(islands) -> Vector{Tuple{Float64,Float64}}
+
+The union of the islands' flattened intervals as disjoint intervals, inward to outward.
+Intervals that touch or overlap merge.
+"""
+function merged_intervals(islands::AbstractVector{FlattenedIsland})
+    spans = sort!([(i.psi_in, i.psi_out) for i in islands if i.psi_out > i.psi_in])
+    merged = Tuple{Float64,Float64}[]
+    for (a, b) in spans
+        if !isempty(merged) && a <= merged[end][2]
+            merged[end] = (merged[end][1], max(merged[end][2], b))
+        else
+            push!(merged, (a, b))
+        end
+    end
+    return merged
+end
+
+"""
+    flattened_pressure(p, intervals, psi) -> Float64
+
+The flattened pressure at `psi` for a pressure profile `p(ψ)` and disjoint flattened
+`intervals`, anchored at the boundary value outside every interval.
+"""
+function flattened_pressure(p, intervals, psi::Real)
+    pm = p(psi)
+    for (a, b) in intervals
+        b > psi && (pm += p(b) - p(max(psi, a)))
+    end
+    return pm
+end
+
+"""
+    island_flattening(equil; widths, jbs_B=nothing) -> IslandFlatteningResult
+
+Flatten the pressure of `equil` across an island of full width `widths["m/n"]` (ψ_N) at
+every rational surface q = m/n (m, n > 0). `jbs_B`, a function of ψ_N giving a bootstrap
+profile, adds the part the flattening removes. The pressure is anchored at the last point of
+the equilibrium profile grid, so for an island reaching the edge the result depends on where
+that grid ends.
+"""
+function island_flattening(equil::PlasmaEquilibrium; widths::AbstractDict, jbs_B=nothing)
+    prof = equil.profiles
+    lo, hi = prof.xs[1], prof.xs[end]
+    p(psi) = prof.P_spline(psi) / mu0
+    dV(psi) = prof.dVdpsi_spline(psi)
+
+    requested = Dict{String,Float64}()
+    islands = FlattenedIsland[]
+    for (label, w) in widths
+        m, n = _mode_numbers(label)
+        w >= 0 || throw(ArgumentError("island width for $label must be non-negative, got $w"))
+        requested[label] = w
+        found = [s.psifac for s in _find_rational_surfaces(equil, n, n) if s.m == m]
+        isempty(found) && @warn "No q = $label surface in the equilibrium; island skipped"
+        for psi in found
+            a, b = psi - w / 2, psi + w / 2
+            push!(islands, FlattenedIsland(m, n, psi, w, max(a, lo), min(b, hi), a < lo || b > hi))
+        end
+    end
+    sort!(islands; by=i -> i.psi)
+    intervals = merged_intervals(islands)
+    delta_p = [p(a) - p(b) for (a, b) in intervals]
+
+    grid = sort!(unique!(vcat(prof.xs, (x for iv in intervals for x in iv)...)))
+    p_eq = p.(grid)
+    p_mod = [flattened_pressure(p, intervals, x) for x in grid]
+
+    # Exact for the piecewise-cubic integrands: nodes on every knot and interval end.
+    x, wq = gausslegendre(4)
+    segment(f, a, b) = (b - a) / 2 * sum(wq[j] * f((b - a) / 2 * x[j] + (a + b) / 2) for j in eachindex(x))
+    integrate(f) = sum(segment(f, grid[i], grid[i+1]) for i in 1:(length(grid)-1))
+    W = 1.5 * integrate(psi -> p(psi) * dV(psi))
+    delta_W = 1.5 * integrate(psi -> (p(psi) - flattened_pressure(p, intervals, psi)) * dV(psi))
+
+    inside(psi) = any(a < psi < b for (a, b) in intervals)
+    if jbs_B === nothing
+        jbs, jbs_mod, lost = Float64[], Float64[], NaN
+    else
+        jbs = Float64.(jbs_B.(grid))
+        jbs_mod = [inside(xi) ? 0.0 : j for (xi, j) in zip(grid, jbs)]
+        total = integrate(psi -> jbs_B(psi) * dV(psi))
+        lost = sum(integrate_between(psi -> jbs_B(psi) * dV(psi), grid, a, b, segment) for (a, b) in intervals; init=0.0) / total
+    end
+    return IslandFlatteningResult(requested, islands, intervals, delta_p, grid, p_eq, p_mod, W, delta_W, sum(delta_p; init=0.0),
+        jbs, jbs_mod, lost)
+end
+
+# Integral of f over [a, b], whose ends are nodes of `grid`.
+function integrate_between(f, grid, a, b, segment)
+    i0, i1 = searchsortedfirst(grid, a), searchsortedfirst(grid, b)
+    return sum(segment(f, grid[i], grid[i+1]) for i in i0:(i1-1); init=0.0)
+end
+
+function _mode_numbers(label::AbstractString)
+    parts = split(label, "/")
+    ok = length(parts) == 2 && all(s -> tryparse(Int, strip(s)) !== nothing, parts)
+    ok || throw(ArgumentError("island widths are keyed by \"m/n\", e.g. \"2/1\"; got \"$label\""))
+    m, n = parse.(Int, strip.(parts))
+    (m > 0 && n > 0) || throw(ArgumentError("island \"$label\" needs positive m and n"))
+    return m, n
+end
diff --git a/src/Tearing/IslandFlattening/IslandFlattening.jl b/src/Tearing/IslandFlattening/IslandFlattening.jl
new file mode 100644
index 0000000000..9e9e3a2949
--- /dev/null
+++ b/src/Tearing/IslandFlattening/IslandFlattening.jl
@@ -0,0 +1,15 @@
+# Static pressure flattening across a set of magnetic islands of given full width, with the
+# stored energy and pressure each set of islands removes. No time evolution.
+
+module IslandFlattening
+
+using FastGaussQuadrature: gausslegendre
+
+using ...Equilibrium: PlasmaEquilibrium, mu0
+using ...ForceFreeStates: _find_rational_surfaces
+
+include("Flattening.jl")
+
+export FlattenedIsland, IslandFlatteningResult, island_flattening, merged_intervals, flattened_pressure
+
+end # module IslandFlattening
diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl
index 0055cb2739..d597b7ee8e 100644
--- a/src/Tearing/Tearing.jl
+++ b/src/Tearing/Tearing.jl
@@ -6,6 +6,8 @@
 #                  extraction (consumes any InnerLayerModel)
 #   CriticalResonantField -- real-Q torque-balance scan for the critical
 #                  resonant field of error-field penetration
+#   IslandFlattening -- static pressure flattening across islands of given
+#                  width, and the stored energy it removes
 #   TearingProblem / TearingConfig / TearingResult -- the public solve,
 #                  `solve(TearingProblem(ffs; ...), model; config)`, with one
 #                  method per inner-layer model (SLAYERModel, GGJModel)
@@ -39,9 +41,11 @@ using .InnerLayer.SLAYER: _slayer_profiles, _layer_zeff, _layer_ion_mass, _layer
 include("LayerOverlap.jl")
 include("Dispersion/Dispersion.jl")
 include("CriticalResonantField/CriticalResonantField.jl")
+include("IslandFlattening/IslandFlattening.jl")
 
 import .Dispersion as Dispersion
 import .CriticalResonantField as CriticalResonantField
+import .IslandFlattening as IslandFlattening
 using .Dispersion: SurfaceCoupling, surface_coupling,
     MultiSurfaceCoupling, multi_surface_coupling,
     ScanResult, brute_force_scan,
@@ -49,6 +53,7 @@ using .Dispersion: SurfaceCoupling, surface_coupling,
     MultiBoxAMRResult, multi_box_amr_scan, as_amr_result,
     GrowthRateResult, find_growth_rates
 using .CriticalResonantField: TorqueBalance, torque_balance_scan
+using .IslandFlattening: IslandFlatteningResult, island_flattening
 
 include("Config.jl")
 include("Result.jl")
@@ -57,12 +62,13 @@ include("SLAYERTearing.jl")
 include("GGJTearing.jl")
 include("HDF5Output.jl")
 
-export InnerLayer, Dispersion, CriticalResonantField
+export InnerLayer, Dispersion, CriticalResonantField, IslandFlattening
 export resistive_layer_overlap, layer_overlap_scan, LayerOverlapScan
 export TearingProblem, TearingConfig, TearingResult
 export run_slayer_from_inputs
 export save_to_hdf5
 export CriticalResonantFieldControl, CriticalResonantFieldScan, CriticalResonantFieldResult
 export run_critical_resonant_field
+export IslandFlatteningResult, island_flattening
 
 end # module Tearing
diff --git a/src/Utilities/Output.jl b/src/Utilities/Output.jl
index bf75c28f85..f255e717d1 100644
--- a/src/Utilities/Output.jl
+++ b/src/Utilities/Output.jl
@@ -28,13 +28,20 @@ function deck_sections end
 """
     record_inputs!(target, result)
 
-Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, replacing any section of
-the same name, so every saved result adds the deck it was computed from.
+Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, so every saved result adds
+the deck it was computed from. A result replaces its section's scalars and the sub-tables it states, and keeps recorded
+sub-tables it does not state (a section stated only as sub-tables keeps the recorded scalars too).
 """
 function record_inputs!(target::Union{HDF5.File,HDF5.Group}, result)
     path = "Input/gpec_toml_raw"
     deck = haskey(target, path) ? TOML.parse(read(target[path])) : Dict{String,Any}()
-    merge!(deck, deck_sections(result))
+    # A result owns its section's scalars and the sub-tables it states; recorded sub-tables it does not state stay.
+    for (k, v) in deck_sections(result)
+        old = get(deck, k, nothing)
+        deck[k] =
+            !(v isa AbstractDict && old isa AbstractDict) ? v :
+            all(x -> x isa AbstractDict, values(v)) ? merge(old, v) : merge(filter(p -> p.second isa AbstractDict, old), v)
+    end
     haskey(target, path) && HDF5.delete_object(target, path)
     target[path] = sprint(TOML.print, deck)
     return nothing
diff --git a/test/runtests_island_flattening.jl b/test/runtests_island_flattening.jl
new file mode 100644
index 0000000000..ae13052986
--- /dev/null
+++ b/test/runtests_island_flattening.jl
@@ -0,0 +1,100 @@
+using TOML
+using HDF5
+
+# Static pressure flattening across islands of given width, on the coarse TJ fixture equilibrium.
+@testset "island flattening" begin
+    GPEC = GeneralizedPerturbedEquilibrium
+    IF = GPEC.Tearing.IslandFlattening
+    template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example")
+    deck = TOML.parsefile(joinpath(template, "gpec.toml"))
+    equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(deck, template), GPEC.analytic_config(deck))
+    prof = equil.profiles
+    p(psi) = prof.P_spline(psi) / GPEC.Equilibrium.mu0
+    dV(psi) = prof.dVdpsi_spline(psi)
+    one = island_flattening(equil; widths=Dict("2/1" => 0.05))
+    psi_s = only(one.islands).psi
+    # Volume inside ψ_N = x, by the midpoint rule.
+    volume_inside(x; n=4000) = sum(dV(prof.xs[1] + (x - prof.xs[1]) * (k - 0.5) / n) for k in 1:n) * (x - prof.xs[1]) / n
+
+    @testset "the Fortran trapezoid rule converges to the closed form" begin
+        # Fortran GPEC: zero dp/dψ on the island, anchor at the edge, integrate inward on a grid.
+        grid = range(prof.xs[1], prof.xs[end]; length=200_001)
+        dp = [abs(x - psi_s) <= 0.025 ? 0.0 : prof.P_deriv(x) / GPEC.Equilibrium.mu0 for x in grid]
+        pm = similar(dp)
+        pm[end] = p(grid[end])
+        for i in (length(grid)-1):-1:1
+            pm[i] = pm[i+1] - 0.5 * (dp[i+1] + dp[i]) * step(grid)
+        end
+        @test pm[1] ≈ IF.flattened_pressure(p, one.intervals, grid[1]) rtol = 1e-5
+        dW = 1.5 * sum((p(x) - y) * dV(x) for (x, y) in zip(grid, pm)) * step(grid)
+        @test dW ≈ one.delta_W rtol = 1e-3
+    end
+
+    @testset "narrow islands: ΔW → 0 and ΔW ∝ w" begin
+        @test island_flattening(equil; widths=Dict("2/1" => 0.0)).delta_W == 0
+        small = [island_flattening(equil; widths=Dict("2/1" => w)).delta_W for w in (1e-3, 2e-3)]
+        @test small[2] / small[1] ≈ 2 rtol = 1e-2
+        # Leading order: the drop −p′w across the island lowers everything inward of it.
+        @test small[1] ≈ 1.5 * (-prof.P_deriv(psi_s) / GPEC.Equilibrium.mu0) * 1e-3 * volume_inside(psi_s) rtol = 2e-2
+    end
+
+    @testset "overlapping islands are flattened once" begin
+        # The 2/1 (ψ_N ≈ 0.48) and 5/2 (≈ 0.79) islands overlap, so they flatten one interval.
+        widths = Dict("2/1" => 0.4, "5/2" => 0.3)
+        two = island_flattening(equil; widths)
+        (a, b) = only(two.intervals)
+        @test (a, b) == (two.islands[1].psi_in, two.islands[2].psi_out)
+        alone = [island_flattening(equil; widths=Dict(k => v)).delta_W for (k, v) in widths]
+        @test two.delta_W < sum(alone)
+        n = 200_000
+        grid = range(prof.xs[1], prof.xs[end]; length=n + 1)
+        mids = (grid[1:end-1] .+ grid[2:end]) ./ 2
+        @test two.delta_W ≈ 1.5 * sum((p(x) - IF.flattened_pressure(p, [(a, b)], x)) * dV(x) for x in mids) * step(grid) rtol = 1e-6
+        @test two.delta_p0 ≈ only(two.delta_p) ≈ p(a) - p(b)
+        # Inside the island the pressure sits at its outer-edge value.
+        @test IF.flattened_pressure(p, two.intervals, (a + b) / 2) ≈ p(b)
+    end
+
+    @testset "edges, bootstrap fraction and bad input" begin
+        wide = island_flattening(equil; widths=Dict("2/1" => 4.0))
+        @test only(wide.islands).clipped
+        # The whole profile flattens to its boundary value.
+        @test wide.delta_W ≈ wide.W - 1.5 * p(prof.xs[end]) * volume_inside(prof.xs[end]) rtol = 1e-8
+        # A uniform ⟨j_bs·B⟩ loses the volume fraction of the flattened interval.
+        frac = island_flattening(equil; widths=Dict("2/1" => 0.05), jbs_B=_ -> 1.0)
+        (a, b) = only(frac.intervals)
+        @test frac.jbs_fraction_lost ≈ (volume_inside(b) - volume_inside(a)) / volume_inside(prof.xs[end]) rtol = 1e-3
+        @test all(iszero, frac.jbs_B_mod[a.<frac.psi.<b])
+        # A profile peaked off the island loses less than its volume share; checked by the midpoint rule.
+        g(x) = exp(-((x - 0.7) / 0.1)^2)
+        peaked = island_flattening(equil; widths=Dict("2/1" => 0.05), jbs_B=g)
+        mid(f, lo, hi; n=20_000) = sum(f(lo + (hi - lo) * (k - 0.5) / n) for k in 1:n) * (hi - lo) / n
+        @test peaked.jbs_fraction_lost ≈ mid(x -> g(x) * dV(x), a, b) / mid(x -> g(x) * dV(x), prof.xs[1], prof.xs[end]) rtol = 1e-6
+        @test peaked.jbs_fraction_lost < frac.jbs_fraction_lost
+        # An island reaching inside the grid's first point is clipped there.
+        core = only(island_flattening(equil; widths=Dict("2/1" => 2 * only(one.islands).psi + 0.1)).islands)
+        @test core.clipped && core.psi_in == prof.xs[1]
+        @test isnan(one.jbs_fraction_lost) && isempty(one.jbs_B)
+        @test_throws ArgumentError island_flattening(equil; widths=Dict("2-1" => 0.05))
+        @test_throws ArgumentError island_flattening(equil; widths=Dict("2/1" => -0.1))
+    end
+
+    @testset "HDF5 output keeps the rest of Tearing/" begin
+        mktempdir() do dir
+            path = joinpath(dir, "gpec.h5")
+            h5open(path, "w") do f
+                f["Tearing/enabled"] = 1
+                f["Input/gpec_toml_raw"] = sprint(TOML.print, Dict("Tearing" => Dict("inner_model" => "SLAYER")))
+            end
+            save_to_hdf5(path, one)
+            h5open(path, "r") do f
+                @test read(f["Tearing/enabled"]) == 1
+                @test read(f["Tearing/IslandFlattening/Delta_W"]) == one.delta_W
+                @test read(f["Tearing/IslandFlattening/Profiles/p_mod"]) == one.p_mod
+                recorded = TOML.parse(read(f["Input/gpec_toml_raw"]))["Tearing"]
+                @test recorded["inner_model"] == "SLAYER"
+                @test recorded["IslandFlattening"]["widths"] == Dict("2/1" => 0.05)
+            end
+        end
+    end
+end
diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl
index f38ada7aa3..1baa914476 100644
--- a/test/runtests_result_struct.jl
+++ b/test/runtests_result_struct.jl
@@ -45,8 +45,8 @@ include("display_check.jl")
             ret = GeneralizedPerturbedEquilibrium.main([dir])
 
             # One field per stage, nothing for the sections the deck does not hold.
-            @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing)
-            @test ret.equil === ret.ffs.equil && all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing))
+            @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing, :island_flattening)
+            @test ret.equil === ret.ffs.equil && all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing, ret.island_flattening))
             ffs = ret.ffs
             forward_result[] = ffs
             @test ffs isa FFS.ForceFreeStatesResult
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index 4f657bbd23..ba9e419471 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -645,5 +645,26 @@ include("display_check.jl")
                 @test minimum(raised.kinetic.total.method_results["fgar"].psi_grid) >= raised.pe.psi_grid[1]
             end
         end
+
+        @testset "island flattening" begin
+            # [Tearing.IslandFlattening] alone runs no layer solve; it is island_flattening on the solve's equilibrium.
+            widths = Dict{String,Any}("2/1" => 0.05, "3/1" => 0.02)
+            extra = Dict{String,Any}("Tearing" => Dict{String,Any}("IslandFlattening" => Dict{String,Any}("widths" => widths)))
+            mktempdir() do dir
+                reference = GPEC.main([_stage_deck(dir; extra_sections=extra)])
+                @test reference.tearing === nothing
+                api = island_flattening(reference.ffs.equil; widths)
+                @test api.delta_W == reference.island_flattening.delta_W
+                @test api.p_mod == reference.island_flattening.p_mod
+                h5open(joinpath(dir, "gpec.h5")) do f
+                    @test read(f["Tearing/IslandFlattening/Delta_W"]) == api.delta_W
+                    replay = GPEC.parse_deck(TOML.parse(read(f["Input/gpec_toml_raw"])), dir)
+                    @test replay.island_flattening == Dict{String,Float64}(widths)
+                    @test replay.tearing === nothing
+                end
+                @test_throws ArgumentError GPEC.island_flattening_from_deck(Dict("Tearing" => Dict("IslandFlattening" => Dict("widths" => 0.05))))
+                @test_throws ArgumentError GPEC.island_flattening_from_deck(Dict("Tearing" => Dict("IslandFlattening" => Dict("width" => widths))))
+            end
+        end
     end
 end
diff --git a/test/suites.jl b/test/suites.jl
index dfe79f9e4f..781be84d26 100644
--- a/test/suites.jl
+++ b/test/suites.jl
@@ -30,6 +30,7 @@ const TEST_SUITES = [
     "runtests_tj_analytic.jl" => "core",
     "runtests_kinetic_profiles.jl" => "core",
     "runtests_resist_eval.jl" => "core",
+    "runtests_island_flattening.jl" => "api",
     "runtests_localstability.jl" => "core",
     "runtests_analysis_plots.jl" => "api",
     "runtests_slayer_params.jl" => "errorfields",