diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 8f20f07d1..aaa414c8d 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -106,4 +106,5 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp
+	opm/autodiff/WellStateFullyImplicitBlackoil.hpp
 	)
diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp
index cc305a78d..2456f9c04 100644
--- a/examples/sim_fibo_ad.cpp
+++ b/examples/sim_fibo_ad.cpp
@@ -39,7 +39,7 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
@@ -200,46 +200,66 @@ try
         param.writeParam(output_dir + "/simulation.param");
     }
 
-#warning "TODO: convert the well handling code to the new parser"
 #if USE_NEW_PARSER
     std::cout << "\n\n================    Starting main simulation loop     ===============\n"
               << std::flush;
 
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
     SimulatorTimer simtimer;
+    std::shared_ptr<EclipseState> eclipseState(new EclipseState(newParserDeck));
 
-    // Create new wells, well_state
-    WellsManager wells(*deck, *grid->c_grid(), props->permeability());
-    // @@@ HACK: we should really make a new well state and
-    // properly transfer old well state to it every epoch,
-    // since number of wells may change etc.
-    well_state.init(wells.c_wells(), state);
+    // initialize variables
+    simtimer.init(timeMap, /*beginReportStepIdx=*/0, /*endReportStepIdx=*/0);
 
-    // Create and run simulator.
-    SimulatorFullyImplicitBlackoil simulator(param,
-                                             *grid->c_grid(),
-                                             *new_props,
-                                             rock_comp->isActive() ? rock_comp.get() : 0,
-                                             wells,
-                                             linsolver,
-                                             grav,
-                                             outputWriter);
-    simtimer.init(timeMap);
-    SimulatorReport report = simulator.run(simtimer, state, well_state);
+    SimulatorReport fullReport;
+    for (size_t episodeIdx = 0; episodeIdx < timeMap->numTimesteps(); ++episodeIdx) {
+        WellsManager wells(eclipseState,
+                           episodeIdx,
+                           *grid->c_grid(),
+                           props->permeability());
 
-    if (output) {
-        report.reportParam(outStream);
-        warnIfUnusedParams(param);
+        if (episodeIdx == 0) {
+            // @@@ HACK: we should really make a new well state and
+            // properly transfer old well state to it every epoch,
+            // since number of wells may change etc.
+            well_state.init(wells.c_wells(), state);
+        }
+
+        simtimer.init(timeMap,
+                      /*beginReportStepIdx=*/episodeIdx,
+                      /*endReportStepIdx=*/episodeIdx + 1);
+
+        if (episodeIdx == 0)
+            outputWriter.writeInit(simtimer, state, well_state.basicWellState());
+
+        // Create and run simulator.
+        SimulatorFullyImplicitBlackoil simulator(param,
+                                                 *grid->c_grid(),
+                                                 *new_props,
+                                                 rock_comp->isActive() ? rock_comp.get() : 0,
+                                                 wells,
+                                                 linsolver,
+                                                 grav,
+                                                 outputWriter);
+        SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
+
+        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+        fullReport += episodeReport;
+
+        if (output) {
+            episodeReport.reportParam(outStream);
+        }
     }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
-    report.report(std::cout);
+    fullReport.report(std::cout);
 
     if (output) {
         std::string filename = output_dir + "/walltime.param";
         std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
-        report.reportParam(tot_os);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
     }
 #else
     std::cout << "\n\n================    Starting main simulation loop     ===============\n"
@@ -248,12 +268,13 @@ try
 
     SimulatorReport rep;
     // With a deck, we may have more epochs etc.
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     int step = 0;
     SimulatorTimer simtimer;
     // Use timer for last epoch to obtain total time.
     deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
     simtimer.init(*deck);
+
     const double total_time = simtimer.totalTime();
     for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
         // Set epoch index.
@@ -285,6 +306,9 @@ try
             well_state.init(wells.c_wells(), state);
         }
 
+        if (epoch == 0)
+            outputWriter.writeInit(simtimer, state, well_state.basicWellState());
+
         // Create and run simulator.
         SimulatorFullyImplicitBlackoil simulator(param,
                                                  *grid->c_grid(),
@@ -294,6 +318,8 @@ try
                                                  linsolver,
                                                  grav,
                                                  outputWriter);
+        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+
         if (epoch == 0) {
             warnIfUnusedParams(param);
         }
diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp
index 0c1e78580..11558b396 100644
--- a/examples/test_implicit_ad.cpp
+++ b/examples/test_implicit_ad.cpp
@@ -23,6 +23,7 @@
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAd.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -37,7 +38,6 @@
 #include <opm/core/utility/Units.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/initState.hpp>
 
 #include <boost/shared_ptr.hpp>
@@ -110,7 +110,7 @@ try
     initStateBasic(*g, props0, param, 0.0, state);
     initBlackoilSurfvol(*g, props0, state);
 
-    Opm::WellState well_state;
+    Opm::WellStateFullyImplicitBlackoil well_state;
     well_state.init(wells.get(), state);
 
     solver.step(1.0, state, well_state);
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp
index bd2982238..4052b2417 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp
@@ -24,12 +24,12 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/WellDensitySegmented.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/well_controls.h>
 
@@ -223,7 +223,7 @@ namespace {
     FullyImplicitBlackoilSolver::
     step(const double   dt,
          BlackoilState& x ,
-         WellState&     xw)
+         WellStateFullyImplicitBlackoil& xw)
     {
         const V pvdt = geo_.poreVolume() / dt;
 
@@ -332,7 +332,7 @@ namespace {
 
     FullyImplicitBlackoilSolver::SolutionState
     FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                               const WellStateFullyImplicitBlackoil& xw)
     {
         const int nc = grid_.number_of_cells;
         const int np = x.numPhases();
@@ -429,7 +429,7 @@ namespace {
 
     FullyImplicitBlackoilSolver::SolutionState
     FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                               const WellStateFullyImplicitBlackoil& xw)
     {
         const int nc = grid_.number_of_cells;
         const int np = x.numPhases();
@@ -615,7 +615,7 @@ namespace {
 
 
     void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state,
-                                                                     const WellState& xw)
+                                                                     const WellStateFullyImplicitBlackoil& xw)
     {
         // 1. Compute properties required by computeConnectionPressureDelta().
         //    Note that some of the complexity of this part is due to the function
@@ -687,7 +687,7 @@ namespace {
     FullyImplicitBlackoilSolver::
     assemble(const V&             pvdt,
              const BlackoilState& x   ,
-             const WellState&     xw  )
+             const WellStateFullyImplicitBlackoil& xw  )
     {
         // Create the primary variables.
         const SolutionState state = variableState(x, xw);
@@ -1148,7 +1148,7 @@ namespace {
 
     void FullyImplicitBlackoilSolver::updateState(const V& dx,
                                                   BlackoilState& state,
-                                                  WellState& well_state)
+                                                  WellStateFullyImplicitBlackoil& well_state)
     {
         const int np = fluid_.numPhases();
         const int nc = grid_.number_of_cells;
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
index 671169bb3..c9c8ac487 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
@@ -33,7 +33,7 @@ namespace Opm {
     class RockCompressibility;
     class LinearSolverInterface;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
 
 
     /// A fully implicit solver for the black-oil problem.
@@ -76,7 +76,7 @@ namespace Opm {
         void
         step(const double   dt    ,
              BlackoilState& state ,
-             WellState&     wstate);
+             WellStateFullyImplicitBlackoil&     wstate);
 
     private:
         // Types and enums
@@ -149,18 +149,18 @@ namespace Opm {
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         SolutionState
         variableState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         void
         computeAccum(const SolutionState& state,
                      const int            aix  );
 
         void computeWellConnectionPressures(const SolutionState& state,
-                                            const WellState& xw);
+                                            const WellStateFullyImplicitBlackoil& xw);
 
         void
         addOldWellEq(const SolutionState& state);
@@ -173,14 +173,14 @@ namespace Opm {
 
         void
         assemble(const V&             dtpv,
-                 const BlackoilState& x   ,
-                 const WellState&     xw  );
+                 const BlackoilState& x,
+                 const WellStateFullyImplicitBlackoil& xw);
 
         V solveJacobianSystem() const;
 
         void updateState(const V& dx,
                          BlackoilState& state,
-                         WellState& well_state);
+                         WellStateFullyImplicitBlackoil& well_state);
 
         std::vector<ADB>
         computePressures(const SolutionState& state) const;
diff --git a/opm/autodiff/SimulatorCompressibleAd.cpp b/opm/autodiff/SimulatorCompressibleAd.cpp
index 42df2fc66..4b954e891 100644
--- a/opm/autodiff/SimulatorCompressibleAd.cpp
+++ b/opm/autodiff/SimulatorCompressibleAd.cpp
@@ -496,13 +496,13 @@ namespace Opm
                       << std::endl;
             std::cout.precision(8);
 
-            watercut.push(timer.currentTime() + timer.currentStepLength(),
+            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
                           produced[0]/(produced[0] + produced[1]),
                           tot_produced[0]/tot_porevol_init);
             if (wells_) {
                 wellreport.push(props_, *wells_,
                                 state.pressure(), state.surfacevol(), state.saturation(),
-                                timer.currentTime() + timer.currentStepLength(),
+                                timer.simulationTimeElapsed() + timer.currentStepLength(),
                                 well_state.bhp(), well_state.perfRates());
             }
             sreport.total_time =  step_timer.secsSinceStart();
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
index 00246f90f..aae39d9ca 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
@@ -29,6 +29,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -48,7 +49,6 @@
 
 #include <opm/core/grid/ColumnExtract.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
 #include <boost/filesystem.hpp>
@@ -77,7 +77,7 @@ namespace Opm
 
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         // Data.
@@ -127,7 +127,7 @@ namespace Opm
 
     SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer,
                                                         BlackoilState& state,
-                                                        WellState& well_state)
+                                                        WellStateFullyImplicitBlackoil& well_state)
     {
         return pimpl_->run(timer, state, well_state);
     }
@@ -198,7 +198,7 @@ namespace Opm
             std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
         }
     }
-    static void outputWellStateMatlab(const Opm::WellState& well_state,
+    static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
                                   const int step,
                                   const std::string& output_dir)
     {
@@ -312,11 +312,8 @@ namespace Opm
 
     SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
                                                               BlackoilState& state,
-                                                              WellState& well_state)
+                                                              WellStateFullyImplicitBlackoil& well_state)
     {
-        eclipseWriter_.writeInit(timer, state, well_state);
-        eclipseWriter_.writeTimeStep(timer, state, well_state);
-
         // Initialisation.
         std::vector<double> porevol;
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -384,13 +381,9 @@ namespace Opm
 
                 // Optionally, check if well controls are satisfied.
                 if (check_well_controls_) {
-                    Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                      well_state.perfRates(),
-                                                      fractional_flows,
-                                                      well_resflows_phase);
                     std::cout << "Checking well conditions." << std::endl;
                     // For testing we set surface := reservoir
-                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
+                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_state.wellRates(), well_state.wellRates());
                     ++well_control_iteration;
                     if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
                         OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
@@ -424,13 +417,13 @@ namespace Opm
                 tstep_os.close();
             }
 
-            // write an output file for later inspection
-            if (output_) {
-                eclipseWriter_.writeTimeStep(timer, state, well_state);
-            }
-
             // advance to next timestep before reporting at this location
             ++timer;
+
+            // write an output file for later inspection
+            if (output_) {
+                eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
+            }
         }
 
         total_timer.stop();
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
index 2bb62de41..26b91b4e3 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
@@ -37,7 +37,7 @@ namespace Opm
     class LinearSolverInterface;
     class SimulatorTimer;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
     struct SimulatorReport;
 
     /// Class collecting all necessary components for a two-phase simulation.
@@ -84,7 +84,7 @@ namespace Opm
         /// \return                    simulation report, with timing data
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         class Impl;
diff --git a/opm/autodiff/SimulatorIncompTwophaseAd.cpp b/opm/autodiff/SimulatorIncompTwophaseAd.cpp
index 9e3f73af5..19cbe93e1 100644
--- a/opm/autodiff/SimulatorIncompTwophaseAd.cpp
+++ b/opm/autodiff/SimulatorIncompTwophaseAd.cpp
@@ -579,12 +579,12 @@ namespace Opm
                     dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
                         .solveGravity(&initial_porevol[0], stepsize, state);
                 }
-                watercut.push(timer.currentTime() + timer.currentStepLength(),
+                watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
                               produced[0]/(produced[0] + produced[1]),
                               tot_produced[0]/tot_porevol_init);
                 if (wells_) {
                     wellreport.push(props_, *wells_, state.saturation(),
-                                    timer.currentTime() + timer.currentStepLength(),
+                                    timer.simulationTimeElapsed() + timer.currentStepLength(),
                                     well_state.bhp(), well_state.perfRates());
                 }
             }
diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp
index 5ba38e87b..669c260c8 100644
--- a/opm/autodiff/WellDensitySegmented.cpp
+++ b/opm/autodiff/WellDensitySegmented.cpp
@@ -19,22 +19,23 @@
 
 #include <opm/autodiff/WellDensitySegmented.hpp>
 #include <opm/core/wells.h>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <numeric>
 #include <cmath>
 
 
-std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
-                                                                              const WellState& wstate,
-                                                                              const PhaseUsage& phase_usage,
-                                                                              const std::vector<double>& b_perf,
-                                                                              const std::vector<double>& rsmax_perf,
-                                                                              const std::vector<double>& rvmax_perf,
-                                                                              const std::vector<double>& z_perf,
-                                                                              const std::vector<double>& surf_dens,
-                                                                              const double gravity)
+std::vector<double>
+Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
+                                                          const WellStateFullyImplicitBlackoil& wstate,
+                                                          const PhaseUsage& phase_usage,
+                                                          const std::vector<double>& b_perf,
+                                                          const std::vector<double>& rsmax_perf,
+                                                          const std::vector<double>& rvmax_perf,
+                                                          const std::vector<double>& z_perf,
+                                                          const std::vector<double>& surf_dens,
+                                                          const double gravity)
 {
     // Verify that we have consistent input.
     const int np = wells.number_of_phases;
@@ -46,7 +47,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
     if (surf_dens.size() != size_t(wells.number_of_phases)) {
         OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage.");
     }
-    if (nperf*np != int(wstate.perfRates().size())) {
+    if (nperf*np != int(wstate.perfPhaseRates().size())) {
         OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
     }
     if (nperf*np != int(b_perf.size())) {
@@ -90,7 +91,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
                     q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
                 }
                 // Subtract outflow through perforation.
-                q_out_perf[perf*np + phase] -= wstate.perfRates()[perf*np + phase];
+                q_out_perf[perf*np + phase] -= wstate.perfPhaseRates()[perf*np + phase];
             }
         }
     }
diff --git a/opm/autodiff/WellDensitySegmented.hpp b/opm/autodiff/WellDensitySegmented.hpp
index 420cfdea9..82a4494b8 100644
--- a/opm/autodiff/WellDensitySegmented.hpp
+++ b/opm/autodiff/WellDensitySegmented.hpp
@@ -27,7 +27,7 @@ struct Wells;
 namespace Opm
 {
 
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
     struct PhaseUsage;
 
 
@@ -51,7 +51,7 @@ namespace Opm
         /// \param[in] surf_dens    surface densities for active components, size P
         /// \param[in] gravity      gravity acceleration constant
         static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
-                                                                  const WellState& wstate,
+                                                                  const WellStateFullyImplicitBlackoil& wstate,
                                                                   const PhaseUsage& phase_usage,
                                                                   const std::vector<double>& b_perf,
                                                                   const std::vector<double>& rsmax_perf,
diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp
new file mode 100644
index 000000000..f2430b1b8
--- /dev/null
+++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp
@@ -0,0 +1,101 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+
+
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
+#include <opm/core/simulator/WellState.hpp>
+#include <vector>
+#include <cassert>
+
+namespace Opm
+{
+
+    /// The state of a set of wells, tailored for use by the fully
+    /// implicit blackoil simulator.
+    class WellStateFullyImplicitBlackoil
+    {
+    public:
+        /// Allocate and initialize if wells is non-null.  Also tries
+        /// to give useful initial values to the bhp(), wellRates()
+        /// and perfPhaseRates() fields, depending on controls
+        template <class State>
+        void init(const Wells* wells, const State& state)
+        {
+            if (wells == 0) {
+                return;
+            }
+
+            // We use the WellState::init() function to do bhp and well rates init.
+            // The alternative would be to copy that function wholesale.
+            basic_well_state_.init(wells, state);
+
+            // Initialize perfphaserates_, which must be done here.
+            const int nw = wells->number_of_wells;
+            const int np = wells->number_of_phases;
+            const int nperf = wells->well_connpos[nw];
+            perfphaserates_.resize(nperf * np, 0.0);
+            for (int w = 0; w < nw; ++w) {
+                assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
+                const WellControls* ctrl = wells->ctrls[w];
+                if (well_controls_well_is_shut(ctrl)) {
+                    // Shut well: perfphaserates_ are all zero.
+                } else {
+                    // Open well: Initialize perfphaserates_ to well
+                    // rates divided by the number of perforations.
+                    const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
+                    for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                        for (int p = 0; p < np; ++p) {
+                            perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
+                        }
+                    }
+                }
+            }
+        }
+
+        /// One bhp pressure per well.
+        std::vector<double>& bhp() { return basic_well_state_.bhp(); }
+        const std::vector<double>& bhp() const { return basic_well_state_.bhp(); }
+
+        /// One rate per well and phase.
+        std::vector<double>& wellRates() { return basic_well_state_.wellRates(); }
+        const std::vector<double>& wellRates() const { return basic_well_state_.wellRates(); }
+
+        /// One rate per phase and well connection.
+        std::vector<double>& perfPhaseRates() { return perfphaserates_; }
+        const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
+
+        /// For interfacing with functions that take a WellState.
+        const WellState& basicWellState() const
+        {
+            return basic_well_state_;
+        }
+
+    private:
+        WellState basic_well_state_;
+        std::vector<double> perfphaserates_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp
index a845c0fde..0ab4e38a5 100644
--- a/tests/test_welldensitysegmented.cpp
+++ b/tests/test_welldensitysegmented.cpp
@@ -28,7 +28,7 @@
 
 #include <opm/autodiff/WellDensitySegmented.hpp>
 #include <opm/core/wells.h>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/utility/Units.hpp>
 
@@ -65,8 +65,8 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
                                   1.0, 0.0, 0.0,
                                   1.0, 0.0, 0.0,
                                   1.0, 0.0, 0.0 };
-    WellState wellstate;
-    wellstate.perfRates() = rates;
+    WellStateFullyImplicitBlackoil wellstate;
+    wellstate.perfPhaseRates() = rates;
     PhaseUsage pu;
     pu.num_phases = 3;
     pu.phase_used[0] = true;