diff --git a/Makefile b/Makefile
index b5dd8acdb..14a1c32e2 100644
--- a/Makefile
+++ b/Makefile
@@ -16,7 +16,7 @@ else
 endif
 
 GENERAL_REGISTRY_PACKAGES := CoordinateConventions FuseExchangeProtocol MillerExtendedHarmonic IMASdd
-FUSE_PACKAGES_MAKEFILE := ADAS BalanceOfPlantSurrogate BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields ThermalSystemModels # XSteam
+FUSE_PACKAGES_MAKEFILE := ADAS BalanceOfPlantSurrogate BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite FRESCO FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields ThermalSystemModels # XSteam
 FUSE_PACKAGES_MAKEFILE := $(sort $(FUSE_PACKAGES_MAKEFILE))
 FUSE_PACKAGES := $(shell echo '$(FUSE_PACKAGES_MAKEFILE)' | awk '{printf("[\"%s\"", $$1); for (i=2; i<=NF; i++) printf(", \"%s\"", $$i); print "]"}')
 DEV_PACKAGES_MAKEFILE := $(shell find ../*/.git/config -exec grep ProjectTorreyPines \{\} /dev/null \; | cut -d'/' -f 2)
@@ -113,6 +113,9 @@ QED:
 FiniteElementHermite:
 	$(call clone_pull_repo,$@)
 
+FRESCO:
+	$(call clone_pull_repo,$@)
+
 CHEASE:
 	$(call clone_pull_repo,$@)
 
diff --git a/Project.toml b/Project.toml
index 1e40d8ad9..fd936d4ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 EPEDNN = "e64856f0-3bb8-4376-b4b7-c03396503991"
 FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
+FRESCO = "f2db6eb9-2c54-4c18-912b-ce467796777a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FuseUtils = "c8e4ee58-0bec-4a6a-b552-9dd24d6374ac"
@@ -101,7 +102,7 @@ OrderedCollections = "1"
 Plots = "1"
 PolygonOps = "0.1"
 ProgressMeter = "1"
-QED = "1"
+QED = "1.4"
 QuadGK = "2"
 RABBIT = "1"
 SimulationParameters = "1.1"
diff --git a/src/FUSE.jl b/src/FUSE.jl
index 4350ff402..f13b143d2 100644
--- a/src/FUSE.jl
+++ b/src/FUSE.jl
@@ -74,7 +74,7 @@ include("actors.jl")
 
 include(joinpath("actors", "noop_actor.jl"))
 
-include(joinpath("actors", "equilibrium", "solovev_actor.jl"))
+include(joinpath("actors", "equilibrium", "fresco_actor.jl"))
 include(joinpath("actors", "equilibrium", "chease_actor.jl"))
 include(joinpath("actors", "equilibrium", "tequila_actor.jl"))
 include(joinpath("actors", "equilibrium", "equilibrium_actor.jl"))
@@ -97,6 +97,7 @@ include(joinpath("actors", "nuclear", "blanket_actor.jl"))
 include(joinpath("actors", "nuclear", "neutronics_actor.jl"))
 
 include(joinpath("actors", "current", "qed_actor.jl"))
+include(joinpath("actors", "current", "qed-coupled_actor.jl"))
 include(joinpath("actors", "current", "steadycurrent_actor.jl"))
 include(joinpath("actors", "current", "current_actor.jl"))
 
diff --git a/src/actors.jl b/src/actors.jl
index 668da39cc..97e3d5a71 100644
--- a/src/actors.jl
+++ b/src/actors.jl
@@ -50,24 +50,24 @@ end
 #  switch_get_from  #
 #= =============== =#
 """
-    switch_get_from(quantity::Symbol)
+    switch_get_from(quantity::Symbol; default::Union{Symbol,Missing}=missing)::Switch{Symbol}
 
-Switch to pick form which IDS `quantity` comes from
+Switch to pick the IDS that `quantity` comes from
 """
-function switch_get_from(quantity::Symbol)::Switch{Symbol}
+function switch_get_from(quantity::Symbol; default::Union{Symbol,Missing}=missing)::Switch{Symbol}
     txt = "Take $quantity from this IDS"
     if quantity == :ip
-        swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule], "-", txt)
+        swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule], "-", txt; default)
     elseif quantity == :vacuum_r0_b0
-        swch = Switch{Symbol}([:equilibrium, :pulse_schedule], "-", txt; default=:pulse_schedule)
+        swch = Switch{Symbol}([:equilibrium, :pulse_schedule], "-", txt; default)
     elseif quantity == :vloop
-        swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule, :controllers__ip], "-", txt)
+        swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule, :controllers__ip], "-", txt; default)
     elseif quantity == :βn
-        swch = Switch{Symbol}([:core_profiles, :equilibrium], "-", txt)
+        swch = Switch{Symbol}([:core_profiles, :equilibrium], "-", txt; default)
     elseif quantity == :ne_ped
-        swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt)
+        swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt; default)
     elseif quantity == :zeff_ped
-        swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt)
+        swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt; default)
     else
         error("`$quantity` not supported in switch_get_from()")
     end
diff --git a/src/actors/current/current_actor.jl b/src/actors/current/current_actor.jl
index 45ac1839d..e0c01211e 100644
--- a/src/actors/current/current_actor.jl
+++ b/src/actors/current/current_actor.jl
@@ -5,7 +5,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorCurrent{T<:Real} <: ParametersAc
     _parent::WeakRef = WeakRef(nothing)
     _name::Symbol = :not_set
     _time::Float64 = NaN
-    model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :none], "-", "Current actor to run"; default=:SteadyStateCurrent)
+    model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :QEDcoupled, :none], "-", "Current actor to run"; default=:SteadyStateCurrent)
     allow_floating_plasma_current::Entry{Bool} = Entry{Bool}("-", "Zero loop voltage if non-inductive fraction exceeds 100% of the target Ip"; default=true)
     #== data flow parameters ==#
     ip_from::Switch{Symbol} = switch_get_from(:ip)
@@ -15,7 +15,7 @@ end
 mutable struct ActorCurrent{D,P} <: CompoundAbstractActor{D,P}
     dd::IMAS.dd{D}
     par::FUSEparameters__ActorCurrent{P}
-    jt_actor::Union{ActorSteadyStateCurrent{D,P},ActorQED{D,P},ActorNoOperation{D,P}}
+    jt_actor::Union{ActorSteadyStateCurrent{D,P},ActorQED{D,P},ActorQEDcoupled{D,P},ActorNoOperation{D,P}}
 end
 
 """
@@ -43,6 +43,8 @@ function ActorCurrent(dd::IMAS.dd, par::FUSEparameters__ActorCurrent, act::Param
         jt_actor = ActorSteadyStateCurrent(dd, act.ActorSteadyStateCurrent; par.ip_from, par.allow_floating_plasma_current)
     elseif par.model == :QED
         jt_actor = ActorQED(dd, act.ActorQED; par.ip_from, par.vloop_from, par.allow_floating_plasma_current)
+    elseif par.model == :QEDcoupled
+        jt_actor = ActorQEDcoupled(dd, act.ActorQEDcoupled)
     end
     return ActorCurrent(dd, par, jt_actor)
 end
diff --git a/src/actors/current/qed-coupled_actor.jl b/src/actors/current/qed-coupled_actor.jl
new file mode 100644
index 000000000..73ecff907
--- /dev/null
+++ b/src/actors/current/qed-coupled_actor.jl
@@ -0,0 +1,230 @@
+import QED
+
+#= ======== =#
+#  ActorQEDcoupled  #
+#= ======== =#
+Base.@kwdef mutable struct FUSEparameters__ActorQEDcoupled{T<:Real} <: ParametersActor{T}
+    _parent::WeakRef = WeakRef(nothing)
+    _name::Symbol = :not_set
+    _time::Float64 = NaN
+    Δt::Entry{Float64} = Entry{Float64}("s", "Evolve for Δt")
+    Nt::Entry{Int} = Entry{Int}("-", "Number of time steps during evolution"; default=100)
+    #== display and debugging parameters ==#
+    debug::Entry{Bool} = Entry{Bool}("-", "Turn on QED debugging/plotting"; default=false)
+end
+
+mutable struct ActorQEDcoupled{D,P} <: SingleAbstractActor{D,P}
+    dd::IMAS.dd{D}
+    par::FUSEparameters__ActorQEDcoupled{P}
+    QO::Union{Nothing,QED.QED_state}
+    build::Union{Nothing,QED.QED_build}
+end
+
+"""
+    ActorQEDcoupled(dd::IMAS.dd, act::ParametersAllActors; kw...)
+
+Evolves the plasma and coil/vessel currents using the QED current diffusion solver.
+
+!!! note
+
+    This actor operates at "dd.global_time", any time advance must be done outside of the actor
+
+        IMAS.new_timeslice!(dd, dd.global_time + Δt)
+        dd.global_time += Δt
+        ActorQEDcoupled(dd, act)
+
+!!! note
+
+    Stores data in `dd.core_profiles.profiles_1d[].j_ohmic`
+"""
+function ActorQEDcoupled(dd::IMAS.dd, act::ParametersAllActors; kw...)
+    actor = ActorQEDcoupled(dd, act.ActorQEDcoupled; kw...)
+    step(actor)
+    finalize(actor)
+    return actor
+end
+
+function ActorQEDcoupled(dd::IMAS.dd, par::FUSEparameters__ActorQEDcoupled; kw...)
+    logging_actor_init(ActorQEDcoupled)
+    par = par(kw...)
+    return ActorQEDcoupled(dd, par, nothing, nothing)
+end
+
+function _step(actor::ActorQEDcoupled)
+    dd = actor.dd
+    par = actor.par
+
+    eqt = dd.equilibrium.time_slice[]
+    cp1d = dd.core_profiles.profiles_1d[]
+
+    # non_inductive contribution
+    B0 = eqt.global_quantities.vacuum_toroidal_field.b0
+    JBni = QED.FE(cp1d.grid.rho_tor_norm, cp1d.j_non_inductive .* B0)
+
+    # initialize QED
+    if actor.QO === nothing
+        actor.QO = qed_init_from_imas(eqt, cp1d; uniform_rho = 33)
+    else
+        actor.QO.JBni = JBni
+    end
+
+    actor.build = qed_build_from_imas(dd, dd.global_time - par.Δt)
+
+    # current diffusion
+    time0 = dd.global_time - 0.5 * par.Δt
+    actor.QO = QED.evolve(actor.QO, η_imas(dd.core_profiles.profiles_1d[time0]), actor.build, par.Δt, par.Nt; par.debug)
+
+    return actor
+end
+
+function _finalize(actor::ActorQEDcoupled)
+    dd = actor.dd
+
+    eqt = dd.equilibrium.time_slice[]
+    B0 = eqt.global_quantities.vacuum_toroidal_field.b0
+
+    cp1d = dd.core_profiles.profiles_1d[]
+    j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0
+
+    if ismissing(cp1d, :j_non_inductive)
+        cp1d.j_ohmic = j_total
+    else
+        cp1d.j_ohmic = j_total .- cp1d.j_non_inductive
+    end
+
+    # NEED TO POPULATE COIL CURRENTS BACK IN dd
+    build = actor.build
+    current_per_turn = build.Ic
+
+    active_coils = dd.pf_active.coil
+    passive_loops = dd.pf_passive.loop
+    for (k, coil) in enumerate(active_coils)
+        IMAS.@ddtime(coil.current.data = current_per_turn[k])
+    end
+    Nactive = length(active_coils)
+    for (k, loop) in enumerate(passive_loops)
+        IMAS.@ddtime(loop.current = current_per_turn[k + Nactive])
+    end
+
+    # MAYBE POPULATE RESISTANCE SO IT'S NOT RECALCULATED?
+
+    return actor
+end
+
+# utils
+"""
+    qed_init_from_imas(dd::IMAS.dd)
+
+Setup QED build from data in IMAS `dd`
+"""
+
+# function loopelement2coil(loop, element)
+#     outline = element.geometry.outline
+#     @assert length(outline.r) == 4 "For the time being passive structures must be composed of quadrilateral elements"
+#     passive_coil = VacuumFields.QuadCoil(outline.r, outline.z)
+#     element_turns = ismissing(element, :turns_with_sign) ? 1 : abs(element.turns_with_sign)
+#     loop_turns = sum(ismissing(elm, :turns_with_sign) ? 1 : abs(elm.turns_with_sign) for elm in loop.element)
+#     Ic = ismissing(loop, :current) ? 0.0 : loop.current / Nturns)
+#     VacuumFields.set_current!(passive_coil, Ic )
+#     passive_coil.resistance = VacuumFields.resistance(passive_coil, loop.resistivity)
+#     return passive_coil
+# end
+
+# elements(), turns(), QuadCoil(), and loop2multi() should all wind up in VacuumFields
+
+function qed_build_from_imas(dd::IMAS.dd{D}, time0::D) where {D <: Real}
+
+    coils = VacuumFields.MultiCoils(dd; load_pf_active=true, load_pf_passive=true);
+
+    # Coil-only quantities
+    Mcc = [VacuumFields.mutual(c1, c2) for c1 in coils, c2 in coils]
+    Ic = [VacuumFields.current(c) / VacuumFields.turns(c) for c in coils] #c urrent per turn
+    Rc = [VacuumFields.resistance(c) for c in coils];
+    Vc = zero(Ic);
+
+    eqt = dd.equilibrium.time_slice[time0]
+    cp1d = dd.core_profiles.profiles_1d[time0]
+    Ip = eqt.global_quantities.ip
+
+    # plasma-coil mutuals
+    image = VacuumFields.Image(eqt)
+    Mpc = [VacuumFields.mutual(image, coil, Ip) for coil in coils]
+    @warn "ActorQEDcoupled doesn't contain dMpc_dt yet"
+    dMpc_dt = zero(Mpc) # How Mpc changes in time (like shape)... to test later
+
+    # self inductance
+    It = IMAS.cumtrapz(cp1d.grid.area, cp1d.j_tor)
+    Wp = 0.5 * IMAS.trapz(cp1d.grid.psi, It)
+    Li = 2 * Wp / Ip^2 # internal
+    ψb = eqt.profiles_1d.psi[end]
+    ψc = sum(Mpc[k] * Ic[k] for k in eachindex(coils))
+    Le = (ψb - ψc) / Ip # external
+    Lp = Li + Le # total self
+
+    # plasma resistance
+    # BCL 9/25/24: from Pohm, which may be wrong
+    Pohm = dd.core_sources.source[:ohmic].profiles_1d[].electrons.power_inside[end]
+    Ini = IMAS.get_time_array(dd.core_profiles.global_quantities, :current_non_inductive, time0, :linear)
+    Iohm = Ip - Ini
+    Rp = Pohm / (Ip * Iohm)
+
+    # non-inductive voltage
+    Vni = Rp * Ini
+
+    # Waveforms
+    # These should come from pulse_schedule
+    #V_waveforms = fill(W0, length(coils));
+    #W = QED.Waveform{D}(t -> dd.pulse_schedule.voltage[1](t))
+
+    Nc = length(coils)
+    V_waveforms = Vector{QED.Waveform{D}}(undef, length(coils))
+    Vactive = Vwaveforms_from_pulse_schedule(dd, time0)
+    @assert length(Vactive) == length(dd.pf_active.coil)
+    Nactive = length(Vactive)
+    V_waveforms[1:Nactive] .= Vactive
+    V_waveforms[Nactive+1:end] .= fill(QED.Waveform{D}(t -> 0.0), Nc - Nactive)
+    return QED.QED_build(Ic, Vc, Rc, Mcc, Vni, Rp, Lp, Mpc, dMpc_dt, V_waveforms)
+end
+
+function Vwaveforms_from_pulse_schedule(dd::IMAS.dd{D}, t_start::D=IMAS.global_time(dd)) where {D<:Real}
+    Nc = length(dd.pf_active.coil)
+    time = dd.pulse_schedule.pf_active.time .- t_start
+    WFs = Vector{QED.Waveform{D}}(undef, Nc)
+    for (k, supply) in enumerate(dd.pulse_schedule.pf_active.supply)
+        if k in (1, 2)
+            WFs[k] = QED.Waveform(time, supply.voltage.reference)
+        elseif k == 3
+            Cu = VacuumFields.MultiCoil(dd.pf_active.coil[3])
+            Cl = VacuumFields.MultiCoil(dd.pf_active.coil[4])
+            Lu = VacuumFields.mutual(Cu, Cu)
+            Ll = VacuumFields.mutual(Cl, Cl)
+            M  = VacuumFields.mutual(Cu, Cl)
+            Lt = (Lu + 2M + Ll)
+
+            facu = (Lu + M) / Lt
+            WFs[3] = QED.Waveform(time , facu .* supply.voltage.reference)
+
+            facl = (Ll + M) / Lt
+            WFs[4] = QED.Waveform(time , facl .* supply.voltage.reference)
+        elseif k < 12
+            coil = dd.pf_active.coil[k+1]
+            WFs[k+1] = QED.Waveform(time, supply.voltage.reference)
+        else
+            Cu = VacuumFields.MultiCoil(dd.pf_active.coil[13])
+            Cl = VacuumFields.MultiCoil(dd.pf_active.coil[14])
+            Lu = VacuumFields.mutual(Cu, Cu)
+            Ll = VacuumFields.mutual(Cl, Cl)
+            M  = -VacuumFields.mutual(Cu, Cl) # oppositely connected
+            Lt = (Lu + 2M + Ll)
+
+            # positive voltage gives negative current in VS3U
+            facu = -(Lu + M) / Lt
+            WFs[13] = QED.Waveform(time , facu .* supply.voltage.reference)
+
+            # positive voltage gives positive current in VS3U
+            facl = (Ll + M) / Lt
+            WFs[14] = QED.Waveform(time , facl .* supply.voltage.reference)
+        end
+    end
+    return WFs
+end
diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl
index 3b04d9887..7861b84ec 100644
--- a/src/actors/current/qed_actor.jl
+++ b/src/actors/current/qed_actor.jl
@@ -95,6 +95,7 @@ function _step(actor::ActorQED)
         for time0 in range(t0 + δt / 2.0, t1 + δt / 2.0, No + 1)[1:end-1]
             if par.solve_for == :ip
                 Ip = IMAS.get_from(dd, Val{:ip}, par.ip_from; time0)
+                Ip_non_inductive = dd.core_profiles.global_quantities.current_non_inductive[time0]
                 Vedge = nothing
                 if par.allow_floating_plasma_current && abs(Ip) < abs(ip_non_inductive)
                     Ip = nothing
@@ -111,6 +112,7 @@ function _step(actor::ActorQED)
         # steady state solution
         if par.solve_for == :ip
             Ip = IMAS.get_from(dd, Val{:ip}, par.ip_from)
+            Ip_non_inductive = dd.core_profiles.global_quantities.current_non_inductive[]
             Vedge = nothing
             if par.allow_floating_plasma_current && abs(Ip) < abs(ip_non_inductive)
                 Ip = nothing
diff --git a/src/actors/equilibrium/chease_actor.jl b/src/actors/equilibrium/chease_actor.jl
index 3e6db8f22..46f2b4011 100644
--- a/src/actors/equilibrium/chease_actor.jl
+++ b/src/actors/equilibrium/chease_actor.jl
@@ -11,8 +11,6 @@ Base.@kwdef mutable struct FUSEparameters__ActorCHEASE{T<:Real} <: ParametersAct
     free_boundary::Entry{Bool} = Entry{Bool}("-", "Convert fixed boundary equilibrium to free boundary one"; default=true)
     clear_workdir::Entry{Bool} = Entry{Bool}("-", "Clean the temporary workdir for CHEASE"; default=true)
     rescale_eq_to_ip::Entry{Bool} = Entry{Bool}("-", "Scale equilibrium to match Ip"; default=true)
-    #== data flow parameters ==#
-    ip_from::Switch{Symbol} = switch_get_from(:ip)
 end
 
 mutable struct ActorCHEASE{D,P} <: SingleAbstractActor{D,P}
@@ -50,7 +48,7 @@ function _step(actor::ActorCHEASE)
 
     # initialize eqt from pulse_schedule and core_profiles
     eqt = dd.equilibrium.time_slice[]
-    eq1d = eqt.profiles_1d
+    eqt1d = eqt.profiles_1d
 
     # boundary
     pr = eqt.boundary.outline.r
@@ -69,9 +67,9 @@ function _step(actor::ActorCHEASE)
     ϵ = eqt.boundary.minor_radius / r_geo
 
     # pressure and j_tor
-    psin = eq1d.psi_norm
-    j_tor = [sign(j) == sign(Ip) ? j : 0.0 for j in eq1d.j_tor]
-    pressure = eq1d.pressure
+    psin = eqt1d.psi_norm
+    j_tor = [sign(j) == sign(Ip) ? j : 0.0 for j in eqt1d.j_tor]
+    pressure = eqt1d.pressure
     rho_pol = sqrt.(psin)
     pressure_sep = pressure[end]
 
@@ -162,8 +160,8 @@ function gEQDSK2IMAS(g::CHEASE.EFIT.GEQDSKFile, eq::IMAS.equilibrium)
     tc = MXHEquilibrium.transform_cocos(1, 11) # chease output is cocos 1 , dd is cocos 11
 
     eqt = eq.time_slice[]
-    eq1d = eqt.profiles_1d
-    eq2d = resize!(eqt.profiles_2d, 1)[1]
+    eqt1d = eqt.profiles_1d
+    eqt2d = resize!(eqt.profiles_2d, 1)[1]
 
     @ddtime(eq.vacuum_toroidal_field.b0 = g.bcentr)
     eq.vacuum_toroidal_field.r0 = g.rcentr
@@ -174,17 +172,17 @@ function gEQDSK2IMAS(g::CHEASE.EFIT.GEQDSKFile, eq::IMAS.equilibrium)
     eqt.global_quantities.magnetic_axis.z = g.zmaxis
     eqt.global_quantities.ip = g.current
 
-    eq1d.psi = g.psi .* tc["PSI"]
-    eq1d.q = g.qpsi
-    eq1d.pressure = g.pres
-    eq1d.dpressure_dpsi = g.pprime .* tc["PPRIME"]
-    eq1d.f = g.fpol .* tc["F"]
-    eq1d.f_df_dpsi = g.ffprim .* tc["F_FPRIME"]
-
-    eq2d.grid_type.index = 1
-    eq2d.grid.dim1 = g.r
-    eq2d.grid.dim2 = g.z
-    eq2d.psi = g.psirz .* tc["PSI"]
+    eqt1d.psi = g.psi .* tc["PSI"]
+    eqt1d.q = g.qpsi
+    eqt1d.pressure = g.pres
+    eqt1d.dpressure_dpsi = g.pprime .* tc["PPRIME"]
+    eqt1d.f = g.fpol .* tc["F"]
+    eqt1d.f_df_dpsi = g.ffprim .* tc["F_FPRIME"]
+
+    eqt2d.grid_type.index = 1
+    eqt2d.grid.dim1 = g.r
+    eqt2d.grid.dim2 = g.z
+    eqt2d.psi = g.psirz .* tc["PSI"]
 
     return nothing
 end
diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl
index 3d13be807..fbee1a702 100644
--- a/src/actors/equilibrium/equilibrium_actor.jl
+++ b/src/actors/equilibrium/equilibrium_actor.jl
@@ -6,12 +6,12 @@ Base.@kwdef mutable struct FUSEparameters__ActorEquilibrium{T<:Real} <: Paramete
     _name::Symbol = :not_set
     _time::Float64 = NaN
     #== actor parameters ==#
-    model::Switch{Symbol} = Switch{Symbol}([:Solovev, :CHEASE, :TEQUILA, :none], "-", "Equilibrium actor to run"; default=:TEQUILA)
+    model::Switch{Symbol} = Switch{Symbol}([:TEQUILA, :FRESCO, :CHEASE, :none], "-", "Equilibrium actor to run"; default=:TEQUILA)
     symmetrize::Entry{Bool} = Entry{Bool}("-", "Force equilibrium up-down symmetry with respect to magnetic axis"; default=false)
     #== data flow parameters ==#
-    j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles)
+    j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles, :pulse_schedule_pp_ffp], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles)
     ip_from::Switch{Symbol} = switch_get_from(:ip)
-    vacuum_r0_b0_from::Switch{Symbol} = switch_get_from(:vacuum_r0_b0)
+    vacuum_r0_b0_from::Switch{Symbol} = switch_get_from(:vacuum_r0_b0; default=:pulse_schedule)
     #== display and debugging parameters ==#
     do_plot::Entry{Bool} = act_common_parameters(; do_plot=false)
 end
@@ -20,7 +20,7 @@ mutable struct ActorEquilibrium{D,P} <: CompoundAbstractActor{D,P}
     dd::IMAS.dd{D}
     par::FUSEparameters__ActorEquilibrium{P}
     act::ParametersAllActors
-    eq_actor::Union{Nothing,ActorSolovev{D,P},ActorCHEASE{D,P},ActorTEQUILA{D,P},ActorNoOperation{D,P}}
+    eq_actor::Union{Nothing,ActorFRESCO{D,P},ActorCHEASE{D,P},ActorTEQUILA{D,P},ActorNoOperation{D,P}}
 end
 
 """
@@ -38,16 +38,16 @@ end
 function ActorEquilibrium(dd::IMAS.dd, par::FUSEparameters__ActorEquilibrium, act::ParametersAllActors; kw...)
     logging_actor_init(ActorEquilibrium)
     par = par(kw...)
-    if par.model == :Solovev
-        eq_actor = ActorSolovev(dd, act.ActorSolovev; par.ip_from)
+    if par.model == :FRESCO
+        eq_actor = ActorFRESCO(dd, act.ActorFRESCO)
     elseif par.model == :CHEASE
-        eq_actor = ActorCHEASE(dd, act.ActorCHEASE; par.ip_from)
+        eq_actor = ActorCHEASE(dd, act.ActorCHEASE)
     elseif par.model == :TEQUILA
         eq_actor = ActorTEQUILA(dd, act.ActorTEQUILA; par.ip_from)
     elseif par.model == :none
         eq_actor = ActorNoOperation(dd, act.ActorNoOperation)
     else
-        error("ActorEquilibrium: model = `$(par.model)` can only be one of [:Solovev, :CHEASE, :TEQUILA, :none]")
+        error("ActorEquilibrium: model = `$(par.model)` can only be one of [:TEQUILA, :FRESCO, :CHEASE, :none]")
     end
     return ActorEquilibrium(dd, par, act, eq_actor)
 end
@@ -74,6 +74,37 @@ function _step(actor::ActorEquilibrium)
         prepare(actor)
     end
 
+    # FRESCO needs p' and f'' to run
+    # We can estimate of p' and ff' based on equilibrium at previous time-slice
+    # if that's not available, we run TEQUILA first
+    if par.model == :FRESCO
+        eqt = dd.equilibrium.time_slice[]
+        idxeqt = IMAS.index(eqt)
+        if false
+            if false && idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip)
+                eqt1d = eqt.profiles_1d
+                eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d
+                psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm)
+                R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm8).(eqt1d.psi_norm)
+                one_R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm9).(eqt1d.psi_norm)
+                one_R2__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm1).(eqt1d.psi_norm)
+                b0 = eqt.global_quantities.vacuum_toroidal_field.b0
+                r0 = eqt.global_quantities.vacuum_toroidal_field.r0
+                dpressure_dpsi, f_df_dpsi, f = IMAS.calc_pprime_ffprim_f(psi__1, R__1, one_R__1, one_R2__1, r0, b0; eqt1d.pressure, eqt1d.j_tor)
+                eqt1d.dpressure_dpsi = dpressure_dpsi
+                eqt1d.f_df_dpsi = f_df_dpsi
+                eqt1d.f = f
+            else
+                # NOTE: using free_boundary=true will also update pf_active
+                #       which is then used as initial guess by FRESCO
+                ActorTEQUILA(dd, actor.act; free_boundary=true)
+                fw = IMAS.first_wall(dd.wall)
+                IMAS.flux_surfaces(eqt, fw.r, fw.z)
+            end
+        end
+
+    end
+
     # step selected equilibrium actor
     step(actor.eq_actor)
 
@@ -137,7 +168,7 @@ Prepare `dd.equilibrium` to run equilibrium actors
   - Clear equilibrium__time_slice
   - Set Ip, Bt from core_profiles, equilibrium, or pulse_schedule
   - Use position control from pulse_schedule
-  - Use j_tor,pressure from core_profiles or equilibrium
+  - Use j_tor,pressure from core_profiles (for self-consistent iterations) or equilibrium (to re-solve equilibrium with different solver)
 """
 function prepare(actor::ActorEquilibrium)
     dd = actor.dd
@@ -153,23 +184,30 @@ function prepare(actor::ActorEquilibrium)
         index = cp1d.grid.psi_norm .> 0.05
         psi0 = cp1d.grid.psi
         rho_tor_norm0 = cp1d.grid.rho_tor_norm
-        rho_pol_norm0 = vcat(-reverse(sqrt.(cp1d.grid.psi_norm[index])), sqrt.(cp1d.grid.psi_norm[index]))
+        rho_pol_norm_sqrt0 = vcat(-reverse(sqrt.(cp1d.grid.psi_norm[index])), sqrt.(cp1d.grid.psi_norm[index]))
         j_tor0 = vcat(reverse(cp1d.j_tor[index]), cp1d.j_tor[index])
         pressure0 = vcat(reverse(cp1d.pressure[index]), cp1d.pressure[index])
-        j_itp = IMAS.interp1d(rho_pol_norm0, j_tor0, :cubic)
-        p_itp = IMAS.interp1d(rho_pol_norm0, pressure0, :cubic)
+        j_itp = IMAS.interp1d(rho_pol_norm_sqrt0, j_tor0, :cubic)
+        p_itp = IMAS.interp1d(rho_pol_norm_sqrt0, pressure0, :cubic)
     elseif par.j_p_from == :equilibrium
         @assert !isempty(dd.equilibrium.time)
         eqt1d = dd.equilibrium.time_slice[].profiles_1d
         psi0 = eqt1d.psi
         rho_tor_norm0 = eqt1d.rho_tor_norm
-        rho_pol_norm0 = sqrt.(eqt1d.psi_norm)
+        rho_pol_norm_sqrt0 = sqrt.(eqt1d.psi_norm)
         j_tor0 = eqt1d.j_tor
         pressure0 = eqt1d.pressure
-        j_itp = IMAS.interp1d(rho_pol_norm0, j_tor0, :cubic)
-        p_itp = IMAS.interp1d(rho_pol_norm0, pressure0, :cubic)
+        j_itp = IMAS.interp1d(rho_pol_norm_sqrt0, j_tor0, :cubic)
+        p_itp = IMAS.interp1d(rho_pol_norm_sqrt0, pressure0, :cubic)
+    elseif par.j_p_from == :pulse_schedule_pp_ffp
+        @assert !isempty(dd.equilibrium.time)
+        eqt1d = dd.equilibrium.time_slice[].profiles_1d
+        psi0 = eqt1d.psi
+        rho_tor_norm0 = eqt1d.rho_tor_norm
+        rho_pol_norm_sqrt0 = sqrt.(eqt1d.psi_norm)
+        pend = eqt1d.pressure[end]
     else
-        @assert par.j_p_from in (:core_profiles, :equilibrium)
+        @assert par.j_p_from in (:core_profiles, :equilibrium, :pulse_schedule_pp_ffp)
     end
 
     # get ip and b0 before wiping eqt in case ip_from=:equilibrium
@@ -179,7 +217,7 @@ function prepare(actor::ActorEquilibrium)
     # add/clear time-slice
     eqt = resize!(dd.equilibrium.time_slice)
     resize!(eqt.profiles_2d, 1)
-    eq1d = dd.equilibrium.time_slice[].profiles_1d
+    eqt1d = dd.equilibrium.time_slice[].profiles_1d
 
     # scalar quantities
     eqt.global_quantities.ip = ip
@@ -230,12 +268,28 @@ function prepare(actor::ActorEquilibrium)
         end
     end
 
-    # set j_tor and pressure, forcing zero derivative on axis
-    eq1d = dd.equilibrium.time_slice[].profiles_1d
-    eq1d.psi = psi0
-    eq1d.rho_tor_norm = rho_tor_norm0
-    eq1d.j_tor = j_itp.(sqrt.(eq1d.psi_norm))
-    eq1d.pressure = p_itp.(sqrt.(eq1d.psi_norm))
+    eqt1d = dd.equilibrium.time_slice[].profiles_1d
+    eqt1d.psi = psi0
+    eqt1d.rho_tor_norm = rho_tor_norm0
+    if par.j_p_from === :pulse_schedule_pp_ffp
+        time0 = dd.global_time
+        profcon = dd.pulse_schedule.profiles_control
+        psin_ps = profcon.psi_norm
+
+        pprime  = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.dpressure_dpsi, :reference, [time0])[:,1], :cubic)
+        ffprime = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.f_df_dpsi, :reference, [time0])[:,1], :cubic)
+        #fpol = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.f, :reference, [time0])[:,1], :cubic)
+
+        psin_eq = eqt1d.psi_norm
+        eqt1d.dpressure_dpsi = pprime.(psin_eq)
+        eqt1d.pressure = IMAS.cumtrapz(eqt1d.psi, eqt1d.dpressure_dpsi)
+        eqt1d.pressure .+= pend .- eqt1d.pressure[end]
+        eqt1d.f_df_dpsi = ffprime.(psin_eq)
+    else
+        # set j_tor and pressure, forcing zero derivative on axis
+        eqt1d.j_tor = j_itp.(sqrt.(eqt1d.psi_norm))
+        eqt1d.pressure = p_itp.(sqrt.(eqt1d.psi_norm))
+    end
 
     return dd
 end
@@ -297,4 +351,4 @@ function IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice)
         (eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z), # Magnetic Axis (raxis,zaxis)
         Int(sign(eqt.profiles_1d.f[end]) * sign(eqt.global_quantities.ip)) # sign(dot(J,B))
     )
-end
\ No newline at end of file
+end
diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl
new file mode 100644
index 000000000..05bee0b71
--- /dev/null
+++ b/src/actors/equilibrium/fresco_actor.jl
@@ -0,0 +1,128 @@
+import FRESCO
+
+#= =========== =#
+#  ActorFRESCO  #
+#= =========== =#
+Base.@kwdef mutable struct FUSEparameters__ActorFRESCO{T<:Real} <: ParametersActor{T}
+    _parent::WeakRef = WeakRef(nothing)
+    _name::Symbol = :not_set
+    _time::Float64 = NaN
+    #== actor parameters ==#
+    control::Switch{Symbol} = Switch{Symbol}([:vertical, :shape], "-", ""; default=:shape)
+    number_of_iterations::Entry{Tuple{Int,Int}} = Entry{Tuple{Int,Int}}("-", "Number of outer and inner iterations"; default=(100, 3))
+    relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.5)
+    tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4)
+    #== data flow parameters ==#
+    fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal)
+    #== display and debugging parameters ==#
+    do_plot::Entry{Bool} = act_common_parameters(; do_plot=false)
+    debug::Entry{Bool} = Entry{Bool}("-", "Print debug information withing FRESCO solve"; default=true)
+    #== IMAS psi grid settings ==#
+    nR::Entry{Int} = Entry{Int}("-", "Grid resolution along R"; default=129)
+    nZ::Entry{Int} = Entry{Int}("-", "Grid resolution along Z"; default=129)
+end
+
+mutable struct ActorFRESCO{D,P} <: SingleAbstractActor{D,P}
+    dd::IMAS.dd{D}
+    par::FUSEparameters__ActorFRESCO{P}
+    canvas::Union{Nothing,FRESCO.Canvas}
+    profile::Union{Nothing,FRESCO.CurrentProfile}
+end
+
+"""
+    ActorFRESCO(dd::IMAS.dd, act::ParametersAllActors; kw...)
+
+Runs the Fixed boundary equilibrium solver FRESCO
+"""
+function ActorFRESCO(dd::IMAS.dd, act::ParametersAllActors; kw...)
+    actor = ActorFRESCO(dd, act.ActorFRESCO; kw...)
+    step(actor)
+    finalize(actor)
+    return actor
+end
+
+function ActorFRESCO(dd::IMAS.dd{D}, par::FUSEparameters__ActorFRESCO{P}; kw...) where {D<:Real,P<:Real}
+    logging_actor_init(ActorFRESCO)
+    par = par(kw...)
+    return ActorFRESCO(dd, par, nothing, nothing)
+end
+
+"""
+    _step(actor::ActorFRESCO)
+
+Runs FRESCO on the r_z boundary, equilibrium pressure and equilibrium j_tor
+"""
+function _step(actor::ActorFRESCO)
+    dd = actor.dd
+    par = actor.par
+    eqt = dd.equilibrium.time_slice[]
+    eqt1d = eqt.profiles_1d
+
+    # Ip_target = eqt.global_quantities.ip
+    # if par.fixed_grid === :poloidal
+    #     rho = sqrt.(eqt1d.psi_norm)
+    #     rho[1] = 0.0
+    #     P = (FRESCO.FE(rho, eqt1d.pressure), :poloidal)
+    #     # don't allow current to change sign
+    #     Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :poloidal)
+    #     Pbnd = eqt1d.pressure[end]
+    # elseif par.fixed_grid === :toroidal
+    #     rho = eqt1d.rho_tor_norm
+    #     P = (FRESCO.FE(rho, eqt1d.pressure), :toroidal)
+    #     # don't allow current to change sign
+    #     Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :toroidal)
+    #     Pbnd = eqt1d.pressure[end]
+    # end
+
+    gpp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.dpressure_dpsi, :cubic)
+    gffp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.f_df_dpsi, :cubic)
+    actor.profile = FRESCO.PprimeFFprime(x -> gpp(x), x -> gffp(x))
+
+    actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ; load_pf_active=true, load_pf_passive=true)
+
+    FRESCO.solve!(actor.canvas, actor.profile, par.number_of_iterations...; par.relax, par.debug, par.control, par.tolerance)
+
+    # using Plots
+    # p1 = Plots.heatmap(Rs, Zs, psi', aspect_ratio=:equal, xrange=(0,12.5), yrange=(-9,9), size=(400,500))
+    # Plots.plot!(p1, coils, label=nothing)
+    # Plots.contour!(p1, Rs, Zs, psi', levels=[C.Ψbnd], color=:white)
+    # display(p1)
+
+    return actor
+end
+
+# finalize by converting FRESCO canvas to dd.equilibrium
+function _finalize(actor::ActorFRESCO)
+    canvas = actor.canvas
+    profile = actor.profile
+    dd = actor.dd
+    eq = dd.equilibrium
+    eqt = eq.time_slice[]
+    eqt1d = eqt.profiles_1d
+    eq2d = resize!(eqt.profiles_2d, 1)[1]
+
+    Raxis, Zaxis, _ = FRESCO.find_axis(canvas)
+    eqt.global_quantities.magnetic_axis.r = Raxis
+    eqt.global_quantities.magnetic_axis.z = Zaxis
+    eqt.global_quantities.psi_boundary = canvas.Ψbnd
+    eqt.global_quantities.psi_axis = canvas.Ψaxis
+    eqt1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eqt1d.psi))
+    # p' doesn't change
+    eqt1d.f_df_dpsi .*= profile.ffp_scale
+
+    fend = dd.equilibrium.vacuum_toroidal_field.b0[] * dd.equilibrium.vacuum_toroidal_field.r0[]
+    f2 = 2 * IMAS.cumtrapz(eqt1d.psi, eqt1d.f_df_dpsi)
+    f2 .= f2 .- f2[end] .+ fend ^ 2
+    eqt1d.f = sign(fend) .* sqrt.(f2)
+
+    pend = eqt1d.pressure[end]
+    eqt1d.pressure = IMAS.cumtrapz(eqt1d.psi, eqt1d.dpressure_dpsi)
+    eqt1d.pressure .+= pend .- eqt1d.pressure[end]
+
+    eq2d.grid_type.index = 1
+    eq2d.grid.dim1 = canvas.Rs
+    eq2d.grid.dim2 = canvas.Zs
+    eq2d.psi = canvas.Ψ
+
+    return actor
+end
diff --git a/src/actors/equilibrium/solovev_actor.jl b/src/actors/equilibrium/solovev_actor.jl
deleted file mode 100644
index 8ce4d0d58..000000000
--- a/src/actors/equilibrium/solovev_actor.jl
+++ /dev/null
@@ -1,199 +0,0 @@
-import MXHEquilibrium
-import Optim
-
-#= ============ =#
-#  ActorSolovev  #
-#= ============ =#
-Base.@kwdef mutable struct FUSEparameters__ActorSolovev{T<:Real} <: ParametersActor{T}
-    _parent::WeakRef = WeakRef(nothing)
-    _name::Symbol = :not_set
-    _time::Float64 = NaN
-    #== actor parameters ==#
-    ngrid::Entry{Int} = Entry{Int}("-", "Grid size (for R, Z follows proportionally to plasma elongation)"; default=129)
-    qstar::Entry{T} = Entry{T}("-", "Initial guess of kink safety factor"; default=1.5)
-    alpha::Entry{T} = Entry{T}("-", "Initial guess of constant relating to pressure"; default=0.0)
-    #== data flow parameters ==#
-    ip_from::Switch{Symbol} = switch_get_from(:ip)
-    #== display and debugging parameters ==#
-    verbose::Entry{Bool} = act_common_parameters(; verbose=false)
-end
-
-mutable struct ActorSolovev{D,P} <: SingleAbstractActor{D,P}
-    dd::IMAS.dd{D}
-    par::FUSEparameters__ActorSolovev{P}
-    S::Union{Nothing,MXHEquilibrium.SolovevEquilibrium}
-end
-
-"""
-    ActorSolovev(dd::IMAS.dd, act::ParametersAllActors; kw...)
-
-Solovev equilibrium actor, based on:
-“One size fits all” analytic solutions to the Grad–Shafranov equation
-Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818
-"""
-function ActorSolovev(dd::IMAS.dd, act::ParametersAllActors; kw...)
-    actor = ActorSolovev(dd, act.ActorSolovev; kw...)
-    step(actor)
-    finalize(actor)
-    return actor
-end
-
-function ActorSolovev(dd::IMAS.dd, par::FUSEparameters__ActorSolovev; kw...)
-    logging_actor_init(ActorSolovev)
-    par = par(kw...)
-    return ActorSolovev(dd, par, nothing)
-end
-
-"""
-    _step(actor::ActorSolovev)
-
-Non-linear optimization to obtain a target `ip` and `pressure_core`
-"""
-function _step(actor::ActorSolovev)
-    dd = actor.dd
-    par = actor.par
-
-    # initialize eqt from pulse_schedule and core_profiles
-    eq = dd.equilibrium
-    eqt = eq.time_slice[]
-
-    # magnetic field
-    B0 = @ddtime eq.vacuum_toroidal_field.b0
-    target_ip = eqt.global_quantities.ip
-    target_pressure_core = eqt.profiles_1d.pressure[1]
-
-    # plasma shape as MXH
-    mxh = IMAS.MXH(eqt.boundary.outline.r, eqt.boundary.outline.z, 4)
-    Z0off = mxh.Z0 # Solovev has a bug for Z!=0.0
-    mxh.Z0 -= Z0off
-    plasma_shape = MXHEquilibrium.MillerExtendedHarmonicShape(mxh.R0, mxh.Z0, mxh.ϵ, mxh.κ, mxh.c0, mxh.c, mxh.s)
-    # check number of x_points to infer symmetry
-    if mod(length(eqt.boundary.x_point), 2) == 0
-        symmetric = true
-    else
-        symmetric = false
-    end
-
-    # add x_point info
-    if length(eqt.boundary.x_point) > 0
-        x_point = (eqt.boundary.x_point[1].r, -abs(eqt.boundary.x_point[1].z) - Z0off)
-    else
-        x_point = nothing
-    end
-
-    # first run of Solovev
-    # display(plot(mxh))
-    # @show abs(B0)
-    # @show plasma_shape
-    # @show par.alpha
-    # @show par.qstar
-    # @show x_point
-    # @show symmetric
-    S0 = MXHEquilibrium.solovev(abs(B0), plasma_shape, par.alpha, par.qstar; B0_dir=Int64(sign(B0)), Ip_dir=1, x_point, symmetric)
-
-    # optimize `alpha` and `qstar` to get target_ip and target_pressure_core
-    function cost(x)
-        S = MXHEquilibrium.solovev(abs(B0), plasma_shape, x[1], x[2]; B0_dir=Int64(sign(B0)), Ip_dir=1, symmetric, x_point)
-        psimag, psibry = MXHEquilibrium.psi_limits(S)
-        pressure_cost = (MXHEquilibrium.pressure(S, psimag) - target_pressure_core) / target_pressure_core
-        ip_cost = (MXHEquilibrium.plasma_current(S) - target_ip) / target_ip
-        return sqrt(pressure_cost^2 + 100.0 * ip_cost^2)
-    end
-    res = Optim.optimize(cost, [S0.alpha, S0.qstar], Optim.NelderMead(), Optim.Options(; g_tol=1E-3))
-    if par.verbose
-        println(res)
-    end
-
-    # final Solovev solution
-    actor.S = MXHEquilibrium.solovev(abs(B0), plasma_shape, res.minimizer[1], res.minimizer[2]; B0_dir=Int64(sign(B0)), Ip_dir=1, symmetric, x_point)
-
-    return actor
-end
-
-"""
-    _finalize(actor::ActorSolovev)
-
-Store ActorSolovev data in IMAS.equilibrium format
-"""
-function _finalize(actor::ActorSolovev)
-    dd = actor.dd
-    par = actor.par
-    mxh_eq = actor.S
-
-    eqt = dd.equilibrium.time_slice[]
-
-    target_ip = eqt.global_quantities.ip
-    target_psi_norm = getproperty(eqt.profiles_1d, :psi_norm, missing)
-    target_pressure = getproperty(eqt.profiles_1d, :pressure, missing)
-    target_j_tor = getproperty(eqt.profiles_1d, :j_tor, missing)
-
-    MXHEquilibrium_to_dd!(dd.equilibrium, mxh_eq, par.ngrid; cocos_in=3)
-
-    # force total plasma current to target_ip to avoid drifting after multiple calls of SolovevActor
-    eqt2d = findfirst(:rectangular, eqt.profiles_2d)
-    eqt2d.psi = (eqt2d.psi .- eqt.profiles_1d.psi[end]) .* (target_ip / eqt.global_quantities.ip) .+ eqt.profiles_1d.psi[end]
-    eqt.profiles_1d.psi = (eqt.profiles_1d.psi .- eqt.profiles_1d.psi[end]) .* (target_ip / eqt.global_quantities.ip) .+ eqt.profiles_1d.psi[end]
-    # match entry target_pressure and target_j_tor as if Solovev could do this
-    if !ismissing(target_pressure)
-        eqt.profiles_1d.pressure = IMAS.interp1d(target_psi_norm, target_pressure).(eqt.profiles_1d.psi_norm)
-    end
-    if !ismissing(target_j_tor)
-        eqt.profiles_1d.j_tor = IMAS.interp1d(target_psi_norm, target_j_tor, :cubic).(eqt.profiles_1d.psi_norm)
-    end
-    IMAS.p_jtor_2_pprime_ffprim_f!(eqt.profiles_1d, mxh_eq.S.R0, mxh_eq.B0)
-
-    # record optimized values of qstar and alpha in `act` for subsequent calls to the same actor
-    actor.par.qstar = actor.S.qstar
-    actor.par.alpha = actor.S.alpha
-
-    return actor
-end
-
-function MXHEquilibrium_to_dd!(eq::IMAS.equilibrium, mxh_eq::MXHEquilibrium.AbstractEquilibrium, ngrid::Int; cocos_in::Int=11)
-    tc = MXHEquilibrium.transform_cocos(cocos_in, 11)
-
-    rlims = MXHEquilibrium.limits(mxh_eq.S; pad=0.3)[1]
-    zlims = MXHEquilibrium.limits(mxh_eq.S; pad=0.3)[2]
-
-    eqt = eq.time_slice[]
-    Ip = eqt.global_quantities.ip
-    sign_Ip = sign(Ip)
-    sign_Bt = sign(eqt.profiles_1d.f[end])
-
-    Z0 = eqt.boundary.geometric_axis.z
-    flip_z = 1.0
-    if mod(length(eqt.boundary.x_point), 2) == 1 && eqt.boundary.x_point[1].z > Z0
-        flip_z = -1.0
-    end
-
-    eq.vacuum_toroidal_field.r0 = mxh_eq.S.R0
-    @ddtime eq.vacuum_toroidal_field.b0 = mxh_eq.B0 * sign_Bt
-
-    t0 = eqt.time
-    empty!(eqt)
-    eqt.time = t0
-
-    eqt.global_quantities.ip = Ip
-    eqt.boundary.geometric_axis.r = mxh_eq.S.R0
-    eqt.boundary.geometric_axis.z = Z0
-    orig_psi = collect(range(MXHEquilibrium.psi_limits(mxh_eq)..., ngrid))
-    eqt.profiles_1d.psi = orig_psi * (tc["PSI"] * sign_Ip)
-
-    eqt.profiles_1d.pressure = MXHEquilibrium.pressure.(mxh_eq, orig_psi)
-    eqt.profiles_1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(mxh_eq, orig_psi) ./ (tc["PSI"] * sign_Ip)
-
-    eqt.profiles_1d.f = MXHEquilibrium.poloidal_current.(mxh_eq, orig_psi) .* (tc["F"] * sign_Bt)
-    eqt.profiles_1d.f_df_dpsi =
-        MXHEquilibrium.poloidal_current.(mxh_eq, orig_psi) .* MXHEquilibrium.poloidal_current_gradient.(mxh_eq, orig_psi) .* (tc["F_FPRIME"] * sign_Bt * sign_Ip)
-
-    eqt2d = resize!(eqt.profiles_2d, 1)[1]
-    eqt2d.grid_type.index = 1
-    eqt2d.grid.dim1 = range(rlims[1], rlims[2], ngrid)
-    eqt2d.grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, Int(ceil(ngrid * mxh_eq.S.κ)))
-    eqt2d.psi = [mxh_eq(rr, flip_z * (zz - Z0)) * (tc["PSI"] * sign_Ip) for rr in eqt2d.grid.dim1, zz in eqt2d.grid.dim2]
-
-    fw = IMAS.first_wall(IMAS.top_dd(eqt).wall)
-    IMAS.flux_surfaces(eqt, fw.r, fw.z)
-
-    return
-end
\ No newline at end of file
diff --git a/src/actors/equilibrium/tequila_actor.jl b/src/actors/equilibrium/tequila_actor.jl
index 3a9c8cce2..1e021de53 100644
--- a/src/actors/equilibrium/tequila_actor.jl
+++ b/src/actors/equilibrium/tequila_actor.jl
@@ -16,7 +16,6 @@ Base.@kwdef mutable struct FUSEparameters__ActorTEQUILA{T<:Real} <: ParametersAc
     relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.25, check=x -> @assert 0.0 <= x <= 1.0 "must be: 0.0 <= relax <= 1.0")
     tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4)
     #== data flow parameters ==#
-    ip_from::Switch{Symbol} = switch_get_from(:ip)
     fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal)
     #== display and debugging parameters ==#
     do_plot::Entry{Bool} = act_common_parameters(; do_plot=false)
@@ -63,27 +62,26 @@ function _step(actor::ActorTEQUILA)
     D = eltype(dd)
     par = actor.par
     eqt = dd.equilibrium.time_slice[]
-    eq1d = eqt.profiles_1d
+    eqt1d = eqt.profiles_1d
 
     # BCL 5/30/23: ψbound should be set time dependently, related to the flux swing of the OH coils
-    #              For now setting to zero as initial eq1d.psi profile from prepare() can be nonsense
+    #              For now setting to zero as initial eqt1d.psi profile from prepare() can be nonsense
     actor.ψbound = D(0.0)
 
     Ip_target = eqt.global_quantities.ip
-
     if par.fixed_grid === :poloidal
-        rho = sqrt.(eq1d.psi_norm)
+        rho = sqrt.(eqt1d.psi_norm)
         rho[1] = D(0.0)
-        P = (TEQUILA.FE(rho, eq1d.pressure), :poloidal)
+        P = (TEQUILA.FE(rho, eqt1d.pressure), :poloidal)
         # don't allow current to change sign
-        Jt = (TEQUILA.FE(rho, D[sign(j) == sign(Ip_target) ? j : D(0.0) for j in eq1d.j_tor]), :poloidal)
-        Pbnd = eq1d.pressure[end]
+        Jt = (TEQUILA.FE(rho, D[sign(j) == sign(Ip_target) ? j : D(0.0) for j in eqt1d.j_tor]), :poloidal)
+        Pbnd = eqt1d.pressure[end]
     elseif par.fixed_grid === :toroidal
-        rho = eq1d.rho_tor_norm
-        P = (TEQUILA.FE(rho, eq1d.pressure), :toroidal)
+        rho = eqt1d.rho_tor_norm
+        P = (TEQUILA.FE(rho, eqt1d.pressure), :toroidal)
         # don't allow current to change sign
-        Jt = (TEQUILA.FE(rho, D[sign(j) == sign(Ip_target) ? j : D(0.0) for j in eq1d.j_tor]), :toroidal)
-        Pbnd = eq1d.pressure[end]
+        Jt = (TEQUILA.FE(rho, D[sign(j) == sign(Ip_target) ? j : D(0.0) for j in eqt1d.j_tor]), :toroidal)
+        Pbnd = eqt1d.pressure[end]
     end
 
     Fbnd = eqt.global_quantities.vacuum_toroidal_field.b0 * eqt.global_quantities.vacuum_toroidal_field.r0
@@ -125,8 +123,8 @@ function _step(actor::ActorTEQUILA)
     catch e
         plot(eqt.boundary.outline.r, eqt.boundary.outline.z; marker=:dot, aspect_ratio=:equal)
         display(plot!(IMAS.MXH(actor.shot.surfaces[:, end])))
-        display(plot(rho, eq1d.pressure; marker=:dot, xlabel="ρ", title="Pressure [Pa]"))
-        display(plot(rho, eq1d.j_tor; marker=:dot, xlabel="ρ", title="Jtor [A]"))
+        display(plot(rho, eqt1d.pressure; marker=:dot, xlabel="ρ", title="Pressure [Pa]"))
+        display(plot(rho, eqt1d.j_tor; marker=:dot, xlabel="ρ", title="Jtor [A]"))
         rethrow(e)
     end
 
@@ -149,7 +147,7 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd{D}, par::FUSEparameters__A
     free_boundary = par.free_boundary
     eq = dd.equilibrium
     eqt = eq.time_slice[]
-    eq1d = eqt.profiles_1d
+    eqt1d = eqt.profiles_1d
 
     R0 = shot.R0fe(1.0)
     Z0 = shot.Z0fe(1.0)
@@ -174,12 +172,12 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd{D}, par::FUSEparameters__A
     psit = shot.C[2:2:end, 1]
     psia = psit[1]
     psib = psit[end]
-    eq1d.psi = range(psia, psib, nψ_grid)
-    rhoi = TEQUILA.ρ.(Ref(shot), eq1d.psi)
-    eq1d.pressure = MXHEquilibrium.pressure.(Ref(shot), eq1d.psi)
-    eq1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(Ref(shot), eq1d.psi)
-    eq1d.f = shot.F.(rhoi)
-    eq1d.f_df_dpsi = TEQUILA.Fpol_dFpol_dψ.(Ref(shot), rhoi; shot.invR, shot.invR2)
+    eqt1d.psi = range(psia, psib, nψ_grid)
+    rhoi = TEQUILA.ρ.(Ref(shot), eqt1d.psi)
+    eqt1d.pressure = MXHEquilibrium.pressure.(Ref(shot), eqt1d.psi)
+    eqt1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(Ref(shot), eqt1d.psi)
+    eqt1d.f = shot.F.(rhoi)
+    eqt1d.f_df_dpsi = TEQUILA.Fpol_dFpol_dψ.(Ref(shot), rhoi; shot.invR, shot.invR2)
 
     resize!(eqt.profiles_2d, 2)
 
@@ -253,7 +251,7 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd{D}, par::FUSEparameters__A
                 eq2d.psi[i, j] = shot(r, z; extrapolate=true) + ψbound
             end
         end
-        eq1d.psi .+= ψbound
+        eqt1d.psi .+= ψbound
     end
 
 end
diff --git a/src/cases/_toksys.jl b/src/cases/_toksys.jl
index 866810e45..4e21bd9dd 100644
--- a/src/cases/_toksys.jl
+++ b/src/cases/_toksys.jl
@@ -61,11 +61,8 @@ function from_TokSys(dd::IMAS.dd, mat::Dict; what...)
         dd.equilibrium.time = eq_time
         dd.equilibrium.vacuum_toroidal_field.r0 = r0
         dd.equilibrium.vacuum_toroidal_field.b0 = b0
-        for k in eachindex(dd.equilibrium.time_slice)
-            IMAS.retime!(dd.equilibrium.time_slice[k], eq_time[k])
-        end
-        resize!(dd.equilibrium.time_slice, length(eq_time))
-        for (k, eqt) in enumerate(dd.equilibrium.time_slice)
+        for (k, time) in enumerate(collect(eq_time))
+            eqt = resize!(dd.equilibrium.time_slice, time)
             eq1d = eqt.profiles_1d
             eq2d = resize!(eqt.profiles_2d, 1)[1]
             eq2d.grid_type.index = 1