diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index fef06c75f18..719b6d62f52 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -1385,9 +1385,9 @@ updateAndCommunicateGroupData(const int reportStepIdx,
             calcResvCoeff(fipnum, pvtreg, this->groupState().production_rates(group.name()), resv_coeff);
             const Scalar efficiencyFactor = well->wellEcl().getEfficiencyFactor() *
                                     ws.efficiency_scaling_factor;
-            // Translate injector type from control to Phase.
-            Scalar group_target = std::numeric_limits<Scalar>::max();
             if (well->isProducer()) {
+                std::pair<WellProducerCMode, Scalar> group_target;
+                group_target.second = std::numeric_limits<Scalar>::max();
                 group_target = wg_helper.getWellGroupTargetProducer(
                     well->name(),
                     well->wellEcl().groupName(),
@@ -1397,6 +1397,9 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                     resv_coeff,
                     deferred_logger
                 );
+                auto& ws_update = this->wellState().well(well->indexOfWell());
+                ws_update.production_cmode_group_translated = group_target.first;
+                ws_update.group_target = group_target.second;
             } else {
                 const auto& well_controls = well->wellEcl().injectionControls(summaryState_);
                 auto injectorType = well_controls.injector_type;
@@ -1420,6 +1423,7 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                 default:
                     throw std::logic_error("MULTI-phase injection is not supported, but was requested for well " + well->name());
                 }
+                std::pair<WellInjectorCMode, Scalar> group_target;
                 group_target = wg_helper.getWellGroupTargetInjector(
                     well->name(),
                     well->wellEcl().groupName(),
@@ -1430,9 +1434,10 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                     resv_coeff,
                     deferred_logger
                 );
+                auto& ws_update = this->wellState().well(well->indexOfWell());
+                ws_update.injection_cmode_group_translated = group_target.first;
+                ws_update.group_target = group_target.second;
             }
-            auto& ws_update = this->wellState().well(well->indexOfWell());
-            ws_update.group_target = group_target;
         }
     }
 }
diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
index b1b585d9126..512d4d0f8c0 100644
--- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
+++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
@@ -681,6 +681,53 @@ getWQTotal() const
     return evaluation_[0][WQTotal];
 }
 
+template<class FluidSystem, class Indices>
+void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
+fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const
+{
+    surface_fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        surface_fractions[oil_pos] = 1.0;
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            surface_fractions[water_pos] = value_[0][WFrac];
+            surface_fractions[oil_pos] -= surface_fractions[water_pos];
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[0][GFrac];
+            surface_fractions[oil_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        surface_fractions[water_pos] = 1.0;
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[0][GFrac];
+            surface_fractions[water_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        surface_fractions[gas_pos] = 1.0;
+    }
+    for (int p = 0; p < well_.numPhases(); ++p) {
+        const Scalar scale = well_.scalingFactor(p);
+        // for injection wells, there should only one non-zero scaling factor
+        if (scale > 0.) {
+            surface_fractions[p] /= scale;
+        } else {
+            // this should only happen to injection wells
+            surface_fractions[p] = 0.;
+        }
+        surface_fractions[p] = std::max(surface_fractions[p], Scalar{0.0});
+    }
+}
+
 template<class FluidSystem, class Indices>
 void
 MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
index 590f87c7b83..3bf897219eb 100644
--- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
+++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
@@ -143,6 +143,9 @@ class MultisegmentWellPrimaryVariables
     void setValue(const int idx, const std::array<Scalar, numWellEq>& val)
     { value_[idx] = val; }
 
+    //! \brief Fetch surface fractions for the well.
+    void fetchWellSurfaceFractions(std::vector<Scalar> &surface_fractions) const;
+
     //! output the segments with pressure close to lower pressure limit for debugging purpose
     void outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const;
 
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index da018bf1569..53a4eb1a798 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -44,6 +44,7 @@
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
 
 #include <algorithm>
 #include <cstddef>
@@ -1638,8 +1639,6 @@ namespace Opm
                 return false;
         }
 
-        updatePrimaryVariables(simulator, well_state, deferred_logger);
-
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<Scalar> measure_history;
         int it = 0;
@@ -1664,7 +1663,8 @@ namespace Opm
         const auto operability_orig = this->operability_status_;
         auto well_status_cur = well_status_orig;
         // don't allow opening wells that has a stopped well status
-        const bool allow_open = well_state.well(this->index_of_well_).status == WellStatus::OPEN;
+        auto& ws = well_state.well(this->index_of_well_); 
+        const bool allow_open = ws.status == WellStatus::OPEN;
         // don't allow switcing for wells under zero rate target or requested fixed status and control
         const bool allow_switching = !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
                                      (!fixed_control || !fixed_status) && allow_open;
@@ -1672,14 +1672,26 @@ namespace Opm
         // well needs to be set operable or else solving/updating of re-opened wells is skipped
         this->operability_status_.resetOperability();
         this->operability_status_.solvable = true;
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+        // keep track of surface rate fractions for checking control feasibility
+        std::vector<Scalar> surface_fractions;
+        bool feasible_constraint = true;
 
         for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
             ++its_since_last_switch;
-            if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
+            if (this->isProducer()) {
+                this->primary_variables_.fetchWellSurfaceFractions(surface_fractions);
+                feasible_constraint = this->isFeasibleProductionConstraint(ws, surface_fractions, ws.production_cmode);
+            }
+            const bool check_controls = allow_switching &&
+                                        its_since_last_switch >= min_its_after_switch &&
+                                        status_switch_count < max_status_switch;
+            // if constraint is not feasible, we force a switch to prevent singularity
+            if (check_controls || !feasible_constraint) {
                 const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
                 bool changed = this->updateWellControlAndStatusLocalIteration(
-                    simulator, wgHelper, inj_controls, prod_controls, wqTotal,
-                    well_state, deferred_logger, fixed_control, fixed_status
+                    simulator, wgHelper, inj_controls, prod_controls, wqTotal, well_state,
+                    surface_fractions, deferred_logger, fixed_control, fixed_status
                 );
                 if (changed) {
                     its_since_last_switch = 0;
diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
index 0894bcd6019..6f7f57fae0a 100644
--- a/opm/simulators/wells/SingleWellState.cpp
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -53,6 +53,7 @@ SingleWellState(const std::string& name_,
     , prev_surface_rates(pu.numActivePhases())
     , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu.numActivePhases())
     , trivial_group_target(false)
+    , prevent_group_control(false)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
@@ -410,7 +411,10 @@ bool SingleWellState<Scalar, IndexTraits>::operator==(const SingleWellState& rhs
            this->production_cmode == rhs.production_cmode &&
            this->alq_state == rhs.alq_state &&
            this->primaryvar == rhs.primaryvar &&
-           this->group_target == rhs.group_target;
+           this->group_target == rhs.group_target &&
+           this->injection_cmode_group_translated == rhs.injection_cmode_group_translated &&
+           this->production_cmode_group_translated == rhs.production_cmode_group_translated &&
+           this->prevent_group_control == rhs.prevent_group_control;
 }
 
 template class SingleWellState<double, BlackOilDefaultFluidSystemIndices>;
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index edff8960896..8214954772b 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -85,6 +85,9 @@ class SingleWellState {
         serializer(primaryvar);
         serializer(alq_state);
         serializer(group_target);
+        serializer(injection_cmode_group_translated);
+        serializer(production_cmode_group_translated);
+        serializer(prevent_group_control);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -125,6 +128,9 @@ class SingleWellState {
     PerfData<Scalar> perf_data;
     bool trivial_group_target;
     std::optional<Scalar> group_target;
+    std::optional<WellInjectorCMode> injection_cmode_group_translated{WellInjectorCMode::CMODE_UNDEFINED};
+    std::optional<WellProducerCMode> production_cmode_group_translated{WellProducerCMode::CMODE_UNDEFINED};
+    bool prevent_group_control;
     SegmentState<Scalar> segments;
     Events events;
     WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};
diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp
index 4aba82cc3ea..b84afd8a912 100644
--- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp
+++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp
@@ -760,6 +760,53 @@ checkFinite(DeferredLogger& deferred_logger) const
     }
 }
 
+template<class FluidSystem, class Indices>
+void StandardWellPrimaryVariables<FluidSystem,Indices>::
+fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const
+{
+    surface_fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        surface_fractions[oil_pos] = 1.0;
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            surface_fractions[water_pos] = value_[WFrac];
+            surface_fractions[oil_pos] -= surface_fractions[water_pos];
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[GFrac];
+            surface_fractions[oil_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        surface_fractions[water_pos] = 1.0;
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[GFrac];
+            surface_fractions[water_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        surface_fractions[gas_pos] = 1.0;
+    }
+    for (int p = 0; p < well_.numPhases(); ++p) {
+        const Scalar scale = well_.scalingFactor(p);
+        // for injection wells, there should be only one non-zero scaling factor
+        if (scale > 0.) {
+            surface_fractions[p] /= scale;
+        } else {
+            // this should only happen to injection wells
+            surface_fractions[p] = 0.;
+        }
+        surface_fractions[p] = std::max(surface_fractions[p], Scalar{0.0});
+    }
+}
+
 #include <opm/simulators/utils/InstantiationIndicesMacros.hpp>
 
 INSTANTIATE_TYPE_INDICES(StandardWellPrimaryVariables, double)
diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp
index 16cf78c0612..fbc82cd8dda 100644
--- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp
+++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp
@@ -148,6 +148,9 @@ class StandardWellPrimaryVariables {
     void setValue(const int idx, const Scalar val)
     { value_[idx] = val; }
 
+    //! \brief Fetch well surface fraction values.
+    void fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const;
+
 private:
     //! \brief Initialize evaluations from values.
     void setEvaluationsFromValues();
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 56d31287e21..01870d58f48 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -37,6 +37,7 @@
 #include <opm/simulators/wells/VFPHelpers.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/wells/WellConvergence.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
 
 #include <algorithm>
 #include <cstddef>
@@ -2444,7 +2445,6 @@ namespace Opm
                                const bool fixed_status /*false*/)
     {
         const auto& group_state = wgHelper.groupState();
-        updatePrimaryVariables(simulator, well_state, deferred_logger);
 
         const int max_iter = this->param_.max_inner_iter_wells_;
         int it = 0;
@@ -2467,7 +2467,8 @@ namespace Opm
         auto well_status_cur = well_status_orig;
         int status_switch_count = 0;
         // don't allow opening wells that has a stopped well status
-        const bool allow_open = well_state.well(this->index_of_well_).status == WellStatus::OPEN;
+        auto& ws = well_state.well(this->index_of_well_);
+        const bool allow_open = ws.status == WellStatus::OPEN;
         // don't allow switcing for wells under zero rate target or requested fixed status and control
         const bool allow_switching =
             !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
@@ -2478,13 +2479,25 @@ namespace Opm
         // well needs to be set operable or else solving/updating of re-opened wells is skipped
         this->operability_status_.resetOperability();
         this->operability_status_.solvable = true;
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+        // keep track of surface rate fractions for checking control feasibility
+        std::vector<Scalar> surface_fractions;
+        bool feasible_constraint = true;
         do {
             its_since_last_switch++;
-            if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
+            if (this->isProducer()) {
+                this->primary_variables_.fetchWellSurfaceFractions(surface_fractions);
+                feasible_constraint = this->isFeasibleProductionConstraint(ws, surface_fractions, ws.production_cmode);
+            }
+            const bool check_controls = allow_switching &&
+                                        its_since_last_switch >= min_its_after_switch &&
+                                        status_switch_count < max_status_switch;
+            // if constraint is not feasible, we force a switch to prevent singularity
+            if (check_controls || !feasible_constraint){
                 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
                 changed = this->updateWellControlAndStatusLocalIteration(
-                    simulator, wgHelper, inj_controls, prod_controls, wqTotal,
-                    well_state, deferred_logger, fixed_control, fixed_status
+                    simulator, wgHelper, inj_controls, prod_controls, wqTotal, well_state,
+                    surface_fractions, deferred_logger, fixed_control, fixed_status
                 );
                 if (changed){
                     its_since_last_switch = 0;
diff --git a/opm/simulators/wells/WellConstraints.cpp b/opm/simulators/wells/WellConstraints.cpp
index 9d73bc43a7a..b5a936933df 100644
--- a/opm/simulators/wells/WellConstraints.cpp
+++ b/opm/simulators/wells/WellConstraints.cpp
@@ -322,6 +322,328 @@ activeProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
     return currentControl;
 }
 
+template<typename Scalar, typename IndexTraits>
+bool
+WellConstraints<Scalar, IndexTraits>::
+isFeasibleProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                               const std::vector<Scalar>& surface_fractions,
+                               const Well::ProducerCMode cmode_input) const
+{
+    auto cmode = cmode_input;
+    if (cmode_input == Well::ProducerCMode::GRUP) {
+        if (ws.production_cmode_group_translated.has_value()) {
+            cmode = ws.production_cmode_group_translated.value();
+        } else {
+            return true; // can't know
+        }
+    }
+    // We only need to check modes that can result in small fractions
+    // i.e. ORAT, WRAT, GRAT, LRAT
+    const bool do_check = (cmode == Well::ProducerCMode::ORAT ||
+                           cmode == Well::ProducerCMode::WRAT ||
+                           cmode == Well::ProducerCMode::GRAT ||
+                           cmode == Well::ProducerCMode::LRAT);
+    if (!do_check) {
+        return true;
+    }
+    // somewhat ad-hoc tolerance
+    Scalar fraction_tolerance = 1e-6;
+    const auto& pu = well_.phaseUsage();
+    // Compute fraction corresponding to the control mode
+    Scalar cmode_frac = 0.0;
+    if (cmode == Well::ProducerCMode::ORAT || cmode == Well::ProducerCMode::LRAT) {
+        const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx);
+        cmode_frac += surface_fractions[oil_pos];
+    }
+    if (cmode == Well::ProducerCMode::WRAT || cmode == Well::ProducerCMode::LRAT) {
+        const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx);
+        cmode_frac += surface_fractions[water_pos];
+    }
+    if (cmode == Well::ProducerCMode::GRAT) {
+        const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx);
+        cmode_frac += surface_fractions[gas_pos];
+        fraction_tolerance *= 100;
+    }
+    // divide by sum just in case our fractions are not normalized
+    cmode_frac /= std::accumulate(surface_fractions.begin(), surface_fractions.end(), 0.0);
+
+    return cmode_frac > fraction_tolerance;
+}
+
+template<typename Scalar, typename IndexTraits>
+bool
+WellConstraints<Scalar, IndexTraits>::
+updateProducerControlMode(SingleWellState<Scalar, IndexTraits>& ws,
+                          const SummaryState& summaryState,
+                          const RateConvFunc& calcReservoirVoidageRates,
+                          const Well::ProductionControls& controls,
+                          const std::vector<Scalar>& surface_fractions,
+                          DeferredLogger& deferred_logger) const
+{
+    // We first estimate strictest individual/group rate constraints, and return true if this is violated. Otherwise,
+    // check if thp or bhp constraints are violated.
+    // NOTE: if a well is under group control, but it's current group target is stricter than previously, this function
+    // will return true (even though the mode does not change).
+    const auto currentControl = ws.production_cmode;
+    // use a relative (ad-hoc) tolerance
+    const Scalar tol = 1e-5;
+
+    // check most strict rate control.
+    const auto [most_strict_control, most_strict_scale] =
+        estimateStrictestProductionRateConstraint(ws, calcReservoirVoidageRates, controls,
+                                                  surface_fractions, /*skip_zero_rate_constraints*/ false);
+    if (most_strict_control != Well::ProducerCMode::CMODE_UNDEFINED && most_strict_scale < 1.0 - tol) {
+        ws.production_cmode = most_strict_control;
+        return true;
+    }
+    // check thp
+    if (well_.wellHasTHPConstraints(summaryState) && currentControl != Well::ProducerCMode::THP) {
+        const auto& thp = well_.getTHPConstraint(summaryState);
+        Scalar current_thp = ws.thp;
+        // For trivial group targets (for instance caused by NETV) we dont want to flip to THP control.
+        const bool dont_check = (currentControl == Well::ProducerCMode::GRUP && ws.trivial_group_target);
+        if (thp > current_thp*(1.0 + tol) && !dont_check) {
+            // If WVFPEXP item 4 is set to YES1 or YES2
+            // switching to THP is prevented if the well will
+            // produce at a higher rate with THP control
+            const auto& wvfpexp = well_.wellEcl().getWVFPEXP();
+            bool rate_less_than_potential = true;
+            if (wvfpexp.prevent()) {
+                for (int p = 0; p < well_.numPhases(); ++p) {
+                    // Currently we use the well potentials here computed before the iterations.
+                    // We may need to recompute the well potentials to get a more
+                    // accurate check here.
+                    rate_less_than_potential = rate_less_than_potential && (-ws.surface_rates[p]) <= ws.well_potentials[p];
+                }
+            }
+            if (!wvfpexp.prevent() || !rate_less_than_potential) {
+                //thp_limit_violated_but_not_switched = false;
+                ws.production_cmode = Well::ProducerCMode::THP;
+                return true;
+            } else {
+                //thp_limit_violated_but_not_switched = true;
+                deferred_logger.info("NOT_SWITCHING_TO_THP",
+                "The THP limit is violated for producer " +
+                well_.name() +
+                ". But the rate will increase if switched to THP. " +
+                "The well is therefore kept at " + WellProducerCMode2String(currentControl));
+            }
+        }
+    }
+    // check bhp
+    if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP) {
+        const Scalar bhp_limit = controls.bhp_limit;
+        Scalar current_bhp = ws.bhp;
+        if (bhp_limit > current_bhp*(1.0 + tol)) {
+            ws.production_cmode = Well::ProducerCMode::BHP;
+            return true;
+        }
+    }
+    // no constraints violated
+    return false;
+}
+
+template<typename Scalar, typename IndexTraits>
+std::pair<Well::ProducerCMode, Scalar>
+WellConstraints<Scalar, IndexTraits>::
+estimateStrictestProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                       const SummaryState& summaryState,
+                                       const RateConvFunc& calcReservoirVoidageRates,
+                                       const Well::ProductionControls& controls,
+                                       const std::vector<Scalar>& surface_fractions,
+                                       const bool skip_zero_rate_constraints,
+                                       DeferredLogger& deferred_logger,
+                                       const std::optional<Scalar> bhp_at_thp_limit) const
+{
+    // Estimate the most strict constraint + corresponding scaling based on current rate fractions:
+    // 1. If bhp_at_thp_limit not given: potential (if available) is used to approximate pressure constraint
+    //    rate-scaling - intended for initial guess
+    // 2. If bhp_at_thp_limit given: we assume a converged well-state with valid ipr - intended for
+    //    use within operability estimates
+
+    const auto rates = ws.surface_rates;
+    const auto tot_rates = std::accumulate(rates.begin(), rates.end(), 0.0);
+    if (std::abs(tot_rates) == 0.0) {
+        deferred_logger.debug("estimateStrictestProductionControl: current surface rates for well " +
+                              ws.name + " are zero. Cannot determine most strict control.");
+        return std::make_pair(Well::ProducerCMode::CMODE_UNDEFINED, 1.0);
+    }
+    Well::ProducerCMode most_strict_control = Well::ProducerCMode::CMODE_UNDEFINED;
+    Scalar most_strict_scale = std::numeric_limits<Scalar>::max();
+
+    if (!bhp_at_thp_limit.has_value() && false) {
+        const auto tot_potential = std::accumulate(ws.well_potentials.begin(), ws.well_potentials.end(), 0.0);
+        if (std::abs(tot_potential) > 0.0) {
+            most_strict_scale = -tot_potential/tot_rates;
+            if (well_.wellHasTHPConstraints(summaryState)) {
+                // not neccessarily true, but most likely
+                most_strict_control = Well::ProducerCMode::THP;
+            } else {
+                most_strict_control = Well::ProducerCMode::BHP;
+            }
+        }
+    } else if (bhp_at_thp_limit.has_value()){
+        if (*bhp_at_thp_limit > static_cast<Scalar>(controls.bhp_limit)) {
+            most_strict_control = Well::ProducerCMode::THP;
+        } else {
+            most_strict_control = Well::ProducerCMode::BHP;
+        }
+        const Scalar most_strict_bhp = std::max(*bhp_at_thp_limit, static_cast<Scalar>(controls.bhp_limit));
+        const Scalar tot_ipr_b = std::accumulate(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.0);
+        const Scalar tot_ipr_a = std::accumulate(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.0);
+        const Scalar tot_rate_at_bhp = tot_ipr_b*most_strict_bhp - tot_ipr_a;
+        deferred_logger.debug("estimateStrictestProductionControl: Well " + ws.name + 
+                              " has rate at bhp limit " + std::to_string(tot_rate_at_bhp) +
+                              " (bhp limit " + std::to_string(most_strict_bhp) + ")" + 
+                              " and current total rate " + std::to_string(tot_rates) + 
+                              " ipr_a " + std::to_string(tot_ipr_a) +
+                              " ipr_b " + std::to_string(tot_ipr_b));
+        most_strict_scale = tot_rate_at_bhp/tot_rates;
+    }
+    // check rate constraints
+    const auto [most_strict_rate_control, most_strict_rate_scale] =
+        estimateStrictestProductionRateConstraint(ws, calcReservoirVoidageRates, controls, surface_fractions, skip_zero_rate_constraints);
+    if (most_strict_rate_control != Well::ProducerCMode::CMODE_UNDEFINED && most_strict_rate_scale < most_strict_scale) {
+        most_strict_scale = most_strict_rate_scale;
+        most_strict_control = most_strict_rate_control;
+    }
+    // If we are computing e.g., well potentials, we may still have CMODE_UNDEFINED at this point. In that case, we scale according to the 
+    // previous rate and keep the current control mode.
+    if (most_strict_control == Well::ProducerCMode::CMODE_UNDEFINED) {
+        const auto tot_prev_rates = std::accumulate(ws.prev_surface_rates.begin(), ws.prev_surface_rates.end(), 0.0);
+        if (std::abs(tot_prev_rates) > 0.0) {
+            most_strict_scale = std::abs(tot_prev_rates/tot_rates);
+            most_strict_control = ws.production_cmode;
+        } else {
+            deferred_logger.debug("estimateStrictestProductionControl: previous surface rates for well " +
+                                  ws.name + " are zero. This is a BUG!.");
+            deferred_logger.debug("  Previous rates are " + std::to_string(ws.prev_surface_rates[0]) + ", " + std::to_string(ws.prev_surface_rates[1]) + ", " + std::to_string(ws.prev_surface_rates[2]) + ". ");
+        }
+    }
+    return std::make_pair(most_strict_control, most_strict_scale);
+}
+
+template<typename Scalar, typename IndexTraits>
+std::pair<Well::ProducerCMode, Scalar>
+WellConstraints<Scalar, IndexTraits>::
+estimateStrictestProductionRateConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                          const RateConvFunc& calcReservoirVoidageRates,
+                                          const Well::ProductionControls& controls,
+                                          const std::vector<Scalar>& surface_fractions,
+                                          const bool skip_zero_rate_constraints) const
+{
+    // Estimate the most strict rate constraint + corresponding scaling based on current rate fractions
+    const auto rates = ws.surface_rates;
+    // If all rates are zero, we can estimate strictest mode from fractions, but no scale
+    const Scalar zero_rates = std::accumulate(rates.begin(), rates.end(), 0.0) == 0.0;
+
+    Well::ProducerCMode most_strict_control = Well::ProducerCMode::CMODE_UNDEFINED;
+    Scalar most_strict_scale = std::numeric_limits<Scalar>::max();
+
+    // start with individual rate constraints
+    const std::array<Well::ProducerCMode, 5> rate_modes = {Well::ProducerCMode::ORAT,
+                                                          Well::ProducerCMode::WRAT,
+                                                          Well::ProducerCMode::GRAT,
+                                                          Well::ProducerCMode::LRAT,
+                                                          Well::ProducerCMode::RESV};
+    for (const auto& mode : rate_modes) {
+        if (!controls.hasControl(mode))
+            continue;
+        const Scalar scale = getProductionControlModeScale(ws, calcReservoirVoidageRates, mode, controls, surface_fractions, skip_zero_rate_constraints);
+        if (scale >= 0.0 && scale < most_strict_scale) {
+            most_strict_scale = scale;
+            most_strict_control = mode;
+        }
+    }
+    // check group constraints if target is given in well-state
+    if (controls.hasControl(Well::ProducerCMode::GRUP) && !ws.prevent_group_control && ws.group_target.has_value() && ws.production_cmode_group_translated.has_value()) {
+        const Scalar scale = getProductionControlModeScale(ws, calcReservoirVoidageRates, ws.production_cmode_group_translated.value(), controls,
+                                                           surface_fractions, skip_zero_rate_constraints, ws.group_target.value());
+        if (scale >= 0.0 && scale < most_strict_scale) {
+            most_strict_scale = scale;
+            most_strict_control = Well::ProducerCMode::GRUP;
+        }
+    }
+    if (zero_rates) {
+        most_strict_scale = std::numeric_limits<Scalar>::max();
+    }
+    return std::make_pair(most_strict_control, most_strict_scale);
+}
+
+template<typename Scalar, typename IndexTraits>
+Scalar
+WellConstraints<Scalar, IndexTraits>::
+getProductionControlModeScale(const SingleWellState<Scalar, IndexTraits>& ws,
+                              const RateConvFunc& calcReservoirVoidageRates,
+                              const Well::ProducerCMode& cmode,
+                              const Well::ProductionControls& controls,
+                              const std::vector<Scalar>& surface_fractions,
+                              const bool skip_zero_rate_constraints,
+                              const std::optional<Scalar> target) const
+{
+    // check if cmode is feasible with current rate fractions
+    if (!isFeasibleProductionConstraint(ws, surface_fractions, cmode)) {
+        return -1.0;
+    }
+    const auto& pu = well_.phaseUsage();
+    Scalar current_rate = 0.0;
+    Scalar target_rate = 0.0;
+    // use fractions if all rates are zero 
+    const bool zero_rates = std::accumulate(ws.surface_rates.begin(), ws.surface_rates.end(), 0.0) == 0.0;
+    const auto rates = zero_rates ? surface_fractions : ws.surface_rates;
+    switch (cmode) {
+        case Well::ProducerCMode::ORAT:
+            current_rate = -rates[pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx)];
+            target_rate = target.has_value() ? *target : controls.oil_rate;
+            break;
+        case Well::ProducerCMode::WRAT:
+            current_rate = -rates[pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx)];
+            target_rate = target.has_value() ? *target : controls.water_rate;
+            break;
+        case Well::ProducerCMode::GRAT:
+            current_rate = -rates[pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx)];
+            target_rate = target.has_value() ? *target : controls.gas_rate;
+            break;
+        case Well::ProducerCMode::LRAT:
+            current_rate = -rates[pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx)];
+            current_rate += -rates[pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx)];
+            target_rate = target.has_value() ? *target : controls.liquid_rate;
+            break;
+        case Well::ProducerCMode::RESV:
+            for (int p = 0; p < well_.numPhases(); ++p) {
+                current_rate += -ws.reservoir_rates[p];
+            }
+            if (controls.prediction_mode || target.has_value()) {
+                target_rate = target.has_value() ? *target : controls.resv_rate;
+            } else {
+                const int fipreg = 0; // not considering the region for now
+                const int np = well_.numPhases();
+
+                std::vector<Scalar> target_surface_rates(np, 0.0);
+                if (pu.phaseIsActive(IndexTraits::waterPhaseIdx))
+                    target_surface_rates[pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx)] = controls.water_rate;
+                if (pu.phaseIsActive(IndexTraits::oilPhaseIdx))
+                    target_surface_rates[pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx)] = controls.oil_rate;
+                if (pu.phaseIsActive(IndexTraits::gasPhaseIdx))
+                    target_surface_rates[pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx)] = controls.gas_rate;
+
+                std::vector<Scalar> target_voidage_rates(np, 0.0);
+                calcReservoirVoidageRates(fipreg, well_.pvtRegionIdx(), target_surface_rates, target_voidage_rates);
+                // Scalar resv_rate = 0.0;
+                for (int p = 0; p < np; ++p) {
+                    target_rate += target_voidage_rates[p];
+                }
+            }
+            break;
+        default:
+            // undefined mode, no scaling applied
+            return -1.0;
+            break;
+    }
+    const bool valid_scale = (!skip_zero_rate_constraints || std::abs(target_rate) > 0.0);
+    return  valid_scale ? std::abs(target_rate/current_rate) : -1.0;
+}
+
 template class WellConstraints<double, BlackOilDefaultFluidSystemIndices>;
 
 #if FLOW_INSTANTIATE_FLOAT
diff --git a/opm/simulators/wells/WellConstraints.hpp b/opm/simulators/wells/WellConstraints.hpp
index d4f062513e6..42707a015ff 100644
--- a/opm/simulators/wells/WellConstraints.hpp
+++ b/opm/simulators/wells/WellConstraints.hpp
@@ -63,6 +63,36 @@ class WellConstraints {
                                const std::optional<Well::InjectionControls>& inj_controls = std::nullopt,
                                const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
 
+    bool
+    isFeasibleProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                   const std::vector<Scalar>& surface_fractions,
+                                   const Well::ProducerCMode cmode) const;
+
+    bool
+    updateProducerControlMode(SingleWellState<Scalar, IndexTraits>& ws,
+                              const SummaryState& summaryState,
+                              const RateConvFunc& calcReservoirVoidageRates,
+                              const Well::ProductionControls& controls,
+                              const std::vector<Scalar>& surface_fractions,
+                              DeferredLogger& deferred_logger) const;
+
+    std::pair<WellProducerCMode, Scalar>
+    estimateStrictestProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                          const SummaryState& summaryState,
+                                          const RateConvFunc& calcReservoirVoidageRates,
+                                          const WellProductionControls& controls,
+                                          const std::vector<Scalar>& surface_fractions,
+                                          const bool skip_zero_rate_constraints,
+                                          DeferredLogger& deferred_logger,
+                                          const std::optional<Scalar> bhp_at_thp_limit = std::nullopt) const;
+
+    std::pair<WellProducerCMode, Scalar>
+    estimateStrictestProductionRateConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                          const RateConvFunc& calcReservoirVoidageRates,
+                                          const Well::ProductionControls& controls,
+                                          const std::vector<Scalar>& surface_fractions,
+                                          const bool skip_zero_rate_constraints) const;
+
 private:
     WellInjectorCMode
     activeInjectionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
@@ -79,6 +109,15 @@ class WellConstraints {
                                DeferredLogger& deferred_logger,
                                const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
 
+    Scalar
+    getProductionControlModeScale(const SingleWellState<Scalar, IndexTraits>& ws,
+                                  const RateConvFunc& calcReservoirVoidageRates,
+                                  const WellProducerCMode& cmode,
+                                  const WellProductionControls& control,
+                                  const std::vector<Scalar>& surface_fractions,
+                                  const bool skip_zero_rate_constraints,
+                                  const std::optional<Scalar> target = std::nullopt) const;
+
     const WellInterfaceGeneric<Scalar, IndexTraits>& well_; //!< Reference to well interface
 };
 
diff --git a/opm/simulators/wells/WellGroupConstraints.cpp b/opm/simulators/wells/WellGroupConstraints.cpp
index c0cc0dd160e..c79612cf65e 100644
--- a/opm/simulators/wells/WellGroupConstraints.cpp
+++ b/opm/simulators/wells/WellGroupConstraints.cpp
@@ -131,7 +131,9 @@ checkGroupConstraints(const WellGroupHelperType& wgHelper,
     const auto& well = well_.wellEcl();
     const int well_index = well_.indexOfWell();
     auto& ws = well_state.well(well_index);
-
+    if (ws.prevent_group_control) {
+        return false;
+    }
     if (well.isInjector()) {
         const auto currentControl = ws.injection_cmode;
 
diff --git a/opm/simulators/wells/WellGroupHelper.cpp b/opm/simulators/wells/WellGroupHelper.cpp
index 0fcbcfc8326..a2494b5d9ab 100644
--- a/opm/simulators/wells/WellGroupHelper.cpp
+++ b/opm/simulators/wells/WellGroupHelper.cpp
@@ -641,7 +641,7 @@ WellGroupHelper<Scalar, IndexTraits>::getProductionGroupRateVector(const std::st
 }
 
 template <typename Scalar, typename IndexTraits>
-Scalar
+std::pair<WellInjectorCMode, Scalar>
 WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::string& name,
                                                                  const std::string& parent,
                                                                  const Group& group,
@@ -663,7 +663,7 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::stri
         || current_group_control == Group::InjectionCMode::NONE) {
         // Return if we are not available for parent group.
         if (!group.injectionGroupControlAvailable(injection_phase)) {
-            return std::numeric_limits<Scalar>::max();
+            return std::make_pair(WellInjectorCMode::CMODE_UNDEFINED, std::numeric_limits<Scalar>::max());
         }
         // Otherwise: check production share of parent's control.
         const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_);
@@ -680,7 +680,7 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::stri
     // This can be false for FLD-controlled groups, we must therefore
     // check for FLD first (done above).
     if (!group.isInjectionGroup()) {
-        return std::numeric_limits<Scalar>::max();
+        return std::make_pair(WellInjectorCMode::CMODE_UNDEFINED, std::numeric_limits<Scalar>::max());
     }
 
     // If we are here, we are at the topmost group to be visited in the recursion.
@@ -771,11 +771,14 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetInjector(const std::stri
         }
     }
     // Avoid negative target rates comming from too large local reductions.
-    return std::max(Scalar(0.0), target / efficiency_factor);
+    target = std::max(Scalar(0.0), target)/efficiency_factor;
+    // Translating injector group control to individual well control requires some care.
+    // We just set it to undefined for now, and update at a later stage if needed.
+    return std::make_pair(WellInjectorCMode::CMODE_UNDEFINED, target);
 }
 
 template <typename Scalar, typename IndexTraits>
-Scalar
+std::pair<WellProducerCMode, Scalar>
 WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::string& name,
                                                                  const std::string& parent,
                                                                  const Group& group,
@@ -798,7 +801,7 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::stri
         || current_group_control == Group::ProductionCMode::NONE) {
         // Return if we are not available for parent group.
         if (!group.productionGroupControlAvailable()) {
-            return std::numeric_limits<Scalar>::max();
+            return std::make_pair(WellProducerCMode::CMODE_UNDEFINED, std::numeric_limits<Scalar>::max());
         }
         // Otherwise: check production share of parent's control.
         const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_);
@@ -814,7 +817,7 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::stri
     // This can be false for FLD-controlled groups, we must therefore
     // check for FLD first (done above).
     if (!group.isProductionGroup()) {
-        return std::numeric_limits<Scalar>::max();
+        return std::make_pair(WellProducerCMode::CMODE_UNDEFINED, std::numeric_limits<Scalar>::max());
     }
 
     // If we are here, we are at the topmost group to be visited in the recursion.
@@ -899,8 +902,30 @@ WellGroupHelper<Scalar, IndexTraits>::getWellGroupTargetProducer(const std::stri
         }
     }
     // Avoid negative target rates comming from too large local reductions.
-    return std::max(Scalar(0.0), target / efficiency_factor);
-    ;
+    target = std::max(Scalar(0.0), target)/efficiency_factor;
+    // Translate group control to well control
+    WellProducerCMode mode;
+    switch (current_group_control) {
+        case Group::ProductionCMode::ORAT:
+            mode = WellProducerCMode::ORAT;
+            break;
+        case Group::ProductionCMode::GRAT:
+            mode = WellProducerCMode::GRAT;
+            break;
+        case Group::ProductionCMode::WRAT:
+            mode = WellProducerCMode::WRAT;
+            break;
+        case Group::ProductionCMode::LRAT:
+            mode = WellProducerCMode::LRAT;
+            break;
+        case Group::ProductionCMode::RESV:
+            mode = WellProducerCMode::RESV;
+            break;
+        default:
+            mode = WellProducerCMode::CMODE_UNDEFINED;
+            break;
+    }
+    return std::make_pair(mode, target);
 }
 
 template <typename Scalar, typename IndexTraits>
diff --git a/opm/simulators/wells/WellGroupHelper.hpp b/opm/simulators/wells/WellGroupHelper.hpp
index 5437dca244f..29fbf0d95af 100644
--- a/opm/simulators/wells/WellGroupHelper.hpp
+++ b/opm/simulators/wells/WellGroupHelper.hpp
@@ -144,22 +144,24 @@ class WellGroupHelper
 
     GuideRate::RateVector getProductionGroupRateVector(const std::string& group_name) const;
 
-    Scalar getWellGroupTargetInjector(const std::string& name,
-                                      const std::string& parent,
-                                      const Group& group,
-                                      const Scalar* rates,
-                                      const Phase injection_phase,
-                                      const Scalar efficiency_factor,
-                                      const std::vector<Scalar>& resv_coeff,
-                                      DeferredLogger& deferred_logger) const;
-
-    Scalar getWellGroupTargetProducer(const std::string& name,
-                                      const std::string& parent,
-                                      const Group& group,
-                                      const Scalar* rates,
-                                      const Scalar efficiency_factor,
-                                      const std::vector<Scalar>& resv_coeff,
-                                      DeferredLogger& deferred_logger) const;
+    std::pair<WellInjectorCMode, Scalar>
+    getWellGroupTargetInjector(const std::string& name,
+                               const std::string& parent,
+                               const Group& group,
+                               const Scalar* rates,
+                               const Phase injection_phase,
+                               const Scalar efficiency_factor,
+                               const std::vector<Scalar>& resv_coeff,
+                               DeferredLogger& deferred_logger) const;
+
+    std::pair<WellProducerCMode, Scalar>
+    getWellGroupTargetProducer(const std::string& name,
+                               const std::string& parent,
+                               const Group& group,
+                               const Scalar* rates,
+                               const Scalar efficiency_factor,
+                               const std::vector<Scalar>& resv_coeff,
+                               DeferredLogger& deferred_logger) const;
 
     GuideRate::RateVector getWellRateVector(const std::string& name) const;
 
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 5798b3cd0ee..f0df426c3eb 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -263,6 +263,7 @@ class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Propertie
                                                   const Well::ProductionControls& prod_controls,
                                                   const Scalar WQTotal,
                                                   WellStateType& well_state,
+                                                  const std::vector<Scalar>& surface_fractions,
                                                   DeferredLogger& deferred_logger,
                                                   const bool fixed_control = false,
                                                   const bool fixed_status = false);
diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
index 6882d6b61ad..c722ab01f57 100644
--- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp
+++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
@@ -331,6 +331,84 @@ zeroGroupRateTarget(const SummaryState& summary_state,
     }
 }
 
+template<typename FluidSystem>
+bool
+WellInterfaceFluidSystem<FluidSystem>::
+isFeasibleProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                               const std::vector<Scalar>& surface_fractions,
+                               const Well::ProducerCMode cmode) const
+{
+    return WellConstraints(*this).
+        isFeasibleProductionConstraint(ws, surface_fractions, cmode);
+}
+
+template<typename FluidSystem>
+bool
+WellInterfaceFluidSystem<FluidSystem>::
+updateProducerControlMode(SingleWellState<Scalar, IndexTraits>& ws,
+                          const SummaryState& summaryState,
+                          const Well::ProductionControls& controls,
+                          const std::vector<Scalar>& surface_fractions,
+                          DeferredLogger& deferred_logger) const
+{   
+    auto rRates = [this](const int fipreg,
+                         const int pvtRegion,
+                         const std::vector<Scalar>& surface_rates,
+                         std::vector<Scalar>& voidage_rates)
+    {
+        return rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegion,
+                                                        surface_rates, voidage_rates);
+    };
+    return WellConstraints(*this).
+            updateProducerControlMode(ws, summaryState, rRates, controls, surface_fractions, deferred_logger);
+}
+
+
+template<typename FluidSystem>
+std::pair<Well::ProducerCMode, typename FluidSystem::Scalar>
+WellInterfaceFluidSystem<FluidSystem>::
+estimateStrictestProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                    const SummaryState& summaryState,
+                                    const Well::ProductionControls& controls,
+                                    const std::vector<Scalar>& surface_fractions,
+                                    const bool skip_zero_rate_constraints,
+                                    DeferredLogger& deferred_logger,
+                                    const std::optional<Scalar> bhp_at_thp_limit) const
+{
+    auto rRates = [this](const int fipreg,
+                         const int pvtRegion,
+                         const std::vector<Scalar>& surface_rates,
+                         std::vector<Scalar>& voidage_rates)
+    {
+        return rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegion,
+                                                        surface_rates, voidage_rates);
+    };
+    return WellConstraints(*this).
+            estimateStrictestProductionConstraint(ws, summaryState, rRates, controls, surface_fractions,
+                                                  skip_zero_rate_constraints, deferred_logger, bhp_at_thp_limit);
+}
+
+template<typename FluidSystem>
+std::pair<Well::ProducerCMode, typename FluidSystem::Scalar>
+WellInterfaceFluidSystem<FluidSystem>::
+estimateStrictestProductionRateConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                          const Well::ProductionControls& controls,
+                                          const std::vector<Scalar>& surface_fractions,
+                                          const bool skip_zero_rate_constraints) const
+{
+    auto rRates = [this](const int fipreg,
+                         const int pvtRegion,
+                         const std::vector<Scalar>& surface_rates,
+                         std::vector<Scalar>& voidage_rates)
+    {
+        return rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegion,
+                                                        surface_rates, voidage_rates);
+    };
+    return WellConstraints(*this).
+            estimateStrictestProductionRateConstraint(ws, rRates, controls, surface_fractions,
+                                                      skip_zero_rate_constraints);
+}
+
 template<class Scalar>
 using FS = BlackOilFluidSystem<Scalar, BlackOilDefaultFluidSystemIndices>;
 
diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.hpp b/opm/simulators/wells/WellInterfaceFluidSystem.hpp
index 2f4654a95c9..9b91b441020 100644
--- a/opm/simulators/wells/WellInterfaceFluidSystem.hpp
+++ b/opm/simulators/wells/WellInterfaceFluidSystem.hpp
@@ -129,6 +129,31 @@ class WellInterfaceFluidSystem : public WellInterfaceGeneric<typename FluidSyste
                              const GroupState<Scalar>& group_state,
                              DeferredLogger& deferredLogger) const;
 
+    bool isFeasibleProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                        const std::vector<Scalar>& surface_fractions,
+                                        const Well::ProducerCMode cmode) const;
+
+    bool updateProducerControlMode(SingleWellState<Scalar, IndexTraits>& ws,
+                                   const SummaryState& summaryState,
+                                   const Well::ProductionControls& controls,
+                                   const std::vector<Scalar>& surface_fractions,
+                                   DeferredLogger& deferred_logger) const;
+
+    std::pair<Well::ProducerCMode, Scalar>
+    estimateStrictestProductionConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                          const SummaryState& summaryState,
+                                          const Well::ProductionControls& controls,
+                                          const std::vector<Scalar>& surface_fractions,
+                                          const bool skip_zero_rate_constraints,
+                                          DeferredLogger& deferred_logger,
+                                          const std::optional<Scalar> bhp_at_thp_limit = std::nullopt) const;
+
+    std::pair<Well::ProducerCMode, Scalar>
+    estimateStrictestProductionRateConstraint(const SingleWellState<Scalar, IndexTraits>& ws,
+                                              const Well::ProductionControls& controls,
+                                              const std::vector<Scalar>& surface_fractions,
+                                              const bool skip_zero_rate_constraintsr) const;
+
     // For the conversion between the surface volume rate and reservoir voidage rate
     const RateConverterType& rateConverter_;
 };
diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp
index 8e678f9e732..3446f335a4f 100644
--- a/opm/simulators/wells/WellInterfaceGeneric.cpp
+++ b/opm/simulators/wells/WellInterfaceGeneric.cpp
@@ -450,7 +450,66 @@ setPrevSurfaceRates(WellState<Scalar, IndexTraits>& well_state,
         ws.prev_surface_rates = ws.surface_rates;
     }
 }
+/*
+template<typename Scalar, typename IndexTraits>
+bool WellInterfaceGeneric<Scalar, IndexTraits>::isFeasibleProductionControl(const  WellState<Scalar, IndexTraits>& well_state,
+                                                                            const Well::ProductionControls& prod_controls,
+                                                                            const std::vector<Scalar>& fractions,
+                                                                            const std::vector<Scalar>& fraction_scaling,
+                                                                            std::optional<WellProducerCMode> cmode_opt) const
+{
+    // Avoid problematic control cases like e.g., WRAT for well with approx zero
+    // water fraction.
+    assert(this->isProducer());
+    auto& ws = well_state.well(this->index_of_well_);
+    const WellProducerCMode cmode = cmode_opt.has_value() ? cmode_opt.value() : ws.production_cmode;
+
+    if (!prod_controls.hasControl(cmode)) {
+        return false; 
+    }
+    if (cmode == Well::ProducerCMode::GRUP) {
+        if (ws.production_cmode_group_translated.has_value()) {
+            cmode = ws.production_cmode_group_translated.value();
+        } else {
+            return true; // can't know
+        }
+    } else {
+        cmode = ws.production_cmode;
+    }
+     // We only need to check ORAT, WRAT, GRAT, LRAT
+    const bool do_check = (gcmode == Well::ProducerCMode::ORAT ||
+                           gcmode == Well::ProducerCMode::WRAT ||
+                           gcmode == Well::ProducerCMode::GRAT ||
+                           gcmode == Well::ProducerCMode::LRAT);
+    if (!do_check) {
+        return true;
+    }
+    const auto& pu = this->phaseUsage();
+    auto scaled_fractions = fractions;
+    for (size_t i = 0; i < scaled_fractions.size(); ++i) {
+        scaled_fractions[i] *= fraction_scaling[i];
+    }
+    const auto scaled_frac_sum = std::accumulate(scaled_fractions.begin(), scaled_fractions.end(), 0.0);
+    // Compute weighted fraction corresponding to the control mode
+    Scalar scaled_cmode_frac = 0.0;
+    if (cmode == Well::ProducerCMode::ORAT || cmode == Well::ProducerCMode::LRAT) {
+        const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx);
+        scaled_cmode_frac += weighted_fractions[oil_pos];
+    }
+    if (cmode == Well::ProducerCMode::WRAT || cmode == Well::ProducerCMode::LRAT) {
+        const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx);
+        scaled_cmode_frac += weighted_fractions[water_pos];
+    }
+    if (cmode == Well::ProducerCMode::GRAT) {
+        const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx);
+        scaled_cmode_frac += weighted_fractions[gas_pos];
+    }
+    const Scalar fraction_cur = scaled_cmode_frac / scaled_frac_sum;
+    const Scalar fraction_tol = this->param_.tolerance_wells_;
 
+    return fraction_cur > fraction_tol; 
+}
+*/
 template<typename Scalar, typename IndexTraits>
 void WellInterfaceGeneric<Scalar, IndexTraits>::
 setGuideRate(const GuideRate* guide_rate_arg)
diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp
index e500ccbe0f1..6b11dfb278a 100644
--- a/opm/simulators/wells/WellInterfaceGeneric.hpp
+++ b/opm/simulators/wells/WellInterfaceGeneric.hpp
@@ -107,6 +107,11 @@ class WellInterfaceGeneric {
     void setVFPProperties(const VFPProperties<Scalar, IndexTraits>* vfp_properties_arg);
     void setPrevSurfaceRates(WellStateType& well_state,
                              const WellStateType& prev_well_state) const;
+    bool isFeasibleProductionControl(const  WellStateType& well_state,
+                                     const Well::ProductionControls& prod_controls,
+                                     const std::vector<Scalar>& fractions,
+                                     const std::vector<Scalar>& fraction_scaling,
+                                     std::optional<WellProducerCMode> cmode_opt) const;
     void setGuideRate(const GuideRate* guide_rate_arg);
     void setWellEfficiencyFactor(const Scalar efficiency_factor);
     void setRepRadiusPerfLength();
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 598bdf60d60..ded13fb9d8d 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -280,6 +280,7 @@ namespace Opm
                                              const Well::ProductionControls& prod_controls,
                                              const Scalar wqTotal,
                                              WellStateType& well_state,
+                                             const std::vector<Scalar>& surface_fractions,
                                              DeferredLogger& deferred_logger,
                                              const bool fixed_control,
                                              const bool fixed_status)
@@ -296,7 +297,7 @@ namespace Opm
         }
         const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
 
-        if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(well_state.well(this->index_of_well_).status == WellStatus::OPEN)) {
+        if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(ws.status == WellStatus::OPEN)) {
            return false;
         }
 
@@ -311,24 +312,27 @@ namespace Opm
                     // Changing to group controls here may lead to inconsistencies in the group handling which in turn
                     // may result in excessive back and forth switching. However, we currently allow this by default.
                     // The switch check_group_constraints_inner_well_iterations_ is a temporary solution.
-
                     const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) :
                                                                       prod_controls.hasControl(Well::ProducerCMode::GRUP);
-                    bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
-                    if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
-                        changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
-                    }
-                    if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
-                        changed = changed || this->checkGroupConstraints(
-                            wgHelper, schedule, summary_state, false, well_state, deferred_logger
-                        );
+                    const bool hasGroupProdTarget = ws.group_target.has_value() && ws.production_cmode_group_translated.has_value();
+                    if (this->isProducer() && (!hasGroupControl || hasGroupProdTarget)) {
+                        changed = this->updateProducerControlMode(ws, summary_state, prod_controls, surface_fractions, deferred_logger);
+                    } else {
+                        bool isGroupControl = ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP;
+                        if (! (isGroupControl && !this->param_.check_group_constraints_inner_well_iterations_)) {
+                            changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
+                        }
+                        if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
+                            changed = changed || this->checkGroupConstraints(
+                                wgHelper, schedule, summary_state, false, well_state, deferred_logger
+                            );
+                        }
                     }
-
                     if (changed) {
                         const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP :
                                                                         ws.production_cmode == Well::ProducerCMode::THP;
                         if (thp_controlled){
-                             ws.thp = this->getTHPConstraint(summary_state);
+                            ws.thp = this->getTHPConstraint(summary_state);
                         } else {
                             // don't call for thp since this might trigger additional local solve
                             updateWellStateWithTarget(simulator, wgHelper, well_state, deferred_logger);
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 703c878a0e7..f6eddab58ed 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -416,6 +416,8 @@ init(const std::vector<Scalar>& cellPressures,
             new_well.reservoir_rates = prev_well.reservoir_rates;
             new_well.well_potentials = prev_well.well_potentials;
             new_well.group_target = prev_well.group_target;
+            new_well.production_cmode_group_translated = prev_well.production_cmode_group_translated;
+            new_well.injection_cmode_group_translated = prev_well.injection_cmode_group_translated;
 
             // perfPhaseRates
             //