diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp
index 390c64c16..f0659a23b 100644
--- a/opm/simulators/wells/MultisegmentWell.hpp
+++ b/opm/simulators/wells/MultisegmentWell.hpp
@@ -258,6 +258,10 @@ namespace Opm
                                  const Simulator& ebos_simulator,
                                  DeferredLogger& deferred_logger) const;
 
+        bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                              WellState& well_state,
+                                              DeferredLogger& deferred_logger) const;
+
         virtual double getRefDensity() const override;
 
         virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 21a7cbb48..3fd7c5180 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -2020,4 +2020,32 @@ namespace Opm
         connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
     }
 
+
+
+
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                     WellState& well_state,
+                                     DeferredLogger& deferred_logger) const
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+
+        auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
+            ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
+        if (bhp_at_thp_limit) {
+            std::vector<double> rates(this->number_of_phases_, 0.0);
+            computeWellRatesWithBhpIterations(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
+            auto& ws = well_state.well(this->name());
+            ws.surface_rates = rates;
+            ws.bhp = *bhp_at_thp_limit;
+            ws.thp = this->getTHPConstraint(summary_state);
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+
 } // namespace Opm
diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp
index 720be9a2b..e819fcba0 100644
--- a/opm/simulators/wells/StandardWell.hpp
+++ b/opm/simulators/wells/StandardWell.hpp
@@ -334,6 +334,9 @@ namespace Opm
             DeferredLogger& deferred_logger,
             const WellState &well_state) const;
 
+        bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                              WellState& well_state,
+                                              DeferredLogger& deferred_logger) const;
 
         virtual double getRefDensity() const override;
 
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 2b322de00..0131bc75a 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -1783,6 +1783,34 @@ namespace Opm
         return potentials;
     }
 
+
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                     WellState& well_state,
+                                     DeferredLogger& deferred_logger) const
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+
+        auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
+            ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
+        if (bhp_at_thp_limit) {
+            std::vector<double> rates(this->number_of_phases_, 0.0);
+            computeWellRatesWithBhp(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
+            auto& ws = well_state.well(this->name());
+            ws.surface_rates = rates;
+            ws.bhp = *bhp_at_thp_limit;
+            ws.thp = this->getTHPConstraint(summary_state);
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+
+
     template<typename TypeTag>
     double
     StandardWell<TypeTag>::
@@ -2363,9 +2391,14 @@ namespace Opm
                                                                                alq_value,
                                                                                this->getTHPConstraint(summary_state),
                                                                                deferred_logger);
-        auto v = frates(*bhpAtLimit);
-        if (bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
-            return bhpAtLimit;
+
+        if (bhpAtLimit) {
+            auto v = frates(*bhpAtLimit);
+            if (std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }) ) {
+                return bhpAtLimit;
+            }
+        }
+
 
         auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
             // Solver the well iterations to see if we are
@@ -2384,9 +2417,15 @@ namespace Opm
                                                                           alq_value,
                                                                           this->getTHPConstraint(summary_state),
                                                                           deferred_logger);
-        v = frates(*bhpAtLimit);
-        if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
-            return bhpAtLimit;
+
+
+        if (bhpAtLimit) {
+            // should we use fratesIter here since fratesIter is used in computeBhpAtThpLimitProd above?
+            auto v = frates(*bhpAtLimit);
+            if (std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }) ) {
+                return bhpAtLimit;
+            }
+        }
 
         // we still don't get a valied solution.
         return std::nullopt;
diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp
index fd0170509..00a6c5a87 100644
--- a/opm/simulators/wells/WellBhpThpCalculator.cpp
+++ b/opm/simulators/wells/WellBhpThpCalculator.cpp
@@ -738,7 +738,7 @@ bruteForceBracket(const std::function<double(const double)>& eq,
     bool finding_bracket = false;
     low = range[0];
     high = range[1];
-    const int sample_number = 100;
+    const int sample_number = 200;
     const double interval = (high - low) / sample_number;
     double eq_low = eq(low);
     double eq_high;
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 341cb52a8..b6c57984d 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -201,6 +201,10 @@ public:
                                            WellState& well_state,
                                            DeferredLogger& deferred_logger) const;
 
+    virtual bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                                  WellState& well_state,
+                                                  DeferredLogger& deferred_logger) const = 0;
+
     enum class IndividualOrGroup { Individual, Group, Both };
     bool updateWellControl(const Simulator& ebos_simulator,
                            const IndividualOrGroup iog,
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 9e94a2ec8..695dfb2e9 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -965,24 +965,26 @@ namespace Opm
             }
             case Well::ProducerCMode::THP:
             {
-                auto rates = ws.surface_rates;
-                this->adaptRatesForVFP(rates);
-                double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
-                                                                             rates,
-                                                                             well,
-                                                                             summaryState,
-                                                                             this->getRefDensity(),
-                                                                             deferred_logger);
-                ws.bhp = bhp;
-                ws.thp = this->getTHPConstraint(summaryState);
+                const bool update_success = updateWellStateWithTHPTargetProd(ebos_simulator, well_state, deferred_logger);
 
-                // if the total rates are negative or zero
-                // we try to provide a better intial well rate
-                // using the well potentials
-                double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
-                if (total_rate <= 0.0){
-                    for (int p = 0; p<np; ++p) {
-                        ws.surface_rates[p] = -ws.well_potentials[p];
+                if (!update_success) {
+                    // the following is the original way of initializing well state with THP constraint
+                    // keeping it for robust reason in case that it fails to get a bhp value with THP constraint
+                    // more sophisticated design might be needed in the future
+                    auto rates = ws.surface_rates;
+                    this->adaptRatesForVFP(rates);
+                    const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
+                        well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
+                    ws.bhp = bhp;
+                    ws.thp = this->getTHPConstraint(summaryState);
+                    // if the total rates are negative or zero
+                    // we try to provide a better initial well rate
+                    // using the well potentials
+                    const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
+                    if (total_rate <= 0.0) {
+                        for (int p = 0; p < this->number_of_phases_; ++p) {
+                            ws.surface_rates[p] = -ws.well_potentials[p];
+                        }
                     }
                 }
                 break;