diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 1f1d63e17..2d6f37946 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -105,6 +105,7 @@ namespace Opm {
             using RateVector = GetPropType<TypeTag, Properties::RateVector>;
             using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
             using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
+            using GasLiftSingleWell = typename WellInterface<TypeTag>::GasLiftSingleWell;
             using GLiftOptWells = typename BlackoilWellModelGeneric::GLiftOptWells;
             using GLiftProdWells = typename BlackoilWellModelGeneric::GLiftProdWells;
             using GLiftWellStateMap =
@@ -398,6 +399,13 @@ namespace Opm {
                 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
                 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map);
 
+            // cannot be const since it accesses the non-const WellState
+            void gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
+                DeferredLogger& deferred_logger,
+                GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
+                GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
+                GLiftSyncGroups& groups_to_sync);
+
             void extractLegacyCellPvtRegionIndex_();
 
             void extractLegacyDepth_();
diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 66d8f878a..3b1b7168b 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -956,11 +956,9 @@ namespace Opm {
                 for (const auto& well : well_container_) {
                     // NOTE: Only the wells in "group_info" needs to be optimized
                     if (group_info.hasWell(well->name())) {
-                        well->gasLiftOptimizationStage1(
-                            this->wellState(), this->groupState(), ebosSimulator_,
-                            deferred_logger,
-                            prod_wells, glift_wells, state_map,
-                            group_info, groups_to_sync, this->glift_debug
+                        gasLiftOptimizationStage1SingleWell(
+                            well.get(), deferred_logger, prod_wells, glift_wells,
+                            group_info, state_map, groups_to_sync
                         );
                     }
                 }
@@ -1028,6 +1026,35 @@ namespace Opm {
         }
     }
 
+    // NOTE: this method cannot be const since it passes this->wellState()
+    //   (see below) to the GasLiftSingleWell constructor which accepts WellState
+    //   as a non-const reference..
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
+        DeferredLogger& deferred_logger,
+        GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
+        GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
+        GLiftSyncGroups& sync_groups)
+    {
+        const auto& summary_state = ebosSimulator_.vanguard().summaryState();
+        std::unique_ptr<GasLiftSingleWell> glift
+            = std::make_unique<GasLiftSingleWell>(
+                *well, ebosSimulator_, summary_state,
+                deferred_logger, this->wellState(), this->groupState(),
+                group_info, sync_groups, this->glift_debug);
+        auto state = glift->runOptimize(
+            ebosSimulator_.model().newtonMethod().numIterations());
+        if (state) {
+            state_map.insert({well->name(), std::move(state)});
+            glift_wells.insert({well->name(), std::move(glift)});
+            return;
+        }
+        prod_wells.insert({well->name(), well});
+    }
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
diff --git a/opm/simulators/wells/GasLiftGroupInfo.cpp b/opm/simulators/wells/GasLiftGroupInfo.cpp
index dde86eb70..5204e7ce1 100644
--- a/opm/simulators/wells/GasLiftGroupInfo.cpp
+++ b/opm/simulators/wells/GasLiftGroupInfo.cpp
@@ -365,6 +365,54 @@ checkNewtonIterationIdxOk_(const std::string &well_name)
     }
 }
 
+void
+GasLiftGroupInfo::
+debugDisplayWellContribution_(
+    const std::string& gr_name, const std::string& well_name,
+    double eff_factor,
+    double well_oil_rate, double well_gas_rate, double well_water_rate,
+    double well_alq,
+    double oil_rate, double gas_rate, double water_rate,
+    double alq
+) const
+{
+    const std::string msg = fmt::format("Group rate for {} : Well {} : "
+        "eff_factor = {}, oil_rate = {}, gas_rate = {}, water_rate = {}, "
+        "alq = {}, New group rates:  oil = {}, gas = {}, water = {}, alq = {}",
+        gr_name, well_name, eff_factor, well_oil_rate, well_gas_rate,
+        well_water_rate, well_alq, oil_rate, gas_rate, water_rate, alq);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugDisplayUpdatedGroupRates(
+      const std::string& name,
+      double oil_rate, double gas_rate, double water_rate, double alq) const
+{
+
+    const std::string msg = fmt::format("Updated group info for {} : "
+        "oil_rate = {}, gas_rate = {}, water_rate = {}, alq = {}",
+        name, oil_rate, gas_rate, water_rate, alq);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugEndInitializeGroup(const std::string& name) const
+{
+    const std::string msg = fmt::format("Finished with group {} ...", name);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugStartInitializeGroup(const std::string& name) const
+{
+    const std::string msg = fmt::format("Initializing group {} ...", name);
+    displayDebugMessage_(msg);
+}
+
 void
 GasLiftGroupInfo::
 displayDebugMessage_(const std::string &msg) const
@@ -417,6 +465,7 @@ GasLiftGroupInfo::
 initializeGroupRatesRecursive_(const Group &group)
 {
     std::array<double,4> rates{};
+    if (this->debug) debugStartInitializeGroup(group.name());
     auto& [oil_rate, water_rate, gas_rate, alq] = rates;
     if (group.wellgroup()) {
         for (const std::string& well_name : group.wells()) {
@@ -439,6 +488,13 @@ initializeGroupRatesRecursive_(const Group &group)
                     gas_rate += (factor * sw_gas_rate);
                     water_rate += (factor * sw_water_rate);
                     alq += (factor * sw_alq);
+                    if (this->debug) {
+                        debugDisplayWellContribution_(
+                            group.name(), well_name, factor,
+                            sw_oil_rate, sw_gas_rate, sw_water_rate, sw_alq,
+                            oil_rate, gas_rate, water_rate, alq
+                        );
+                    }
                 }
             }
         }
@@ -459,6 +515,7 @@ initializeGroupRatesRecursive_(const Group &group)
             alq += (gefac * sg_alq);
         }
     }
+    if (this->debug) debugEndInitializeGroup(group.name());
     std::optional<double> oil_target, gas_target, water_target, liquid_target, max_total_gas, max_alq;
     const auto controls = group.productionControls(this->summary_state_);
     if (group.has_control(Group::ProductionCMode::LRAT)) {
@@ -482,6 +539,10 @@ initializeGroupRatesRecursive_(const Group &group)
         updateGroupIdxMap_(group.name());
         this->group_rate_map_.try_emplace(group.name(),
             oil_rate, gas_rate, water_rate, alq, oil_target, gas_target, water_target, liquid_target, max_total_gas, max_alq);
+        if (this->debug) {
+            debugDisplayUpdatedGroupRates(
+                group.name(), oil_rate, gas_rate, water_rate, alq);
+        }
     }
     return std::make_tuple(oil_rate, gas_rate, water_rate, alq);
 }
diff --git a/opm/simulators/wells/GasLiftGroupInfo.hpp b/opm/simulators/wells/GasLiftGroupInfo.hpp
index 5bddc13d9..de60eda58 100644
--- a/opm/simulators/wells/GasLiftGroupInfo.hpp
+++ b/opm/simulators/wells/GasLiftGroupInfo.hpp
@@ -114,6 +114,18 @@ public:
 protected:
     bool checkDoGasLiftOptimization_(const std::string& well_name);
     bool checkNewtonIterationIdxOk_(const std::string& well_name);
+    void debugDisplayWellContribution_(
+        const std::string& gr_name, const std::string& well_name,
+        double eff_factor,
+        double well_oil_rate, double well_gas_rate, double well_water_rate,
+        double well_alq,
+        double oil_rate, double gas_rate, double water_rate,
+        double alq
+    ) const;
+    void debugDisplayUpdatedGroupRates(const std::string& name,
+        double oil_rate, double gas_rate, double water_rate, double alq) const;
+    void debugEndInitializeGroup(const std::string& name) const;
+    void debugStartInitializeGroup(const std::string& name) const;
     void displayDebugMessage_(const std::string& msg) const override;
     void displayDebugMessage_(const std::string& msg, const std::string& well_name);
     std::tuple<double, double, double> getProducerWellRates_(const int index);
diff --git a/opm/simulators/wells/GasLiftSingleWell.hpp b/opm/simulators/wells/GasLiftSingleWell.hpp
index 6dd32dfc4..19652008d 100644
--- a/opm/simulators/wells/GasLiftSingleWell.hpp
+++ b/opm/simulators/wells/GasLiftSingleWell.hpp
@@ -23,16 +23,9 @@
 #include <opm/models/utils/propertysystem.hh>
 #include <opm/models/utils/parametersystem.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
-
-// NOTE: StandardWell.hpp includes ourself (GasLiftSingleWell.hpp), so we need
-//   to forward declare StandardWell for it to be defined in this file.
-namespace Opm {
-    template<typename TypeTag> class StandardWell;
-}
-#include <opm/simulators/wells/StandardWell.hpp>
 #include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
 #include <opm/simulators/wells/GasLiftGroupInfo.hpp>
-
+#include <opm/simulators/wells/WellInterface.hpp>
 
 #include <optional>
 #include <vector>
@@ -45,12 +38,11 @@ namespace Opm
     class GasLiftSingleWell : public GasLiftSingleWellGeneric
     {
         using Simulator = GetPropType<TypeTag, Properties::Simulator>;
-        using StdWell = StandardWell<TypeTag>;
         using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
 
     public:
         GasLiftSingleWell(
-            const StdWell &std_well,
+            const WellInterface<TypeTag> &well,
             const Simulator &ebos_simulator,
             const SummaryState &summary_state,
             DeferredLogger &deferred_logger,
@@ -60,7 +52,7 @@ namespace Opm
             GLiftSyncGroups &sync_groups,
             bool glift_debug
         );
-        const WellInterfaceGeneric &getStdWell() const override { return std_well_; }
+        const WellInterfaceGeneric &getWell() const override { return well_; }
 
     private:
         std::optional<double> computeBhpAtThpLimit_(double alq) const override;
@@ -69,7 +61,7 @@ namespace Opm
         void setAlqMaxRate_(const GasLiftOpt::Well& well);
 
         const Simulator &ebos_simulator_;
-        const StdWell &std_well_;
+        const WellInterface<TypeTag> &well_;
     };
 
 } // namespace Opm
diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
index 1e77480bd..196c07a6f 100644
--- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
+++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
@@ -285,7 +285,13 @@ checkThpControl_() const
     const int well_index = this->well_state_.index(this->well_name_).value();
     const Well::ProducerCMode& control_mode =
                          this->well_state_.well(well_index).production_cmode;
-    return control_mode == Well::ProducerCMode::THP;
+    bool thp_control = control_mode == Well::ProducerCMode::THP;
+    if (this->debug) {
+        if (!thp_control) {
+            displayDebugMessage_("Well is not under THP control, skipping iteration..");
+        }
+    }
+    return thp_control;
 }
 
 
@@ -431,6 +437,18 @@ debugShowLimitingTargets_(const LimitedRates& rates) const
     }
 }
 
+void
+GasLiftSingleWellGeneric::
+debugShowProducerControlMode() const
+{
+    const int well_index = this->well_state_.index(this->well_name_).value();
+    const Well::ProducerCMode& control_mode =
+                         this->well_state_.well(well_index).production_cmode;
+    const std::string msg = fmt::format("Current control mode is: {}",
+        Well::ProducerCMode2String(control_mode));
+    displayDebugMessage_(msg);
+}
+
 void
 GasLiftSingleWellGeneric::
 debugShowStartIteration_(double alq, bool increase, double oil_rate)
@@ -1139,6 +1157,7 @@ std::unique_ptr<GasLiftWellState>
 GasLiftSingleWellGeneric::
 runOptimizeLoop_(bool increase)
 {
+    if (this->debug) debugShowProducerControlMode();
     std::unique_ptr<GasLiftWellState> ret_value; // nullptr initially
     auto rates = getInitialRatesWithLimit_();
     if (!rates) return ret_value;
@@ -1160,7 +1179,11 @@ runOptimizeLoop_(bool increase)
     }
 
     OptimizeState state {*this, increase};
-    if (!checkThpControl_()) {
+    auto temp_alq = cur_alq;
+    if (checkThpControl_()) {
+        if (this->debug) debugShowStartIteration_(temp_alq, increase, new_rates.oil);
+    }
+    else {
         // If the well is not under THP control, we can still use the previous
         //   initial adjustment of ALQ by using the well's THP limit to calculate
         //   BHP and then well rates from that.
@@ -1169,8 +1192,6 @@ runOptimizeLoop_(bool increase)
         //   Then gaslift can be reduced while still keeping the group target.
         state.stop_iteration = true;
     }
-    auto temp_alq = cur_alq;
-    if (this->debug) debugShowStartIteration_(temp_alq, increase, new_rates.oil);
     while (!state.stop_iteration && (++state.it <= this->max_iterations_)) {
         if (state.checkRatesViolated(new_rates)) break;
         if (state.checkAlqOutsideLimits(temp_alq, new_rates.oil)) break;
@@ -1524,11 +1545,11 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate)
     else { // we are decreasing lift gas
         if ( alq == 0 ) {
             ss << "ALQ is zero, cannot decrease further. Stopping iteration.";
-            return true;
+            result = true;
         }
         else if ( alq < 0 ) {
             ss << "Negative ALQ: " << alq << ". Stopping iteration.";
-            return true;
+            result = true;
         }
         // NOTE: A negative min_alq_ means: allocate at least enough lift gas
         //  to enable the well to flow, see WLIFTOPT item 5.
diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
index 10d8534ba..f4117096c 100644
--- a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
+++ b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
@@ -96,7 +96,7 @@ public:
 
     std::unique_ptr<GasLiftWellState> runOptimize(const int iteration_idx);
 
-    virtual const WellInterfaceGeneric& getStdWell() const = 0;
+    virtual const WellInterfaceGeneric& getWell() const = 0;
 
 protected:
     GasLiftSingleWellGeneric(
@@ -260,6 +260,7 @@ protected:
     void debugShowAlqIncreaseDecreaseCounts_();
     void debugShowBhpAlqTable_();
     void debugShowLimitingTargets_(const LimitedRates& rates) const;
+    void debugShowProducerControlMode() const;
     void debugShowStartIteration_(double alq, bool increase, double oil_rate);
     void debugShowTargets_();
     void displayDebugMessage_(const std::string& msg) const override;
diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp
index 76c7e08b9..40fcb864f 100644
--- a/opm/simulators/wells/GasLiftSingleWell_impl.hpp
+++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp
@@ -21,7 +21,7 @@ namespace Opm {
 
 template<typename TypeTag>
 GasLiftSingleWell<TypeTag>::
-GasLiftSingleWell(const StdWell &std_well,
+GasLiftSingleWell(const WellInterface<TypeTag> &well,
                   const Simulator &ebos_simulator,
                   const SummaryState &summary_state,
                   DeferredLogger &deferred_logger,
@@ -37,17 +37,17 @@ GasLiftSingleWell(const StdWell &std_well,
         deferred_logger,
         well_state,
         group_state,
-        std_well.wellEcl(),
+        well.wellEcl(),
         summary_state,
         group_info,
-        std_well.phaseUsage(),
+        well.phaseUsage(),
         ebos_simulator.vanguard().schedule(),
         ebos_simulator.episodeIndex(),
         sync_groups,
         glift_debug
     )
    , ebos_simulator_{ebos_simulator}
-   , std_well_{std_well}
+   , well_{well}
 {
     const auto& gl_well = *gl_well_;
     if(useFixedAlq_(gl_well)) {
@@ -59,7 +59,7 @@ GasLiftSingleWell(const StdWell &std_well,
         this->optimize_ = true;
     }
 
-    const auto& pu = std_well_.phaseUsage();
+    const auto& pu = well_.phaseUsage();
     this->oil_pos_ = pu.phase_pos[Oil];
     this->gas_pos_ = pu.phase_pos[Gas];
     this->water_pos_ = pu.phase_pos[Water];
@@ -103,7 +103,7 @@ GasLiftSingleWell<TypeTag>::
 computeWellRates_( double bhp, bool bhp_is_limited, bool debug_output ) const
 {
     std::vector<double> potentials(this->num_phases_, 0.0);
-    this->std_well_.computeWellRatesWithBhp(
+    this->well_.computeWellRatesWithBhp(
         this->ebos_simulator_, bhp, potentials, this->deferred_logger_);
     if (debug_output) {
         const std::string msg = fmt::format("computed well potentials given bhp {}, "
@@ -128,7 +128,7 @@ std::optional<double>
 GasLiftSingleWell<TypeTag>::
 computeBhpAtThpLimit_(double alq) const
 {
-    auto bhp_at_thp_limit = this->std_well_.computeBhpAtThpLimitProdWithAlq(
+    auto bhp_at_thp_limit = this->well_.computeBhpAtThpLimitProdWithAlq(
         this->ebos_simulator_,
         this->summary_state_,
         this->deferred_logger_,
@@ -167,7 +167,7 @@ setAlqMaxRate_(const GasLiftOpt::Well &well)
         // According to the manual for WLIFTOPT, item 3:
         //   The default value should be set to the largest ALQ
         //   value in the well's VFP table
-        const auto& table = std_well_.vfp_properties_->getProd()->getTable(
+        const auto& table = well_.vfpProperties()->getProd()->getTable(
                 this->controls_.vfp_table_number);
         const auto& alq_values = table.getALQAxis();
         // Assume the alq_values are sorted in ascending order, so
diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp
index bf702da33..9ea089b73 100644
--- a/opm/simulators/wells/GasLiftStage2.cpp
+++ b/opm/simulators/wells/GasLiftStage2.cpp
@@ -349,7 +349,7 @@ getCurrentWellRates_(const std::string &well_name, const std::string &group_name
     std::string debug_info;
     if (this->stage1_wells_.count(well_name) == 1) {
         GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(well_name).get());
-        const WellInterfaceGeneric &well = gs_well.getStdWell();
+        const WellInterfaceGeneric &well = gs_well.getWell();
         well_ptr = &well;
         GasLiftWellState &state = *(this->well_state_map_.at(well_name).get());
         std::tie(oil_rate, gas_rate) = state.getRates();
@@ -358,7 +358,7 @@ getCurrentWellRates_(const std::string &well_name, const std::string &group_name
     }
     else if (this->prod_wells_.count(well_name) == 1) {
         well_ptr = this->prod_wells_.at(well_name);
-        std::tie(oil_rate, gas_rate) = getStdWellRates_(*well_ptr);
+        std::tie(oil_rate, gas_rate) = getWellRates_(*well_ptr);
         success = true;
         if ( this->debug) debug_info = "(B)";
     }
@@ -413,7 +413,7 @@ getCurrentWellRates_(const std::string &well_name, const std::string &group_name
 
 std::pair<double, double>
 GasLiftStage2::
-getStdWellRates_(const WellInterfaceGeneric &well)
+getWellRates_(const WellInterfaceGeneric &well)
 {
     const int well_index = well.indexOfWell();
     const auto& ws = this->well_state_.well(well_index);
@@ -890,7 +890,7 @@ checkAtLeastTwoWells(std::vector<GasLiftSingleWell *> &wells)
 {
     int numberOfwells = 0;
     for (auto well : wells){
-        int index_of_wells = well->getStdWell().indexOfWell();
+        int index_of_wells = well->getWell().indexOfWell();
         if (!this->parent.well_state_.wellIsOwned(index_of_wells, well->name()))
             continue;
         numberOfwells++;
@@ -1107,7 +1107,7 @@ updateRates(const std::string &well_name)
         const GradInfo &gi = this->parent.dec_grads_.at(well_name);
         GasLiftWellState &state = *(this->parent.well_state_map_.at(well_name).get());
         GasLiftSingleWell &gs_well = *(this->parent.stage1_wells_.at(well_name).get());
-        const WellInterfaceGeneric &well = gs_well.getStdWell();
+        const WellInterfaceGeneric &well = gs_well.getWell();
         // only get deltas for wells owned by this rank
         if (this->parent.well_state_.wellIsOwned(well.indexOfWell(), well_name)) {
             const auto &well_ecl = well.wellEcl();
diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp
index dfeaa8680..63562aaa0 100644
--- a/opm/simulators/wells/GasLiftStage2.hpp
+++ b/opm/simulators/wells/GasLiftStage2.hpp
@@ -96,7 +96,7 @@ protected:
         const Group& group);
     void getGroupGliftWellsRecursive_(
         const Group& group, std::vector<GasLiftSingleWell *>& wells);
-    std::pair<double, double> getStdWellRates_(const WellInterfaceGeneric& well);
+    std::pair<double, double> getWellRates_(const WellInterfaceGeneric& well);
     void optimizeGroup_(const Group& group);
     void optimizeGroupsRecursive_(const Group& group);
     void recalculateGradientAndUpdateData_(
diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp
index 699331f36..656e1f02d 100644
--- a/opm/simulators/wells/MultisegmentWell.hpp
+++ b/opm/simulators/wells/MultisegmentWell.hpp
@@ -52,11 +52,6 @@ namespace Opm
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
         using typename Base::FluidState;
-        using typename Base::GasLiftSingleWell;
-        using typename Base::GLiftProdWells;
-        using typename Base::GLiftOptWells;
-        using typename Base::GLiftWellStateMap;
-        using typename Base::GLiftSyncGroups;
 
         using Base::has_solvent;
         using Base::has_polymer;
@@ -99,21 +94,6 @@ namespace Opm
 
         virtual void initPrimaryVariablesEvaluation() const override;
 
-        virtual void gasLiftOptimizationStage1 (
-            WellState&,
-            const GroupState&,
-            const Simulator&,
-            DeferredLogger&,
-            GLiftProdWells &,
-            GLiftOptWells &,
-            GLiftWellStateMap &,
-            GasLiftGroupInfo &,
-            GLiftSyncGroups &,
-            bool
-        ) const override {
-            // Not implemented yet
-        }
-
         /// updating the well state based the current control mode
         virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
                                                const GroupState& group_state,
@@ -173,6 +153,12 @@ namespace Opm
                                     double* connII,
                                     DeferredLogger& deferred_logger) const;
 
+        virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
+            const Simulator& ebos_simulator,
+            const SummaryState& summary_state,
+            DeferredLogger& deferred_logger,
+            double alq_value) const override;
+
     protected:
         int number_segments_;
 
@@ -258,19 +244,20 @@ namespace Opm
                                         std::vector<double>& well_flux,
                                         DeferredLogger& deferred_logger) const;
 
-        void computeWellRatesWithBhp(const Simulator& ebosSimulator,
-                                     const Scalar& bhp,
+        virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
+                                     const double& bhp,
                                      std::vector<double>& well_flux,
-                                     DeferredLogger& deferred_logger) const;
+                                     DeferredLogger& deferred_logger) const override;
 
         void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
                                                const Scalar& bhp,
                                                std::vector<double>& well_flux,
                                                DeferredLogger& deferred_logger) const;
 
-        std::vector<double>
-        computeWellPotentialWithTHP(const Simulator& ebos_simulator,
-                                    DeferredLogger& deferred_logger) const;
+        std::vector<double> computeWellPotentialWithTHP(
+                                 const WellState& well_state,
+                                 const Simulator& ebos_simulator,
+                                 DeferredLogger& deferred_logger) const;
 
         virtual double getRefDensity() const override;
 
@@ -306,9 +293,12 @@ namespace Opm
         bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
 
 
-        std::optional<double> computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
-                                                       const SummaryState& summary_state,
-                                                       DeferredLogger& deferred_logger) const;
+
+        std::optional<double> computeBhpAtThpLimitProd(
+            const WellState& well_state,
+            const Simulator& ebos_simulator,
+            const SummaryState& summary_state,
+            DeferredLogger& deferred_logger) const;
 
         std::optional<double> computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
                                                       const SummaryState& summary_state,
diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp
index 252fffffe..48fa91054 100644
--- a/opm/simulators/wells/MultisegmentWellGeneric.cpp
+++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp
@@ -436,11 +436,13 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
 template<typename Scalar>
 std::optional<double>
 MultisegmentWellGeneric<Scalar>::
-computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
-                         const SummaryState& summary_state,
-                         const double maxPerfPress,
-                         const double rho,
-                         DeferredLogger& deferred_logger) const
+computeBhpAtThpLimitProdWithAlq(
+                   const std::function<std::vector<double>(const double)>& frates,
+                   const SummaryState& summary_state,
+                   const double maxPerfPress,
+                   const double rho,
+                   DeferredLogger& deferred_logger,
+                   double alq_value) const
 {
     // Given a VFP function returning bhp as a function of phase
     // rates and thp:
@@ -468,10 +470,11 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
     const double vfp_ref_depth = table.getDatumDepth();
     const double thp_limit = baseif_.getTHPConstraint(summary_state);
     const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
-    auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
+    auto fbhp = [this, &controls, thp_limit, dp, alq_value](
+                                  const std::vector<double>& rates) {
         assert(rates.size() == 3);
         return baseif_.vfpProperties()->getProd()
-        ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp;
+        ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - dp;
     };
 
     // Make the flo() function.
diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp
index d68c51c4a..9c8d2b47e 100644
--- a/opm/simulators/wells/MultisegmentWellGeneric.hpp
+++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp
@@ -70,11 +70,13 @@ protected:
                                                   const double rho,
                                                   DeferredLogger& deferred_logger) const;
 
-    std::optional<double> computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
-                                                   const SummaryState& summary_state,
-                                                   const double maxPerfPress,
-                                                   const double rho,
-                                                   DeferredLogger& deferred_logger) const;
+    std::optional<double> computeBhpAtThpLimitProdWithAlq(
+        const std::function<std::vector<double>(const double)>& frates,
+        const SummaryState& summary_state,
+        const double maxPerfPress,
+        const double rho,
+        DeferredLogger& deferred_logger,
+        double alq_value) const;
 
     std::optional<double> bhpMax(const std::function<double(const double)>& fflo,
                                  const double bhp_limit,
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 9a11be6dc..893385b50 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -300,7 +300,8 @@ namespace Opm
         if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
             computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
         } else {
-            well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
+            well_potentials = computeWellPotentialWithTHP(
+                well_state, ebosSimulator, deferred_logger);
         }
         deferred_logger.debug("Cost in iterations of finding well potential for well "
                               + this->name() + ": " + std::to_string(debug_cost_counter_));
@@ -342,7 +343,7 @@ namespace Opm
     void
     MultisegmentWell<TypeTag>::
     computeWellRatesWithBhp(const Simulator& ebosSimulator,
-                            const Scalar& bhp,
+                            const double& bhp,
                             std::vector<double>& well_flux,
                             DeferredLogger& deferred_logger) const
     {
@@ -448,8 +449,10 @@ namespace Opm
     template<typename TypeTag>
     std::vector<double>
     MultisegmentWell<TypeTag>::
-    computeWellPotentialWithTHP(const Simulator& ebos_simulator,
-                                DeferredLogger& deferred_logger) const
+    computeWellPotentialWithTHP(
+          const WellState& well_state,
+          const Simulator& ebos_simulator,
+          DeferredLogger& deferred_logger) const
     {
         std::vector<double> potentials(this->number_of_phases_, 0.0);
         const auto& summary_state = ebos_simulator.vanguard().summaryState();
@@ -472,7 +475,8 @@ namespace Opm
                 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
             }
         } else {
-            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
+            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
+                  well_state, ebos_simulator, summary_state, deferred_logger);
             if (bhp_at_thp_limit) {
                 const auto& controls = well.productionControls(summary_state);
                 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
@@ -1293,11 +1297,16 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellState& /*well_state*/, DeferredLogger& deferred_logger)
+    checkOperabilityUnderTHPLimit(
+             const Simulator& ebos_simulator,
+             const WellState& well_state,
+             DeferredLogger& deferred_logger)
     {
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
-        const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger)
-        : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
+        const auto obtain_bhp = this->isProducer()
+            ? computeBhpAtThpLimitProd(
+                        well_state, ebos_simulator, summaryState, deferred_logger)
+            : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
 
         if (obtain_bhp) {
             this->operability_status_.can_obtain_bhp_with_thp_limit = true;
@@ -1734,15 +1743,30 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    std::optional<double>
+    MultisegmentWell<TypeTag>::
+    computeBhpAtThpLimitProd(const WellState& well_state,
+                             const Simulator& ebos_simulator,
+                             const SummaryState& summary_state,
+                             DeferredLogger& deferred_logger) const
+    {
+        return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
+                                               ebos_simulator,
+                                               summary_state,
+                                               deferred_logger,
+                                               this->getALQ(well_state));
+    }
 
 
 
     template<typename TypeTag>
     std::optional<double>
     MultisegmentWell<TypeTag>::
-    computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
-                             const SummaryState& summary_state,
-                             DeferredLogger& deferred_logger) const
+    computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
+                                    const SummaryState& summary_state,
+                                    DeferredLogger& deferred_logger,
+                                    double alq_value) const
     {
         // Make the frates() function.
         auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -1757,11 +1781,12 @@ namespace Opm
         };
 
         auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
-               computeBhpAtThpLimitProd(frates,
+               computeBhpAtThpLimitProdWithAlq(frates,
                                         summary_state,
                                         maxPerfPress(ebos_simulator),
                                         getRefDensity(),
-                                        deferred_logger);
+                                        deferred_logger,
+                                        alq_value);
 
        if(bhpAtLimit)
            return bhpAtLimit;
@@ -1776,11 +1801,12 @@ namespace Opm
        };
 
        return this->MultisegmentWellGeneric<Scalar>::
-              computeBhpAtThpLimitProd(fratesIter,
+              computeBhpAtThpLimitProdWithAlq(fratesIter,
                                        summary_state,
                                        maxPerfPress(ebos_simulator),
                                        getRefDensity(),
-                                       deferred_logger);
+                                       deferred_logger,
+                                       alq_value);
 
     }
 
diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp
index 3c3101e10..add189b0d 100644
--- a/opm/simulators/wells/StandardWell.hpp
+++ b/opm/simulators/wells/StandardWell.hpp
@@ -31,8 +31,6 @@
 #include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
-#include <opm/simulators/wells/GasLiftSingleWell.hpp>
-#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
 
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
@@ -83,11 +81,6 @@ namespace Opm
         using typename Base::SparseMatrixAdapter;
         using typename Base::FluidState;
         using typename Base::RateVector;
-        using typename Base::GasLiftSingleWell;
-        using typename Base::GLiftOptWells;
-        using typename Base::GLiftProdWells;
-        using typename Base::GLiftWellStateMap;
-        using typename Base::GLiftSyncGroups;
 
         using Base::has_solvent;
         using Base::has_zFraction;
@@ -202,19 +195,6 @@ namespace Opm
             return this->param_.matrix_add_well_contributions_;
         }
 
-        virtual void gasLiftOptimizationStage1 (
-            WellState& well_state,
-            const GroupState& group_state,
-            const Simulator& ebosSimulator,
-            DeferredLogger& deferred_logger,
-            GLiftProdWells &prod_wells,
-            GLiftOptWells &glift_wells,
-            GLiftWellStateMap &state_map,
-            GasLiftGroupInfo &group_info,
-            GLiftSyncGroups &sync_groups,
-            bool glift_debug
-        ) const override;
-
         /* returns BHP */
         double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
                                const SummaryState &summary_state,
@@ -229,19 +209,17 @@ namespace Opm
             std::vector<double> &potentials,
             double alq) const;
 
-        // NOTE: Cannot be protected since it is used by GasLiftRuntime
-        std::optional<double> computeBhpAtThpLimitProdWithAlq(
+        virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
             const Simulator& ebos_simulator,
             const SummaryState& summary_state,
             DeferredLogger& deferred_logger,
-            double alq_value) const;
+            double alq_value) const override;
 
-        // NOTE: Cannot be protected since it is used by GasLiftRuntime
-        void computeWellRatesWithBhp(
+        virtual void computeWellRatesWithBhp(
             const Simulator& ebosSimulator,
             const double& bhp,
             std::vector<double>& well_flux,
-            DeferredLogger& deferred_logger) const;
+            DeferredLogger& deferred_logger) const override;
 
         // NOTE: These cannot be protected since they are used by GasLiftRuntime
         using Base::phaseUsage;
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index eee8c5cb0..c717c92ad 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -1812,38 +1812,6 @@ namespace Opm
                                              alq);
     }
 
-    template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
-    gasLiftOptimizationStage1(
-                       WellState& well_state,
-                       const GroupState& group_state,
-                       const Simulator& ebos_simulator,
-                       DeferredLogger& deferred_logger,
-                       GLiftProdWells &prod_wells,
-                       GLiftOptWells &glift_wells,
-                       GLiftWellStateMap &glift_state_map,
-                       GasLiftGroupInfo &group_info,
-                       GLiftSyncGroups &sync_groups,
-                       bool glift_debug
-    ) const
-    {
-        const auto& summary_state = ebos_simulator.vanguard().summaryState();
-        std::unique_ptr<GasLiftSingleWell> glift
-            = std::make_unique<GasLiftSingleWell>(
-                *this, ebos_simulator, summary_state,
-                deferred_logger, well_state, group_state, group_info,
-                sync_groups, glift_debug);
-        auto state = glift->runOptimize(
-            ebos_simulator.model().newtonMethod().numIterations());
-        if (state) {
-            glift_state_map.insert({this->name(), std::move(state)});
-            glift_wells.insert({this->name(), std::move(glift)});
-            return;
-        }
-        prod_wells.insert({this->name(), this});
-    }
-
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 25871e13e..7d35989b9 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -57,6 +57,7 @@ namespace Opm {
 #include <opm/material/densead/Evaluation.hpp>
 
 #include <opm/simulators/wells/WellInterfaceIndices.hpp>
+#include <opm/simulators/timestepping/ConvergenceReport.hpp>
 
 #include <cassert>
 #include <vector>
@@ -159,17 +160,18 @@ public:
                         const GroupState& group_state,
                         DeferredLogger& deferred_logger);
 
-    virtual void gasLiftOptimizationStage1 (
-        WellState& well_state,
-        const GroupState& group_state,
+    virtual void computeWellRatesWithBhp(
         const Simulator& ebosSimulator,
+        const double& bhp,
+        std::vector<double>& well_flux,
+        DeferredLogger& deferred_logger
+    ) const = 0;
+
+    virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
+        const Simulator& ebos_simulator,
+        const SummaryState& summary_state,
         DeferredLogger& deferred_logger,
-        GLiftProdWells& prod_wells,
-        GLiftOptWells& glift_wells,
-        GLiftWellStateMap& state_map,
-        GasLiftGroupInfo &group_info,
-        GLiftSyncGroups &sync_groups,
-        bool glift_debug
+        double alq_value
     ) const = 0;
 
     /// using the solution x to recover the solution xw for wells and applying