diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 241049189..8a1c84b06 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -69,5 +69,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$
diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl
index e258fb457..828f1bbb0 100644
--- a/benchmarks/benchmark_q_vs_iota_edge.jl
+++ b/benchmarks/benchmark_q_vs_iota_edge.jl
@@ -20,11 +20,12 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp
 
 # Dense ldp reference equilibrium: treat its q(ψ) as ground truth
 function reference_q()
-    _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
-    eq_config.grid_type = "ldp"
-    eq_config.mpsi = 1024
+    inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
+    equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "ldp", "mpsi" => 1024))
+    eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR)
     equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
-    return equil, eq_config
+    # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable.
+    return equil, equil_dict
 end
 
 # Rational surface ψ_s(q = m/n) and q'(ψ_s) recovered from a fitted spline by bisection
@@ -44,7 +45,7 @@ function rational_locations(q_itp, dq_itp, psis, q_targets)
 end
 
 function main()
-    equil, eq_config = reference_q()
+    equil, equil_dict = reference_q()
     xs_ref = equil.profiles.xs
     q_ref_itp = equil.profiles.q_spline
     dq_ref_itp = equil.profiles.q_deriv
@@ -64,8 +65,7 @@ function main()
 
     for N in (16, 32, 64, 128, 256)
         # Same log_asymptotic knot layout the auto grid uses for fixed mpsi
-        cfg = deepcopy(eq_config)
-        cfg.mpsi = N
+        cfg = GPE.Equilibrium.EquilibriumConfig(merge(equil_dict, Dict{String,Any}("mpsi" => N)), EXAMPLE_DIR)
         knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh)
         q_nodes = [q_ref_itp(p) for p in knots]
 
diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml
index 7065439c5..f5d5be92b 100644
--- a/examples/DIIID-like_SLAYER_example/gpec.toml
+++ b/examples/DIIID-like_SLAYER_example/gpec.toml
@@ -14,7 +14,7 @@ psihigh = 0.9995               # Upper limit of normalized flux coordinate
 mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
 psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-10                   # Error tolerance for equilibrium solver
 force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
index cbb702db5..efa77bb14 100644
--- a/examples/DIIID-like_gal_resistive_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -14,7 +14,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-10                   # Error tolerance for equilibrium solver
 force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
index 82001c077..cff5472ff 100644
--- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -13,7 +13,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-10                   # Error tolerance for equilibrium solver
 force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 068432501..b84e12872 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -13,7 +13,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-10                   # Error tolerance for equilibrium solver
 force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml
index 14c192bd9..a54725908 100644
--- a/examples/Solovev_ideal_example/gpec.toml
+++ b/examples/Solovev_ideal_example/gpec.toml
@@ -11,7 +11,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index 4d9c89086..7a0e40200 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml
index 5286b5302..a82c3bcc0 100644
--- a/examples/Solovev_ideal_example_multi_n/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/gpec.toml
@@ -11,7 +11,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
index 8b2b75007..59760b231 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
index c7a51c723..a80ac7741 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml
index 8d3f90d3f..e52b0d537 100644
--- a/examples/Solovev_kinetic_NTV_example/gpec.toml
+++ b/examples/Solovev_kinetic_NTV_example/gpec.toml
@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid points (0 = auto-compute from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
index 760ba8e7d..7d3a13059 100644
--- a/examples/Solovev_kinetic_calculated_example/gpec.toml
+++ b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -11,7 +11,7 @@ psilow = 1e-4                           # Lower limit of normalized flux coordin
 psihigh = 0.9995                        # Upper limit of normalized flux coordinate
 mpsi = 16                               # Number of radial grid points
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml
index c6a6e9639..fd8881122 100644
--- a/examples/a10_kinetic_example/gpec.toml
+++ b/examples/a10_kinetic_example/gpec.toml
@@ -15,7 +15,7 @@ psilow = 1e-3                  # Lower limit of normalized poloidal flux
 psihigh = 0.99                 # Upper limit of normalized poloidal flux
 mpsi = 16                      # Number of radial grid points (low value for fast iteration)
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                    # Error tolerance for equilibrium solver
 
 [Wall]
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index ed6582131..961a0f7f5 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -478,7 +478,7 @@ robustness.
     psio = raw_profile.psio
     mtheta = equil_params.mtheta
     psilow = equil_params.psilow
-    psihigh = equil_params.psihigh
+    psihigh = raw_profile.psihigh_resolved
 
     # Locate the magnetic axis and separatrix for the field-line integrations
     ro, zo, _, rs2 = direct_position!(raw_profile)
@@ -545,15 +545,18 @@ robustness.
     if q0 <= 0.0
         @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use newq0 to override."
     end
-    if equil_params.newq0 == -1
-        equil_params.newq0 = -q0
+    # The -1 sentinel means "flip the extrapolated q0"; resolve it into a local so the
+    # config stays the user's request (matches equilibrium_solver(::InverseRunInput)).
+    newq0 = equil_params.newq0
+    if newq0 == -1
+        newq0 = -q0
     end
-    if equil_params.newq0 != 0.0
-        @info "Revising q-profile for newq0 = $(@sprintf("%.3f", equil_params.newq0))"
+    if newq0 != 0.0
+        @info "Revising q-profile for newq0 = $(@sprintf("%.3f", newq0))"
         f0 = profiles.F_spline.y[1] - profiles.F_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1]
-        f0fac = f0^2 * ((equil_params.newq0 / q0)^2 - 1.0)
+        f0fac = f0^2 * ((newq0 / q0)^2 - 1.0)
         for i in 1:(mpsi+1)
-            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(equil_params.newq0)
+            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(newq0)
             sq_nodes[i, 1] *= ffac
             sq_nodes[i, 4] *= ffac
             rzphi_nodes[i, :, 3] .*= ffac
@@ -651,6 +654,7 @@ robustness.
 
     params = EquilibriumParameters()
     params.bt_sign = raw_profile.bt_sign
+    params.psihigh_resolved = psihigh
 
     return PlasmaEquilibrium(raw_profile.config, params, profiles, geometry,
         rzphi_xs, rzphi_ys,
diff --git a/src/Equilibrium/DirectEquilibriumByInversion.jl b/src/Equilibrium/DirectEquilibriumByInversion.jl
index 252ae515b..136d8e6ac 100644
--- a/src/Equilibrium/DirectEquilibriumByInversion.jl
+++ b/src/Equilibrium/DirectEquilibriumByInversion.jl
@@ -73,12 +73,12 @@ end
 """
     clamp_psihigh_to_separatrix(raw_profile) -> (clamped_psihigh, was_adjusted)
 
-Binary-searches for the highest psihigh ≤ raw_profile.config.psihigh at which the
+Binary-searches for the highest psihigh ≤ `raw_profile.psihigh_resolved` at which the
 ψ level set is still a closed curve in the EFIT grid. Returns the safe value and a
 Bool indicating whether any clamping occurred.
 """
 function clamp_psihigh_to_separatrix(raw_profile::DirectRunInput)
-    psihigh = raw_profile.config.psihigh
+    psihigh = raw_profile.psihigh_resolved
     ψ_coarse = raw_profile.psi_in.nodal_derivs.partials[1, :, :]
 
     has_closed_contour(ψ_high) = any(
@@ -99,6 +99,22 @@ function clamp_psihigh_to_separatrix(raw_profile::DirectRunInput)
     return (lo, true)
 end
 
+"""
+    resolve_psihigh!(raw_profile) -> raw_profile
+
+Clamp `raw_profile.psihigh_resolved` to the outermost closed flux surface, warning if the
+requested `psihigh` had none. Mutates and returns the input; `config.psihigh` keeps the request.
+"""
+function resolve_psihigh!(raw_profile::DirectRunInput)
+    psihigh_safe, adjusted = clamp_psihigh_to_separatrix(raw_profile)
+    if adjusted
+        @warn "psihigh=$(raw_profile.psihigh_resolved) has no closed flux surface in EFIT grid; " *
+              "clamped to $(round(psihigh_safe; sigdigits=7))"
+        raw_profile.psihigh_resolved = psihigh_safe
+    end
+    return raw_profile
+end
+
 """
     resample_contour_to_theta_grid!(R_out, Z_out, curve, ro, zo, theta_grid)
 
@@ -408,7 +424,7 @@ function equilibrium_solver_by_inversion(
     psio = raw_profile.psio
     mtheta = equil_params.mtheta
     psilow = equil_params.psilow
-    psihigh = equil_params.psihigh
+    psihigh = raw_profile.psihigh_resolved
 
     # Locate the magnetic axis and separatrix for the contour tracing
     ro, zo, _, rs2 = direct_position!(raw_profile)
@@ -672,7 +688,7 @@ function equilibrium_solver_by_inversion(
     # Intermediate inverse input; the captured DirectIngest rides on the original eq_input,
     # which setup_equilibrium forwards onto the equilibrium, so this one carries ingest=nothing.
     inv_input = InverseRunInput(raw_profile.config, raw_profile.sq_in,
-        rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing)
+        rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing, raw_profile.psihigh_resolved)
 
     pe = equilibrium_solver(inv_input; override_psi_nodes)
 
diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl
index cf08d8a4f..6070e3359 100644
--- a/src/Equilibrium/Equilibrium.jl
+++ b/src/Equilibrium/Equilibrium.jl
@@ -92,27 +92,16 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin
     if additional_input isa DirectRunInput
         eq_input = additional_input
         eq_input.config = eq_config
-        # Re-run the separatrix clamp for efit-family replays so an overridden
-        # psihigh from the rerun TOML is re-validated against the closed flux region.
-        if eq_type in EFIT_KINDS
-            psihigh_safe, adjusted = clamp_psihigh_to_separatrix(eq_input)
-            if adjusted
-                @warn "psihigh=$(eq_input.config.psihigh) has no closed flux surface in EFIT grid; " *
-                      "clamped to $(round(psihigh_safe; sigdigits=7))"
-                eq_input.config.psihigh = psihigh_safe
-            end
-        end
+        # Reset before re-resolving: a pass-1 clamped value must not shadow a psihigh
+        # overridden in the rerun TOML.
+        eq_input.psihigh_resolved = eq_config.psihigh
+        eq_type in EFIT_KINDS && resolve_psihigh!(eq_input)
     elseif additional_input isa InverseRunInput
         eq_input = additional_input
         eq_input.config = eq_config
+        eq_input.psihigh_resolved = eq_config.psihigh
     elseif eq_type in EFIT_KINDS
-        eq_input = read_efit(eq_config)
-        psihigh_safe, adjusted = clamp_psihigh_to_separatrix(eq_input)
-        if adjusted
-            @warn "psihigh=$(eq_input.config.psihigh) has no closed flux surface in EFIT grid; " *
-                  "clamped to $(round(psihigh_safe; sigdigits=7))"
-            eq_input.config.psihigh = psihigh_safe
-        end
+        eq_input = resolve_psihigh!(read_efit(eq_config))
     elseif eq_type in ["chease2", "chease_ascii"]
         eq_input = read_chease_ascii(eq_config)
     elseif eq_type in ["chease", "chease_binary"]
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index 4bb29f8f2..344f5a60f 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -5,8 +5,8 @@ end
 """
     EquilibriumConfig(...)
 
-A mutable struct containing configuration parameters for equilibrium reconstruction.
-Bundles all necessary settings originally specified in the equil fortran namelists.
+An immutable struct containing configuration parameters for equilibrium reconstruction
+specified in the input.
 
 ## Fields
 
@@ -27,7 +27,8 @@ Bundles all necessary settings originally specified in the equil fortran namelis
     refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform";
     "log_asymptotic" is a legacy alias for "auto")
   - `psilow::Float64` - Lower limit of normalized flux coordinate
-  - `psihigh::Float64` - Upper limit of normalized flux coordinate
+  - `psihigh::Float64` - Requested upper limit of normalized flux coordinate; the value the
+    equilibrium is actually formed on is `EquilibriumParameters.psihigh_resolved`.
   - `mpsi::Int` - Number of radial grid intervals; 0 with grid_type="auto" selects the
     two-pass auto grid: the main driver forms a coarse pass-1 equilibrium, measures its curvature,
     pins knots on rational surfaces, and re-forms on the refined grid. Standalone `setup_equilibrium`
@@ -35,12 +36,13 @@ Bundles all necessary settings originally specified in the equil fortran namelis
   - `psi_accuracy::Float64` - Target relative accuracy τ of splined profile derivatives for the
     two-pass auto grid (knot count scales as τ^(-1/3))
   - `mtheta::Int` - Number of poloidal grid points
-  - `newq0::Int` - Override for on-axis safety factor (0 = use input value)
+  - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
+    meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
   - `etol::Float64` - Error tolerance for equilibrium solver
   - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations)
   - `use_galgrid::Bool` - Use the same grid as galerkin method
 """
-@kwdef mutable struct EquilibriumConfig
+@kwdef struct EquilibriumConfig
     eq_type::String = "efit"
     eq_filename::String = "mypath"
     r0exp::Float64 = 1.0
@@ -64,7 +66,7 @@ Bundles all necessary settings originally specified in the equil fortran namelis
     psi_accuracy::Float64 = 0.001
     mtheta::Int = 512
 
-    newq0::Int = 0
+    newq0::Float64 = 0.0
     etol::Float64 = 1e-10
 
     force_termination::Bool = false
@@ -176,15 +178,14 @@ function EquilibriumConfig(equil_dict::Dict{String,Any}, base_path::String="./")
         end
     end
 
-    # Construct validated struct
-    config = EquilibriumConfig(; symbolize_keys(config_data)...)
-    # Only resolve `eq_filename` against `base_path` if the user actually
-    # supplied one (otherwise leave the kwdef sentinel for the embedded path).
-    if haskey(config_data, "eq_filename") && !isabspath(config.eq_filename)
-        config.eq_filename = normpath(joinpath(base_path, config.eq_filename))
+    # Only resolve `eq_filename` against `base_path` if the user actually supplied one
+    # (otherwise leave the kwdef sentinel for the embedded path). The empty string is the
+    # rerun path's "no input file" marker and must stay empty, not become `base_path`.
+    if haskey(config_data, "eq_filename") && !isempty(config_data["eq_filename"]) && !isabspath(config_data["eq_filename"])
+        config_data["eq_filename"] = normpath(joinpath(base_path, config_data["eq_filename"]))
     end
 
-    return config
+    return EquilibriumConfig(; symbolize_keys(config_data)...)
 end
 
 """
@@ -208,12 +209,7 @@ function EquilibriumConfig(path::String)
     end
 
     # Construct validated struct
-    config = EquilibriumConfig(; symbolize_keys(config_data)...)
-    if !isabspath(config.eq_filename)
-        config.eq_filename = normpath(joinpath(dirname(path), config.eq_filename))
-    end
-
-    return config
+    return EquilibriumConfig(Dict{String,Any}(config_data), dirname(path))
 end
 
 """
@@ -444,6 +440,12 @@ raw equilibrium data and preparing the initial splines.
   - `bt_sign::Int` — Sign of the toroidal field (+1 or -1); read from fpol sign in EFIT g-files
   - `ingest::EquilibriumIngest` — captured raw arrays for the `gpec.h5` rerun snapshot
     (a [`DirectIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria)
+  - `psihigh_resolved::Float64` — outer flux limit the equilibrium is formed on: `config.psihigh`
+    clamped to the outermost closed flux surface by [`resolve_psihigh!`](@ref). Defaults to
+    `config.psihigh` and only differs for efit-family equilibria whose requested limit falls
+    outside the closed-flux region. The solvers build their ψ grid from this field. IMAS
+    equilibria are read into this struct but are not in `EFIT_KINDS`, so they are never
+    clamped and always keep the request.
 """
 mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND}
     config::EquilibriumConfig
@@ -458,8 +460,16 @@ mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:
     psio::Float64    # The total flux difference |ψ_axis - ψ_boundary| [Weber / radian].
     bt_sign::Int     # Sign of the toroidal field: +1 or -1 (from fpol sign in g-file)
     ingest::EquilibriumIngest
+    psihigh_resolved::Float64
 end
 
+# Readers construct without a resolved psihigh; it starts at the request and `resolve_psihigh!`
+# clamps it for efit-family equilibria.
+DirectRunInput(config::EquilibriumConfig, sq_in, psi_in, psi_in_xs, psi_in_ys,
+    rmin, rmax, zmin, zmax, psio, bt_sign, ingest) =
+    DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys,
+        rmin, rmax, zmin, zmax, psio, bt_sign, ingest, config.psihigh)
+
 """
     InverseRunInput(...)
 
@@ -478,6 +488,9 @@ A container struct for inputs to the `inverse_run` function.
   - `psio::Float64` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad]
   - `ingest::EquilibriumIngest` - captured raw arrays for the `gpec.h5` rerun snapshot
     (an [`InverseIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria)
+  - `psihigh_resolved::Float64` - outer flux limit the equilibrium is formed on; see
+    [`DirectRunInput`](@ref). Equals `config.psihigh` for every inverse reader (CHEASE,
+    analytic); `efit_by_inversion` forwards the clamped value from its `DirectRunInput`.
 """
 mutable struct InverseRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND}
     config::EquilibriumConfig
@@ -490,8 +503,14 @@ mutable struct InverseRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<
     zo::Float64                 # Z axis location
     psio::Float64               # Total flux difference |psi_axis - psi_boundary|
     ingest::EquilibriumIngest
+    psihigh_resolved::Float64
 end
 
+InverseRunInput(config::EquilibriumConfig, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z,
+    ro, zo, psio, ingest) =
+    InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z,
+        ro, zo, psio, ingest, config.psihigh)
+
 """
     EquilibriumParameters
 
@@ -502,6 +521,8 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags
   - `ro::Union{Nothing,Float64}` - R-coordinate of the magnetic axis [m]
   - `zo::Union{Nothing,Float64}` - Z-coordinate of the magnetic axis [m]
   - `psio::Union{Nothing,Float64}` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad]
+  - `psihigh_resolved::Union{Nothing,Float64}` - Outer flux limit the equilibrium was formed on
+    (the outermost ψ node); the plasma edge downstream of `setup_equilibrium`.
   - `rsep::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma boundary [m]
   - `zsep::Union{Nothing,Vector{Float64}}` - Z-coordinates of the plasma boundary [m]
   - `rext::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma edge [m]
@@ -554,6 +575,7 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags
     ro::Union{Nothing,Float64} = nothing # R-coordinate of the magnetic axis [m]
     zo::Union{Nothing,Float64} = nothing # Z-coordinate of the magnetic axis [m]
     psio::Union{Nothing,Float64} = nothing # Total flux difference |ψ_axis - ψ_boundary| [Wb/rad]
+    psihigh_resolved::Union{Nothing,Float64} = nothing # Outer flux limit actually formed on (clamped config.psihigh)
     rsep::Union{Nothing,Vector{Float64}} = nothing # R-coordinates of the plasma boundary [m]
     zsep::Union{Nothing,Vector{Float64}} = nothing # Z-coordinates of the plasma boundary [m]
     rext::Union{Nothing,Vector{Float64}} = nothing # R-coordinates of the plasma edge [m]
diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl
index 725391b15..4a8d48bb6 100644
--- a/src/Equilibrium/InverseEquilibrium.jl
+++ b/src/Equilibrium/InverseEquilibrium.jl
@@ -59,7 +59,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
     mpsi = config.mpsi
     mtheta = config.mtheta
     psilow = config.psilow
-    psihigh = config.psihigh
+    psihigh = input.psihigh_resolved
     newq0 = config.newq0
 
     # c-----------------------------------------------------------------------
@@ -393,7 +393,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
 
     return PlasmaEquilibrium(
         input.config,
-        EquilibriumParameters(),
+        EquilibriumParameters(; psihigh_resolved=psihigh),
         profiles,
         geometry,
         rzphi_xs, rzphi_ys,
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index 50eed7da0..69fdd5f98 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -240,7 +240,6 @@ gpec.toml.
   - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (PENTRC — not yet implemented)
   - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components.
   - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil
-  - `reform_eq_with_psilim::Bool` - Reform equilibrium with computed psilim (not yet implemented)
   - `psiedge::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [psiedge, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below).
   - `truncate_at_dW_peak::Bool` - When `true` and `psiedge < psilim`, the edge-dW scan's peak location is adopted as the new physical plasma edge — `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, AND the FM Δ' chunks/propagators are made self-consistent with the new boundary (the chunk that straddles the peak is rebuilt + re-integrated; any chunks past the peak are dropped). This reproduces the spirit of the original ode_record_edge heuristic from Fortran STRIDE while keeping Δ' and δW well-defined at the new boundary. The Δ' metric is still physically dependent on where the peak falls in the edge band, so use this flag deliberately when you mean to scan against the peak-defined edge (e.g. for studying edge-mode regimes); leave at `false` (default) for the full-domain Δ' at `qhigh` / `psihigh` / `dmlim`.
   - `diagnose::Bool` - Enable diagnostic output (not yet implemented)
@@ -278,7 +277,6 @@ gpec.toml.
     kinetic_source::String = "fixed"
     kinetic_factor::Float64 = 0.0
     qlow::Float64 = 0.0
-    reform_eq_with_psilim::Bool = false
     psiedge::Float64 = 0.99
     truncate_at_dW_peak::Bool = false   # Edge-dW peak becomes new physical edge; Δ' BVP made self-consistent. See docstring.
     diagnose::Bool = false
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 13d400128..f41e6dd2f 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -115,7 +115,7 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     # Initial guesses based on equilibrium
     intr.qlim = min(equil.params.qmax, ctrl.qhigh) # equilibrium solve only goes up to qmax, so we're capped there
     intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1))
-    intr.psilim = equil.config.psihigh
+    intr.psilim = equil.params.psihigh_resolved
 
     # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). The cutoff reads
     # the *resolved* toroidal range on `intr`, so callers must assign intr.nlow / intr.nhigh
diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl
index de4e16801..3bd938b43 100644
--- a/src/ForcingTerms/ForcingTerms.jl
+++ b/src/ForcingTerms/ForcingTerms.jl
@@ -28,7 +28,7 @@ Coil settings (used when `forcing_data_format = "coil"`):
   - `nzeta_coil::Int` - Toroidal grid resolution; 0 = auto (32 × n)
   - `coil_sets_raw::Vector{Dict{String,Any}}` - Parsed `[[ForcingTerms.coil_set]]` TOML blocks
 """
-Base.@kwdef mutable struct ForcingTermsControl
+@kwdef struct ForcingTermsControl
     # Forcing data file settings (ascii/hdf5 formats)
     forcing_data_file::String = "forcing.dat"
     forcing_data_format::String = "ascii"
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 0de2b07a5..7b17683cf 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -75,7 +75,7 @@ using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions!
 using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate
 
-const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag",
+const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim",
                               "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc")
 
@@ -327,15 +327,6 @@ function main_from_inputs(
     # Determine psilim and qlim (where we will integrate to)
     sing_lim!(intr, ctrl, equil)
 
-    # If truncating before psihigh, reform equilibrium if desired
-    if intr.psilim != equil.config.psihigh && ctrl.reform_eq_with_psilim
-        @warn "Reforming equilibrium splines from psihigh to psilim not implemented yet. Proceeding with psihigh = $(equil.config.psihigh)."
-        # JMH - Nik please put the logic we discussed here
-        # something like ?
-        # equil.config.psihigh = intr.psilim
-        # equil = set_up_equilibrium(equil.config)
-    end
-
     # Compute local stability (if desired). `locstab` holds `D_I` from the ballooning
     # coefficient system and the local ballooning result; `nothing` when not computed.
     locstab = nothing
@@ -567,13 +558,12 @@ function main_from_inputs(
         if "ForcingTerms" in keys(inputs)
             forcing_raw = inputs["ForcingTerms"]
             # [[ForcingTerms.coil_set]] becomes a Vector{Dict} — must be excluded from
-            # kwarg splatting and handled separately via coil_sets_raw field
+            # kwarg splatting and passed as the explicit coil_sets_raw keyword
             coil_sets_raw = Vector{Dict{String,Any}}(get(forcing_raw, "coil_set", Dict{String,Any}[]))
             scalar_forcing = filter(p -> p.first != "coil_set", forcing_raw)
             ft_ctrl = ForcingTerms.ForcingTermsControl(;
-                (Symbol(k) => v for (k, v) in scalar_forcing)...
+                (Symbol(k) => v for (k, v) in scalar_forcing)..., coil_sets_raw=coil_sets_raw
             )
-            ft_ctrl.coil_sets_raw = coil_sets_raw
         else
             ft_ctrl = ForcingTerms.ForcingTermsControl()  # Use defaults
         end
diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl
index a8beacffa..a1a47ef44 100644
--- a/src/HDF5Schema.jl
+++ b/src/HDF5Schema.jl
@@ -65,6 +65,7 @@ const MAIN_H5_ANNOTATIONS = [
     "Equilibrium/R_axis" => (; long_name="R-coordinate of the magnetic axis", units="m"),
     "Equilibrium/Z_axis" => (; long_name="Z-coordinate of the magnetic axis", units="m"),
     "Equilibrium/psi_total" => (; long_name="total poloidal flux difference |ψ_axis − ψ_boundary|", units="Wb/rad"),
+    "Equilibrium/psihigh_resolved" => (; long_name="normalized poloidal flux at the outermost formed flux surface (requested psihigh clamped to the outermost closed surface)"),
     "Equilibrium/R_midplane" => (; long_name="R of the boundary at the inboard and outboard midplane crossings", units="m"),
     "Equilibrium/R_extremum" => (; long_name="R of the boundary at its upper and lower Z extrema", units="m"),
     "Equilibrium/Z_extremum" => (; long_name="Z of the boundary at its upper and lower Z extrema", units="m"),
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 6df40e6a8..0d1b46464 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -60,7 +60,7 @@ Constructed via keyword arguments or from a TOML dict:
 ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...)
 ```
 """
-@kwdef mutable struct KineticForcesControl
+@kwdef struct KineticForcesControl
     # Moment type
     moment::String = "pressure"     # "heat" or "pressure"
 
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
index 96f982386..1ce0895d6 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
@@ -29,7 +29,7 @@ Regularization:
     # High Priority (MWE)
   - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0.
 """
-@kwdef mutable struct PerturbedEquilibriumControl
+@kwdef struct PerturbedEquilibriumControl
     # High Priority (MWE)
     fixed_boundary::Bool = false
     output_eigenmodes::Bool = true
diff --git a/src/Rerun.jl b/src/Rerun.jl
index 7d932addf..03c19bea6 100644
--- a/src/Rerun.jl
+++ b/src/Rerun.jl
@@ -267,9 +267,10 @@ function build_inputs_from_h5(args::Vector{String})
           "  output: $(abspath(joinpath(output_dir, output_name)))\n$_BANNER"
 
     _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium")
-    eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], output_dir)
-    # Clear eq_filename: unused on replay, and a stale absolute path could mislead downstream code.
-    eq_config.eq_filename = ""
+    # Clear eq_filename on a copy: unused on replay, a stale absolute path could mislead
+    # downstream code, and `inputs` itself is re-serialized into the rerun's gpec_toml_raw.
+    equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("eq_filename" => ""))
+    eq_config = Equilibrium.EquilibriumConfig(equil_dict, output_dir)
 
     # Analytic kinds regenerate from their TOML section; file-based kinds rebuild splines from
     # the stored ingest. A file-based run with no ingest can only come from a pre-ingest gpec.h5.
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index dd1e49ce5..af3cd5743 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -68,6 +68,78 @@
         @test all(>(0), B_nodes)
     end
 
+    @testset "Resolved psihigh" begin
+        # The config holds the user's request and is never written to; the value the
+        # equilibrium is actually formed on rides on params.psihigh_resolved.
+        for eq in (plasma_eq_efit, plasma_eq_arclength, plasma_eq_inversion)
+            @test eq.params.psihigh_resolved == eq.rzphi_xs[end]
+            @test eq.params.psihigh_resolved <= eq.config.psihigh
+        end
+        # 0.994 sits inside the closed-flux region, so nothing is clamped
+        @test plasma_eq_efit.params.psihigh_resolved == 0.994
+
+        # Requesting the separatrix itself must leave the request intact on the config
+        edge_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(;
+            eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"),
+            eq_type="efit",
+            jac_type="boozer",
+            grid_type="ldp",
+            psilow=0.01,
+            psihigh=1.0
+        )
+        plasma_eq_edge = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(edge_config)
+        @test edge_config.psihigh == 1.0
+        @test plasma_eq_edge.params.psihigh_resolved <= 1.0
+        @test plasma_eq_edge.params.psihigh_resolved == plasma_eq_edge.rzphi_xs[end]
+    end
+
+    @testset "newq0 q-profile revision" begin
+        # newq0 is a target q(0), not an index: non-integer requests must survive the config,
+        # and neither the -1 sentinel nor an explicit target may write back to the frozen config.
+        newq0_config(newq0) = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(;
+            eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"),
+            eq_type="efit",
+            jac_type="boozer",
+            grid_type="ldp",
+            mpsi=32,
+            psilow=0.01,
+            psihigh=0.994,
+            newq0=newq0
+        )
+
+        base_config = newq0_config(0)
+        eq_base = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(base_config)
+
+        # -1 means "flip the sign of the axis extrapolation": f0fac vanishes, so ffac is
+        # exactly -1 and the revised q-profile is the negated baseline.
+        flip_config = newq0_config(-1)
+        eq_flip = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(flip_config)
+        @test flip_config.newq0 == -1.0
+        @test eq_flip.profiles.q_spline.y == -eq_base.profiles.q_spline.y
+
+        # A non-integer target used to throw InexactError at construction (newq0 was ::Int)
+        target_config = newq0_config(1.05)
+        @test target_config.newq0 == 1.05
+        eq_target = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(target_config)
+        @test all(isfinite, eq_target.profiles.q_spline.y)
+        # Extrapolated q(0) meets the target; not exact because ffac uses the first node's F
+        # rather than the extrapolated axis value.
+        xs = eq_target.profiles.xs
+        q_axis = eq_target.profiles.q_spline.y[1] - eq_target.profiles.q_deriv(xs[1]; hint=Ref(1)) * xs[1]
+        @test isapprox(q_axis, 1.05; rtol=0.02)
+    end
+
+    @testset "Deprecated TOML keys are dropped, not fatal" begin
+        # Removed control knobs must keep old gpec.toml decks (and older gpec.h5 replays,
+        # whose stored TOML blob goes through the same path) parsing with a warning.
+        ffs_table = Dict{String,Any}("reform_eq_with_psilim" => false, "nn_low" => 1)
+        @test_logs (:warn, r"reform_eq_with_psilim") GeneralizedPerturbedEquilibrium._drop_deprecated_keys!(
+            ffs_table, GeneralizedPerturbedEquilibrium._DEPRECATED_FFS_KEYS, "ForceFreeStates")
+        @test !haskey(ffs_table, "reform_eq_with_psilim")
+        @test GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
+            (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl
+    end
+
     @testset "EFIT Method Consistency" begin
         # All three methods solve the same equilibrium — q-profiles should broadly agree.
         # Tolerance is 10% to allow for method-specific discretisation differences.
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index f23d10315..29955c2c7 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -411,8 +411,8 @@
         # Hit maxevals: error above tolerance
         @test_logs (:warn, r"maxevals_psi") KF.check_psi_quadrature_convergence(total, 0.5, ctrl, "fgar")
         # Nonzero atol_psi dominating a small torque: the silent-garbage scenario
-        ctrl.atol_psi = 1e-2
-        @test_logs (:warn, r"atol_psi") KF.check_psi_quadrature_convergence(1e-3 + 0.0im, 1e-3, ctrl, "fgar")
+        ctrl_atol = KF.KineticForcesControl(; atol_psi=1e-2)
+        @test_logs (:warn, r"atol_psi") KF.check_psi_quadrature_convergence(1e-3 + 0.0im, 1e-3, ctrl_atol, "fgar")
     end
 
     @testset "METHOD_REGISTRY" begin
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index 0536bf175..83b85a765 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -175,12 +175,12 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal()
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test isapprox(intr.qlim, equil.params.qmax; atol=1e-12)
-        @test isapprox(intr.psilim, equil.config.psihigh; atol=1e-12)
+        @test isapprox(intr.psilim, equil.params.psihigh_resolved; atol=1e-12)
 
         ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=max(equil.params.qmin + 0.1, equil.params.qmax - 0.5), set_psilim_via_dmlim=false)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test intr.qlim < equil.params.qmax + 1e-12
-        @test intr.psilim <= equil.config.psihigh
+        @test intr.psilim <= equil.params.psihigh_resolved
         q_at_psilim = equil.profiles.q_spline(intr.psilim)
         @test isapprox(q_at_psilim, intr.qlim; atol=1e-6)
         ctrl_dmlim = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=equil.params.qmax, set_psilim_via_dmlim=true)
diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml
index 6b9d537e6..bf3ceb3df 100644
--- a/test/test_data/regression_solovev_ideal_example/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
index 62e8569fb..d74b2884d 100644
--- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
index cca19bfff..38b1a2394 100644
--- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml
index d8349716b..60b8f67fa 100644
--- a/test/test_data/regression_solovev_kinetic_example/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
index 4527c7826..6e82349dc 100644
--- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
@@ -10,7 +10,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 16                               # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
 
diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
index ff4bc38ae..a3cd0efe6 100644
--- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
@@ -11,7 +11,7 @@ psilow = 1e-4                           # Lower limit of normalized flux coordin
 psihigh = 0.9995                        # Upper limit of normalized flux coordinate
 mpsi = 16                               # Number of radial grid points
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
 etol = 1e-7                             # Error tolerance for equilibrium solver
 force_termination = false               # Terminate after equilibrium setup (skip stability calculations)