diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index f0b3cd2f7..d18971f0c 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -47,6 +47,9 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/VFPProdProperties.cpp
 	opm/autodiff/VFPInjProperties.cpp
 	opm/autodiff/WellMultiSegment.cpp
+        opm/autodiff/BlackoilSolventState.cpp
+        opm/polymer/PolymerState.cpp
+        opm/polymer/PolymerBlackoilState.cpp
 	opm/polymer/CompressibleTpfaPolymer.cpp
 	opm/polymer/IncompTpfaPolymer.cpp
 	opm/polymer/PolymerInflow.cpp
diff --git a/examples/sim_2p_incomp_ad.cpp b/examples/sim_2p_incomp_ad.cpp
index e4e203228..8e671a6ee 100644
--- a/examples/sim_2p_incomp_ad.cpp
+++ b/examples/sim_2p_incomp_ad.cpp
@@ -105,8 +105,9 @@ try
     std::unique_ptr<GridManager> grid;
     std::unique_ptr<IncompPropertiesInterface> props;
     std::unique_ptr<RockCompressibility> rock_comp;
+    std::unique_ptr<TwophaseState> state;
     EclipseStateConstPtr eclipseState;
-    TwophaseState state;
+
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
@@ -118,19 +119,25 @@ try
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // check_well_controls = param.getDefault("check_well_controls", false);
-        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(deck, eclipseState));
-        // Gravity.
-        gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
-        // Init state variables (saturation and pressure).
-        if (param.has("init_saturation")) {
-            initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init
+            props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+
+            state.reset( new TwophaseState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(deck, eclipseState));
+            // Gravity.
+            gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+            // Init state variables (saturation and pressure).
+            if (param.has("init_saturation")) {
+                initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
         }
     } else {
         // Grid init.
@@ -141,14 +148,20 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+
+            state.reset( new TwophaseState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+        }
     }
 
     // Warn if gravity but no density difference.
@@ -170,7 +183,7 @@ try
         // terms of total pore volume.
         std::vector<double> porevol;
         if (rock_comp->isActive()) {
-            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
         } else {
             computePorevolume(*grid->c_grid(), props->porosity(), porevol);
         }
@@ -236,8 +249,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more report steps etc.
         WellState well_state;
@@ -255,7 +268,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             simtimer.setCurrentStepNum(reportStepIdx);
@@ -273,7 +286,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
             if (output) {
                 epoch_rep.reportParam(epoch_os);
             }
diff --git a/examples/sim_poly2p_comp_reorder.cpp b/examples/sim_poly2p_comp_reorder.cpp
index 175c7019b..b4fd352ce 100644
--- a/examples/sim_poly2p_comp_reorder.cpp
+++ b/examples/sim_poly2p_comp_reorder.cpp
@@ -92,7 +92,7 @@ try
     boost::scoped_ptr<RockCompressibility> rock_comp;
     Opm::DeckConstPtr deck;
     EclipseStateConstPtr eclipseState;
-    PolymerBlackoilState state;
+    std::unique_ptr<PolymerBlackoilState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -106,23 +106,29 @@ try
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // check_well_controls = param.getDefault("check_well_controls", false);
-        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(deck, eclipseState));
-        // Gravity.
-        gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
-        // Init state variables (saturation and pressure).
-        if (param.has("init_saturation")) {
-            initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init
+            props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, ug_grid));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+
+            state.reset( new PolymerBlackoilState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(deck, eclipseState));
+            // Gravity.
+            gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+            // Init state variables (saturation and pressure).
+            if (param.has("init_saturation")) {
+                initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
+            initBlackoilSurfvol(ug_grid, *props, *state);
+            // Init polymer properties.
+            poly_props.readFromDeck(deck, eclipseState);
         }
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
-        // Init polymer properties.
-        poly_props.readFromDeck(deck, eclipseState);
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
@@ -132,31 +138,43 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
-        // Init Polymer state
-        if (param.has("poly_init")) {
-            double poly_init = param.getDefault("poly_init", 0.0);
-            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-                double smin[2], smax[2];
-                props->satRange(1, &cell, smin, smax);
-                if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                    state.concentration()[cell] = poly_init;
-                    state.maxconcentration()[cell] = poly_init;
-                } else {
-                    state.saturation()[2*cell + 0] = 0.;
-                    state.saturation()[2*cell + 1] = 1.;
-                    state.concentration()[cell] = 0.;
-                    state.maxconcentration()[cell] = 0.;
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new BlackoilPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+            state.reset( new PolymerBlackoilState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            initBlackoilSurfvol(ug_grid, *props, *state);
+            // Init Polymer state
+
+            if (param.has("poly_init")) {
+                double poly_init = param.getDefault("poly_init", 0.0);
+                for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
+                    double smin[2], smax[2];
+
+                    auto& saturation = state->saturation();
+                    auto& concentration = state->getCellData( state->CONCENTRATION );
+                    auto& max_concentration = state->getCellData( state->CMAX );
+
+                    props->satRange(1, &cell, smin, smax);
+                    if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                        concentration[cell] = poly_init;
+                        max_concentration[cell] = poly_init;
+                    } else {
+                        saturation[2*cell + 0] = 0.;
+                        saturation[2*cell + 1] = 1.;
+                        concentration[cell] = 0.;
+                        max_concentration[cell] = 0.;
+                    }
                 }
             }
+
         }
         // Init polymer properties.
         // Setting defaults to provide a simple example case.
@@ -240,8 +258,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more epochs etc.
         WellState well_state;
@@ -284,7 +302,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             // Create and run simulator.
@@ -300,7 +318,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
 
             // Update total timing report and remember step number.
             rep += epoch_rep;
diff --git a/examples/sim_poly2p_incomp_reorder.cpp b/examples/sim_poly2p_incomp_reorder.cpp
index cbb234eea..3826a06d8 100644
--- a/examples/sim_poly2p_incomp_reorder.cpp
+++ b/examples/sim_poly2p_incomp_reorder.cpp
@@ -92,7 +92,7 @@ try
     boost::scoped_ptr<IncompPropertiesInterface> props;
     boost::scoped_ptr<RockCompressibility> rock_comp;
     EclipseStateConstPtr eclipseState;
-    PolymerState state;
+    std::unique_ptr<PolymerState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -106,22 +106,27 @@ try
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // check_well_controls = param.getDefault("check_well_controls", false);
-        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(deck, eclipseState));
-        // Gravity.
-        gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
-        // Init state variables (saturation and pressure).
-        if (param.has("init_saturation")) {
-            initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+            // Rock and fluid init
+            props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid ));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+            state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(deck, eclipseState));
+            // Gravity.
+            gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+            // Init state variables (saturation and pressure).
+            if (param.has("init_saturation")) {
+                initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
+            // Init polymer properties.
+            poly_props.readFromDeck(deck, eclipseState);
         }
-        // Init polymer properties.
-        poly_props.readFromDeck(deck, eclipseState);
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
@@ -131,28 +136,38 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        // Init Polymer state
-        if (param.has("poly_init")) {
-            double poly_init = param.getDefault("poly_init", 0.0);
-            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-                double smin[2], smax[2];
-                props->satRange(1, &cell, smin, smax);
-                if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                    state.concentration()[cell] = poly_init;
-                    state.maxconcentration()[cell] = poly_init;
-                } else {
-                    state.saturation()[2*cell + 0] = 0.;
-                    state.saturation()[2*cell + 1] = 1.;
-                    state.concentration()[cell] = 0.;
-                    state.maxconcentration()[cell] = 0.;
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));;
+            state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            // Init Polymer state
+            if (param.has("poly_init")) {
+                double poly_init = param.getDefault("poly_init", 0.0);
+                for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
+                    double smin[2], smax[2];
+
+                    auto& saturation = state->saturation();
+                    auto& concentration = state->getCellData( state->CONCENTRATION );
+                    auto& max_concentration = state->getCellData( state->CMAX );
+
+                    props->satRange(1, &cell, smin, smax);
+                    if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                        concentration[cell] = poly_init;
+                        max_concentration[cell] = poly_init;
+                    } else {
+                        saturation[2*cell + 0] = 0.;
+                        saturation[2*cell + 1] = 1.;
+                        concentration[cell] = 0.;
+                        max_concentration[cell] = 0.;
+                    }
                 }
             }
         }
@@ -212,7 +227,7 @@ try
         // terms of total pore volume.
         std::vector<double> porevol;
         if (rock_comp->isActive()) {
-            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
         } else {
             computePorevolume(*grid->c_grid(), props->porosity(), porevol);
         }
@@ -276,8 +291,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more epochs etc.
 
@@ -321,7 +336,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             // Create and run simulator.
@@ -339,7 +354,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
 
             // Update total timing report and remember step number.
             rep += epoch_rep;
diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp
index eb266ef6a..4b17c59c7 100644
--- a/examples/sim_poly_fi2p_comp_ad.cpp
+++ b/examples/sim_poly_fi2p_comp_ad.cpp
@@ -115,7 +115,7 @@ try
     std::shared_ptr<BlackoilPropertiesInterface> props;
     std::shared_ptr<BlackoilPropsAdInterface> new_props;
     std::shared_ptr<RockCompressibility> rock_comp;
-    PolymerBlackoilState state;
+    std::unique_ptr<PolymerBlackoilState> state;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
@@ -187,6 +187,8 @@ try
                                                 Opm::UgGridHelpers::globalCell(cGrid),
                                                 Opm::UgGridHelpers::cartDims(cGrid),
                                                 param));
+
+    state.reset( new PolymerBlackoilState( Opm::UgGridHelpers::numCells(cGrid), Opm::UgGridHelpers::numFaces(cGrid), 2));
     new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, materialLawManager, cGrid));
     PolymerProperties polymer_props(deck, eclipseState);
     PolymerPropsAd polymer_props_ad(polymer_props);
@@ -199,10 +201,10 @@ try
 
     // Init state variables (saturation and pressure).
     if (param.has("init_saturation")) {
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
+        initStateBasic(*grid->c_grid(), *props, param, gravity[2], *state);
+        initBlackoilSurfvol(*grid->c_grid(), *props, *state);
     } else {
-        initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], *state);
     }
 
     bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
@@ -254,7 +256,7 @@ try
                   deck,
                   *fis_solver,
                   grav);
-    fullReport= simulator.run(simtimer, state);
+    fullReport= simulator.run(simtimer, *state);
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     fullReport.report(std::cout);
diff --git a/opm/autodiff/BackupRestore.hpp b/opm/autodiff/BackupRestore.hpp
index 70fe4e9b6..0885257a3 100644
--- a/opm/autodiff/BackupRestore.hpp
+++ b/opm/autodiff/BackupRestore.hpp
@@ -21,8 +21,9 @@
 #define OPM_BACKUPRESTORE_HEADER_INCLUDED
 
 #include <iostream>
+#include <opm/common/data/SimulationDataContainer.hpp>
+
 
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
@@ -159,7 +160,7 @@ namespace Opm {
                BlackoilStateId = 2,
                WellStateFullyImplicitBackoilId = 3 };
 
-        inline int objectId( const SimulatorState& /* state */) {
+        inline int objectId( const SimulationDataContainer& /* state */) {
             return SimulatorStateId;
         }
         inline int objectId( const WellState& /* state */) {
@@ -184,9 +185,9 @@ namespace Opm {
         }
     }
 
-    // SimulatorState
+    // SimulationDataContainer
     inline
-    std::ostream& operator << (std::ostream& out, const SimulatorState& state )
+    std::ostream& operator << (std::ostream& out, const SimulationDataContainer& state )
     {
         // write id of object to stream
         writeValue( out, objectId( state ) );
@@ -205,7 +206,7 @@ namespace Opm {
     }
 
     inline
-    std::istream& operator >> (std::istream& in, SimulatorState& state )
+    std::istream& operator >> (std::istream& in, SimulationDataContainer& state )
     {
         // check id of stored object
         checkObjectId( in, state );
@@ -213,7 +214,7 @@ namespace Opm {
         int numPhases = 0;
         readValue( in, numPhases );
 
-        if( numPhases != state.numPhases() )
+        if( numPhases != (int) state.numPhases() )
             OPM_THROW(std::logic_error,"num phases wrong");
 
         // read variables
@@ -234,7 +235,7 @@ namespace Opm {
         writeValue( out, objectId( state ) );
 
         // backup simulator state
-        const SimulatorState& simstate = static_cast< const SimulatorState& > (state);
+        const SimulationDataContainer& simstate = static_cast< const SimulationDataContainer& > (state);
         out << simstate;
 
         // backup additional blackoil state variables
@@ -252,7 +253,7 @@ namespace Opm {
         checkObjectId( in, state );
 
         // restore simulator state
-        SimulatorState& simstate = static_cast< SimulatorState& > (state);
+        SimulationDataContainer& simstate = static_cast< SimulationDataContainer& > (state);
         in >> simstate;
 
         // restore additional blackoil state variables
diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp
index 9eceb3565..2551bd7fb 100644
--- a/opm/autodiff/BlackoilModelBase.hpp
+++ b/opm/autodiff/BlackoilModelBase.hpp
@@ -45,7 +45,7 @@ namespace Opm {
     class RockCompressibility;
     class NewtonIterationBlackoilInterface;
     class VFPProperties;
-
+    class SimulationDataContainer;
 
     /// Struct for containing iteration variables.
     struct DefaultBlackoilSolutionState
@@ -207,7 +207,7 @@ namespace Opm {
 
         /// \brief compute the relative change between to simulation states
         //  \return || u^n+1 - u^n || / || u^n+1 ||
-        double relativeChange( const SimulatorState& previous, const SimulatorState& current ) const;
+        double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const;
 
         /// The size (number of unknowns) of the nonlinear system of equations.
         int sizeNonLinear() const;
diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index c0b3fefad..ad8a32562 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -48,6 +48,7 @@
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
 
+#include <opm/common/data/SimulationDataContainer.hpp>
 #include <cassert>
 #include <cmath>
 #include <iostream>
@@ -2520,8 +2521,8 @@ namespace detail {
     template <class Grid, class Implementation>
     double
     BlackoilModelBase<Grid, Implementation>::
-    relativeChange(const SimulatorState& previous,
-                   const SimulatorState& current ) const
+    relativeChange(const SimulationDataContainer& previous,
+                   const SimulationDataContainer& current ) const
     {
         std::vector< double > p0  ( previous.pressure() );
         std::vector< double > sat0( previous.saturation() );
diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp
index 3f703a045..fb4688caf 100644
--- a/opm/autodiff/BlackoilSolventModel_impl.hpp
+++ b/opm/autodiff/BlackoilSolventModel_impl.hpp
@@ -137,11 +137,14 @@ namespace Opm {
 
         // Initial polymer concentration.
         if (has_solvent_) {
-            assert (not x.solvent_saturation().empty());
-            const int nc = x.solvent_saturation().size();
-            const V ss = Eigen::Map<const V>(&x.solvent_saturation()[0], nc);
+            const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL );
+            const int nc = solvent_saturation.size();
+            const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
+
+            // This is some insanely detailed flickety flackety code;
             // Solvent belongs after other reservoir vars but before well vars.
             auto solvent_pos = vars0.begin() + fluid_.numPhases();
+            assert (not solvent_saturation.empty());
             assert(solvent_pos == vars0.end() - 2);
             vars0.insert(solvent_pos, ss);
         }
@@ -548,9 +551,10 @@ namespace Opm {
             Base::updateState(modified_dx, reservoir_state, well_state);
 
             // Update solvent.
-            const V ss_old = Eigen::Map<const V>(&reservoir_state.solvent_saturation()[0], nc, 1);
+            auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
+            const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
             const V ss = (ss_old - dss).max(zero);
-            std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin());
+            std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
 
             // adjust oil saturation
             const Opm::PhaseUsage& pu = fluid_.phaseUsage();
diff --git a/opm/autodiff/BlackoilSolventState.cpp b/opm/autodiff/BlackoilSolventState.cpp
new file mode 100644
index 000000000..681bd1ca2
--- /dev/null
+++ b/opm/autodiff/BlackoilSolventState.cpp
@@ -0,0 +1,34 @@
+/*
+  Copyright 2015 IRIS AS, 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/>.
+*/
+
+
+#include <opm/autodiff/BlackoilSolventState.hpp>
+
+namespace Opm
+{
+    const std::string BlackoilSolventState::SSOL = "SSOL";
+
+    BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases)
+        : BlackoilState( number_of_cells , number_of_faces , number_of_phases)
+    {
+        registerCellData( SSOL , 1 );
+    };
+
+} // namespace Opm
+
diff --git a/opm/autodiff/BlackoilSolventState.hpp b/opm/autodiff/BlackoilSolventState.hpp
index cf5246443..0103e0904 100644
--- a/opm/autodiff/BlackoilSolventState.hpp
+++ b/opm/autodiff/BlackoilSolventState.hpp
@@ -20,9 +20,9 @@
 #ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
 #define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
 
-#include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/grid.h>
-#include <vector>
+#include <string>
+
+ #include <opm/core/simulator/BlackoilState.hpp>
 
 
 namespace Opm
@@ -33,22 +33,9 @@ namespace Opm
     class BlackoilSolventState : public BlackoilState
     {
     public:
-        void init(const UnstructuredGrid& g, int num_phases)
-        {
-            this->init(g.number_of_cells, g.number_of_faces, num_phases);
-        }
+        static const std::string SSOL;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            BlackoilState::init(number_of_cells, number_of_faces, num_phases);
-            solventId_  = SimulatorState::registerCellData( "SSOL", 1 );
-        }
-
-        std::vector<double>& solvent_saturation()    { return cellData()[ solventId_ ]; }
-        const std::vector<double>& solvent_saturation() const    { return cellData()[ solventId_ ]; }
-
-    private:
-        int solventId_;
+        BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases);
     };
 
 } // namespace Opm
diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp
index b28703ffc..b8af7d6a4 100644
--- a/opm/autodiff/FlowMain.hpp
+++ b/opm/autodiff/FlowMain.hpp
@@ -183,7 +183,8 @@ namespace Opm
         bool use_local_perm_ = true;
         std::unique_ptr<DerivedGeology> geoprops_;
         // setupState()
-        ReservoirState state_;
+        std::unique_ptr<ReservoirState> state_;
+
         std::vector<double> threshold_pressures_;
         // distributeData()
         boost::any parallel_information_;
@@ -441,8 +442,13 @@ namespace Opm
                                               Opm::UgGridHelpers::cartDims(grid),
                                               param_);
 
+
             // Init state variables (saturation and pressure).
             if (param_.has("init_saturation")) {
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases() ));
+
                 initStateBasic(Opm::UgGridHelpers::numCells(grid),
                                Opm::UgGridHelpers::globalCell(grid),
                                Opm::UgGridHelpers::cartDims(grid),
@@ -451,26 +457,35 @@ namespace Opm
                                Opm::UgGridHelpers::beginFaceCentroids(grid),
                                Opm::UgGridHelpers::beginCellCentroids(grid),
                                Opm::UgGridHelpers::dimensions(grid),
-                               props, param_, gravity_[2], state_);
+                               props, param_, gravity_[2], *state_);
 
-                initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state_);
+                initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_);
 
                 enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
                 if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
                     const int numPhases = props.numPhases();
                     const int numCells  = Opm::UgGridHelpers::numCells(grid);
+
+                    // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState.
+                    auto& gor = state_->getCellData( BlackoilState::GASOILRATIO );
+                    const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL );
                     for (int c = 0; c < numCells; ++c) {
-                        state_.gasoilratio()[c] = state_.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
-                            / state_.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
+                        // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field.
+                        gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]];
                     }
                 }
             } else if (deck_->hasKeyword("EQUIL") && props.numPhases() == 3) {
-                state_.init(Opm::UgGridHelpers::numCells(grid),
-                           Opm::UgGridHelpers::numFaces(grid),
-                           props.numPhases());
-                initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], state_);
-                state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
+                // Which state class are we really using - what a f... mess?
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases()));
+
+                initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], *state_);
+                //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
             } else {
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases()));
                 initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
                                           Opm::UgGridHelpers::globalCell(grid),
                                           Opm::UgGridHelpers::numFaces(grid),
@@ -478,12 +493,12 @@ namespace Opm
                                           Opm::UgGridHelpers::beginFaceCentroids(grid),
                                           Opm::UgGridHelpers::beginCellCentroids(grid),
                                           Opm::UgGridHelpers::dimensions(grid),
-                                          props, deck_, gravity_[2], state_);
+                                          props, deck_, gravity_[2], *state_);
             }
 
             // Threshold pressures.
             std::map<std::pair<int, int>, double> maxDp;
-            computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]);
+            computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]);
             threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp);
             std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipse_state_, geoprops_->nnc(), maxDp);
             threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
@@ -493,9 +508,9 @@ namespace Opm
                 const int numCells = Opm::UgGridHelpers::numCells(grid);
                 std::vector<int> cells(numCells);
                 for (int c = 0; c < numCells; ++c) { cells[c] = c; }
-                std::vector<double> pc = state_.saturation();
-                props.capPress(numCells, state_.saturation().data(), cells.data(), pc.data(), nullptr);
-                fluidprops_->setSwatInitScaling(state_.saturation(), pc);
+                std::vector<double> pc = state_->saturation();
+                props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr);
+                fluidprops_->setSwatInitScaling(state_->saturation(), pc);
             }
         }
 
@@ -518,7 +533,7 @@ namespace Opm
             // and initilialize new properties and states for it.
             if (must_distribute_) {
                 distributeGridAndData(grid_init_->grid(), deck_, eclipse_state_,
-                                      state_, *fluidprops_, *geoprops_,
+                                      *state_, *fluidprops_, *geoprops_,
                                       material_law_manager_, threshold_pressures_,
                                       parallel_information_, use_local_perm_);
             }
@@ -617,7 +632,7 @@ namespace Opm
                               << std::flush;
                 }
 
-                SimulatorReport fullReport = simulator_->run(simtimer, state_);
+                SimulatorReport fullReport = simulator_->run(simtimer, *state_);
 
                 if (output_cout_) {
                     std::cout << "\n\n================    End of simulation     ===============\n\n";
diff --git a/opm/autodiff/ParallelDebugOutput.hpp b/opm/autodiff/ParallelDebugOutput.hpp
index 2398a6f28..dabe9b4be 100644
--- a/opm/autodiff/ParallelDebugOutput.hpp
+++ b/opm/autodiff/ParallelDebugOutput.hpp
@@ -19,8 +19,10 @@
 #ifndef OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
 #define OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+
 #include <opm/core/grid.h>
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/wells/WellsManager.hpp>
 
@@ -35,6 +37,7 @@ namespace Opm
 {
 
     class ParallelDebugOutputInterface
+
     {
     protected:
         ParallelDebugOutputInterface () {}
@@ -42,11 +45,11 @@ namespace Opm
         virtual ~ParallelDebugOutputInterface() {}
 
         // gather solution to rank 0 for EclipseWriter
-        virtual bool collectToIORank( const SimulatorState& localReservoirState,
+        virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
                                       const int reportStep ) = 0;
 
-        virtual const SimulatorState& globalReservoirState() const = 0 ;
+        virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
         virtual const WellState& globalWellState() const = 0 ;
         virtual bool isIORank() const = 0;
         virtual bool isParallel() const = 0;
@@ -60,7 +63,7 @@ namespace Opm
     protected:
         const GridImpl& grid_;
 
-        const SimulatorState* globalState_;
+        const SimulationDataContainer* globalState_;
         const WellState*      wellState_;
 
     public:
@@ -71,7 +74,7 @@ namespace Opm
             : grid_( grid ) {}
 
         // gather solution to rank 0 for EclipseWriter
-        virtual bool collectToIORank( const SimulatorState& localReservoirState,
+        virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
                                       const int /* reportStep */)
         {
@@ -80,7 +83,7 @@ namespace Opm
             return true ;
         }
 
-        virtual const SimulatorState& globalReservoirState() const { return *globalState_; }
+        virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; }
         virtual const WellState& globalWellState() const { return *wellState_; }
         virtual bool isIORank () const { return true; }
         virtual bool isParallel () const { return false; }
@@ -243,7 +246,7 @@ namespace Opm
                     Dune::CpGrid& globalGrid = *grid_;
 
                     // initialize global state with correct sizes
-                    globalReservoirState_.init( globalGrid.numCells(), globalGrid.numFaces(), numPhases );
+                    globalReservoirState_.reset( new SimulationDataContainer( globalGrid.numCells(), globalGrid.numFaces(), numPhases ));
 
                     // copy global cartesian index
                     globalIndex_ = globalGrid.globalCell();
@@ -306,18 +309,18 @@ namespace Opm
             }
         }
 
-        class PackUnPackSimulatorState : public P2PCommunicatorType::DataHandleInterface
+        class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface
         {
-            const SimulatorState& localState_;
-            SimulatorState& globalState_;
+            const SimulationDataContainer& localState_;
+            SimulationDataContainer& globalState_;
             const WellState& localWellState_;
             WellState& globalWellState_;
             const IndexMapType& localIndexMap_;
             const IndexMapStorageType& indexMaps_;
 
         public:
-            PackUnPackSimulatorState( const SimulatorState& localState,
-                                      SimulatorState& globalState,
+            PackUnPackSimulationDataContainer( const SimulationDataContainer& localState,
+                                      SimulationDataContainer& globalState,
                                       const WellState& localWellState,
                                       WellState& globalWellState,
                                       const IndexMapType& localIndexMap,
@@ -333,12 +336,11 @@ namespace Opm
                 if( isIORank )
                 {
                     // add missing data to global state
-                    for( size_t i=globalState_.cellData().size();
-                         i<localState.cellData().size(); ++i )
-                    {
-                        const size_t components = localState.cellData()[ i ].size() / localState.numCells();
-                        assert( components * localState.numCells() == localState.cellData()[ i ].size() );
-                        globalState_.registerCellData( localState.cellDataNames()[ i ], components );
+   		    for (const auto& pair : localState.cellData()) {
+		        const std::string& key = pair.first;
+                        if (!globalState_.hasCellData( key )) {
+			    globalState_.registerCellData( key , localState.numCellDataComponents( key ));
+                        }
                     }
 
                     MessageBufferType buffer;
@@ -357,13 +359,10 @@ namespace Opm
                 }
 
                 // write all cell data registered in local state
-                const size_t numCells = localState_.numCells();
-                const size_t numCellData = localState_.cellData().size();
-                for( size_t d=0; d<numCellData; ++d )
-                {
-                    const std::vector< double >& data = localState_.cellData()[ d ];
-                    const size_t stride = data.size() / numCells ;
-                    assert( numCells * stride == data.size() );
+		for (const auto& pair : localState_.cellData()) {
+		    const std::string& key = pair.first;
+                    const auto& data = pair.second;
+                    const size_t stride = localState_.numCellDataComponents( key );
 
                     for( size_t i=0; i<stride; ++i )
                     {
@@ -378,22 +377,20 @@ namespace Opm
 
             void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
             {
-                // read all cell data registered in local state
-                const size_t numCells = globalState_.numCells();
-                const size_t numCellData = globalState_.cellData().size();
-                for( size_t d=0; d<numCellData; ++d )
-                {
-                    std::vector< double >& data = globalState_.cellData()[ d ];
-                    const size_t stride = data.size() / numCells ;
-                    assert( numCells * stride == data.size() );
+                // write all cell data registered in local state
+	        for (auto& pair : globalState_.cellData()) {
+                    const std::string& key = pair.first;
+		    auto& data = pair.second;
+                    const size_t stride = globalState_.numCellDataComponents( key );
 
                     for( size_t i=0; i<stride; ++i )
                     {
-                        // write all data from local state to buffer
-                        read( buffer, indexMap, data, i, stride );
+                        //write all data from local state to buffer
+		        read( buffer, indexMap, data, i, stride );
                     }
                 }
 
+
                 // read well data from buffer
                 readWells( buffer );
             }
@@ -516,7 +513,7 @@ namespace Opm
         };
 
         // gather solution to rank 0 for EclipseWriter
-        bool collectToIORank( const SimulatorState& localReservoirState,
+        bool collectToIORank( const SimulationDataContainer& localReservoirState,
                               const WellState& localWellState,
                               const int reportStep )
         {
@@ -536,10 +533,10 @@ namespace Opm
                                            false);
 
                 const Wells* wells = wells_manager.c_wells();
-                globalWellState_.init(wells, globalReservoirState_, globalWellState_ );
+                globalWellState_.init(wells, *globalReservoirState_, globalWellState_ );
             }
 
-            PackUnPackSimulatorState packUnpack( localReservoirState, globalReservoirState_,
+            PackUnPackSimulationDataContainer packUnpack( localReservoirState, *globalReservoirState_,
                                                  localWellState, globalWellState_,
                                                  localIndexMap_, indexMaps_, isIORank() );
 
@@ -552,7 +549,7 @@ namespace Opm
             return isIORank();
         }
 
-        const SimulatorState& globalReservoirState() const { return globalReservoirState_; }
+        const SimulationDataContainer& globalReservoirState() const { return *globalReservoirState_; }
         const WellState& globalWellState() const { return globalWellState_; }
 
         bool isIORank() const
@@ -573,18 +570,18 @@ namespace Opm
         }
 
     protected:
-        std::unique_ptr< Dune::CpGrid > grid_;
-        Opm::EclipseStateConstPtr       eclipseState_;
-        const double*                   permeability_;
-        P2PCommunicatorType             toIORankComm_;
-        IndexMapType                    globalIndex_;
-        IndexMapType                    localIndexMap_;
-        IndexMapStorageType             indexMaps_;
-        SimulatorState                  globalReservoirState_;
+        std::unique_ptr< Dune::CpGrid >           grid_;
+        Opm::EclipseStateConstPtr                 eclipseState_;
+        const double*                             permeability_;
+        P2PCommunicatorType                       toIORankComm_;
+        IndexMapType                              globalIndex_;
+        IndexMapType                              localIndexMap_;
+        IndexMapStorageType                       indexMaps_;
+        std::unique_ptr<SimulationDataContainer>  globalReservoirState_;
         // this needs to be revised
-        WellStateFullyImplicitBlackoil  globalWellState_;
+        WellStateFullyImplicitBlackoil            globalWellState_;
         // true if we are on I/O rank
-        const bool                      isIORank_;
+        const bool                                isIORank_;
     };
 #endif // #if HAVE_DUNE_CORNERPOINT
 
diff --git a/opm/autodiff/RedistributeDataHandles.hpp b/opm/autodiff/RedistributeDataHandles.hpp
index e401f00d0..e2c5724c9 100644
--- a/opm/autodiff/RedistributeDataHandles.hpp
+++ b/opm/autodiff/RedistributeDataHandles.hpp
@@ -232,7 +232,7 @@ public:
     {
         assert( T::codimension == 0);
 
-        for ( int i=0; i<sendState_.numPhases(); ++i )
+        for ( size_t i=0; i<sendState_.numPhases(); ++i )
         {
             buffer.write(sendState_.surfacevol()[e.index()*sendState_.numPhases()+i]);
         }
@@ -240,7 +240,7 @@ public:
         buffer.write(sendState_.rv()[e.index()]);
         buffer.write(sendState_.pressure()[e.index()]);
         buffer.write(sendState_.temperature()[e.index()]);
-        for ( int i=0; i<sendState_.numPhases(); ++i )
+        for ( size_t i=0; i<sendState_.numPhases(); ++i )
         {
             buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
         }
@@ -257,11 +257,11 @@ public:
     void scatter(B& buffer, const T& e, std::size_t size_arg)
     {
         assert( T::codimension == 0);
-        assert( int(size_arg) == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
+        assert( size_arg == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
         static_cast<void>(size_arg);
 
         double val;
-        for ( int i=0; i<recvState_.numPhases(); ++i )
+        for ( size_t i=0; i<recvState_.numPhases(); ++i )
         {
             buffer.read(val);
             recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
@@ -274,7 +274,7 @@ public:
         recvState_.pressure()[e.index()]=val;
         buffer.read(val);
         recvState_.temperature()[e.index()]=val;
-        for ( int i=0; i<recvState_.numPhases(); ++i )
+        for ( size_t i=0; i<recvState_.numPhases(); ++i )
         {
             buffer.read(val);
             recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
@@ -411,8 +411,6 @@ void distributeGridAndData( Dune::CpGrid& grid,
                             boost::any& parallelInformation,
                             const bool useLocalPerm)
 {
-    BlackoilState distributed_state;
-
     Dune::CpGrid global_grid ( grid );
     global_grid.switchToGlobalView();
 
@@ -441,8 +439,8 @@ void distributeGridAndData( Dune::CpGrid& grid,
     BlackoilPropsAdFromDeck distributed_props(properties,
                                               distributed_material_law_manager,
                                               grid.numCells());
-    distributed_state.init(grid.numCells(), grid.numFaces(), state.numPhases());
-    // init does not resize surfacevol. Do it manually.
+    BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
+    // construction does not resize surfacevol. Do it manually.
     distributed_state.surfacevol().resize(grid.numCells()*state.numPhases(),
                                           std::numeric_limits<double>::max());
     BlackoilStateDataHandle state_handle(global_grid, grid,
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
index 4a393fc28..f0055b268 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
@@ -23,6 +23,11 @@
 
 #include "SimulatorFullyImplicitBlackoilOutput.hpp"
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/utility/DataMap.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/common/ErrorMacros.hpp>
@@ -32,9 +37,6 @@
 #include <opm/autodiff/GridHelpers.hpp>
 #include <opm/autodiff/BackupRestore.hpp>
 
-#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
-
-
 #include <sstream>
 #include <iomanip>
 #include <fstream>
@@ -52,7 +54,7 @@ namespace Opm
 
 
     void outputStateVtk(const UnstructuredGrid& grid,
-                        const SimulatorState& state,
+                        const SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir)
     {
@@ -89,16 +91,15 @@ namespace Opm
 
 
     void outputStateMatlab(const UnstructuredGrid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir)
     {
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
-        for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
-        {
-            const std::string& name = state.cellDataNames()[ i ];
+        for (const auto& pair : state.cellData()) {
+            const std::string& name = pair.first;
             std::string key;
             if( name == "SURFACEVOL" ) {
                 key = "surfvolume";
@@ -113,7 +114,7 @@ namespace Opm
                 continue;
             }
             // set data to datmap
-            dm[ key ] = &state.cellData()[ i ];
+            dm[ key ] = &pair.second;
         }
 
         std::vector<double> cell_velocity;
@@ -206,7 +207,7 @@ namespace Opm
 
 #ifdef HAVE_DUNE_CORNERPOINT
     void outputStateVtk(const Dune::CpGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir)
     {
@@ -296,7 +297,7 @@ namespace Opm
     void
     BlackoilOutputWriter::
     writeTimeStep(const SimulatorTimerInterface& timer,
-                  const SimulatorState& localState,
+                  const SimulationDataContainer& localState,
                   const WellState& localWellState,
                   bool substep)
     {
@@ -312,7 +313,7 @@ namespace Opm
             isIORank = parallelOutput_->collectToIORank( localState, localWellState, timer.reportStepNum() );
         }
 
-        const SimulatorState& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
+        const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
         const WellState& wellState  = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
 
         // serial output is only done on I/O rank
@@ -379,7 +380,35 @@ namespace Opm
                 }
                 catch ( const std::bad_cast& e )
                 {
+                    // store report step
+                    lastBackupReportStep_ = reportStep;
+                    // write resport step number
+                    backupfile_.write( (const char *) &reportStep, sizeof(int) );
 
+                    /*
+                    try {
+                        const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state);
+                        backupfile_ << boState;
+
+                        const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
+                        backupfile_ << boWellState;
+                    }
+                    catch ( const std::bad_cast& e )
+                    {
+
+                    }
+                    */
+
+                    /*
+                    const WellStateFullyImplicitBlackoil* boWellState =
+                        dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState);
+                    if( boWellState ) {
+                        backupfile_ << (*boWellState);
+                    }
+                    else
+                        OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed");
+                    */
+                    backupfile_ << std::flush;
                 }
 
                 /*
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
index e0208bff2..e1ed6915c 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
@@ -20,7 +20,6 @@
 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
 #define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
 #include <opm/core/grid.h>
-#include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/DataMap.hpp>
 #include <opm/common/ErrorMacros.hpp>
@@ -54,14 +53,17 @@
 namespace Opm
 {
 
+    class SimulationDataContainer;
+    class BlackoilState;
+
     void outputStateVtk(const UnstructuredGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir);
 
 
     void outputStateMatlab(const UnstructuredGrid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir);
 
@@ -70,23 +72,23 @@ namespace Opm
                                const std::string& output_dir);
 #ifdef HAVE_DUNE_CORNERPOINT
     void outputStateVtk(const Dune::CpGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir);
 #endif
 
     template<class Grid>
     void outputStateMatlab(const Grid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir)
     {
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
-        for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
+        for (const auto& pair : state.cellData()) 
         {
-            const std::string& name = state.cellDataNames()[ i ];
+            const std::string& name = pair.first;
             std::string key;
             if( name == "SURFACEVOL" ) {
                 key = "surfvolume";
@@ -101,7 +103,7 @@ namespace Opm
                 continue;
             }
             // set data to datmap
-            dm[ key ] = &state.cellData()[ i ];
+            dm[ key ] = &pair.second;
         }
 
         std::vector<double> cell_velocity;
@@ -155,7 +157,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& state,
+                           const SimulationDataContainer& state,
                            const WellState&,
                            bool /*substep*/ = false)
         {
@@ -185,7 +187,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& reservoirState,
+                           const SimulationDataContainer& reservoirState,
                            const WellState& wellState,
                            bool /*substep*/ = false)
         {
@@ -215,7 +217,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& reservoirState,
+                           const SimulationDataContainer& reservoirState,
                            const Opm::WellState& wellState,
                            bool substep = false);
 
@@ -242,7 +244,7 @@ namespace Opm
         void initFromRestartFile(const PhaseUsage& phaseusage,
                                  const double* permeability,
                                  const Grid& grid,
-                                 SimulatorState& simulatorstate,
+                                 SimulationDataContainer& simulatorstate,
                                  WellStateFullyImplicitBlackoil& wellstate);
 
         bool isRestart() const;
@@ -325,7 +327,7 @@ namespace Opm
     initFromRestartFile( const PhaseUsage& phaseusage,
                          const double* permeability,
                          const Grid& grid,
-                         SimulatorState& simulatorstate,
+                         SimulationDataContainer& simulatorstate,
                          WellStateFullyImplicitBlackoil& wellstate)
     {
         WellsManager wellsmanager(eclipseState_,
diff --git a/opm/autodiff/TransportSolverTwophaseAd.cpp b/opm/autodiff/TransportSolverTwophaseAd.cpp
index 65cc55b30..61ba63ab0 100644
--- a/opm/autodiff/TransportSolverTwophaseAd.cpp
+++ b/opm/autodiff/TransportSolverTwophaseAd.cpp
@@ -22,6 +22,7 @@
 #include <opm/autodiff/TransportSolverTwophaseAd.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/props/IncompPropertiesInterface.hpp>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/common/ErrorMacros.hpp>
diff --git a/opm/polymer/CompressibleTpfaPolymer.cpp b/opm/polymer/CompressibleTpfaPolymer.cpp
index 67feef107..604ba0334 100644
--- a/opm/polymer/CompressibleTpfaPolymer.cpp
+++ b/opm/polymer/CompressibleTpfaPolymer.cpp
@@ -88,8 +88,8 @@ namespace Opm
                                   PolymerBlackoilState& state,
                                   WellState& well_state)
     {
-        c_ = &state.concentration();
-        cmax_ = &state.maxconcentration();
+        c_ = &state.getCellData( state.CONCENTRATION );
+        cmax_ = &state.getCellData( state.CMAX );
         CompressibleTpfa::solve(dt, state, well_state);
     }
 
diff --git a/opm/polymer/IncompTpfaPolymer.cpp b/opm/polymer/IncompTpfaPolymer.cpp
index 138b12499..61905d87b 100644
--- a/opm/polymer/IncompTpfaPolymer.cpp
+++ b/opm/polymer/IncompTpfaPolymer.cpp
@@ -20,6 +20,8 @@
 
 #include <config.h>
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
 #include <opm/polymer/IncompTpfaPolymer.hpp>
 
 #include <opm/core/props/IncompPropertiesInterface.hpp>
@@ -99,8 +101,8 @@ namespace Opm
                                   PolymerState& state,
                                   WellState& well_state)
     {
-        c_ = &state.concentration();
-        cmax_ = &state.maxconcentration();
+        c_ = &state.getCellData( state.CONCENTRATION );
+        cmax_ = &state.getCellData( state.CMAX) ;
         if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) {
             solveRockComp(dt, state, well_state);
         } else {
@@ -116,7 +118,7 @@ namespace Opm
 
     /// Compute per-solve dynamic properties.
     void IncompTpfaPolymer::computePerSolveDynamicData(const double /*dt*/,
-                                                       const SimulatorState& state,
+                                                       const SimulationDataContainer& state,
                                                        const WellState& /*well_state*/)
     {
         // Computed here:
diff --git a/opm/polymer/IncompTpfaPolymer.hpp b/opm/polymer/IncompTpfaPolymer.hpp
index 8f61e5dc9..9feeb9bbe 100644
--- a/opm/polymer/IncompTpfaPolymer.hpp
+++ b/opm/polymer/IncompTpfaPolymer.hpp
@@ -37,6 +37,7 @@ namespace Opm
     class LinearSolverInterface;
     class PolymerState;
     class WellState;
+    class SimulationDataContainer;
 
     /// Encapsulating a tpfa pressure solver for the incompressible-fluid case with polymer.
     /// Supports gravity, wells controlled by bhp or reservoir rates,
@@ -96,7 +97,7 @@ namespace Opm
 
     private:
         virtual void computePerSolveDynamicData(const double dt,
-                                                const SimulatorState& state,
+                                                const SimulationDataContainer& state,
                                                 const WellState& well_state);
     private:
         // ------ Data that will remain unmodified after construction. ------
diff --git a/opm/polymer/PolymerBlackoilState.cpp b/opm/polymer/PolymerBlackoilState.cpp
new file mode 100644
index 000000000..e89d16c29
--- /dev/null
+++ b/opm/polymer/PolymerBlackoilState.cpp
@@ -0,0 +1,42 @@
+/*
+  Copyright 2012 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/>.
+*/
+
+
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/polymer/PolymerBlackoilState.hpp>
+
+namespace Opm
+{
+    const std::string PolymerBlackoilState::CONCENTRATION = "CONCENTRATION";
+    const std::string PolymerBlackoilState::CMAX = "CMAX";
+
+    PolymerBlackoilState::PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases) :
+        BlackoilState( number_of_cells , number_of_faces, num_phases)
+    {
+        registerCellData(CONCENTRATION , 1 , 0 );
+        registerCellData(CMAX , 1 , 0 );
+    }
+
+
+} // namespace Opm
+
+
+
+
diff --git a/opm/polymer/PolymerBlackoilState.hpp b/opm/polymer/PolymerBlackoilState.hpp
index 6fb1c36e5..521c284f4 100644
--- a/opm/polymer/PolymerBlackoilState.hpp
+++ b/opm/polymer/PolymerBlackoilState.hpp
@@ -33,27 +33,10 @@ namespace Opm
     class PolymerBlackoilState : public BlackoilState
     {
     public:
-        void init(const UnstructuredGrid& g, int num_phases)
-        {
-            this->init(g.number_of_cells, g.number_of_faces, num_phases);
-        }
+        static const std::string CONCENTRATION;
+        static const std::string CMAX;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            BlackoilState::init(number_of_cells, number_of_faces, num_phases);
-            concentration_.resize(number_of_cells, 0.0);
-            cmax_.resize(number_of_cells, 0.0);
-        }
-
-        std::vector<double>& concentration()    { return concentration_; }
-        std::vector<double>& maxconcentration() { return cmax_; }
-
-        const std::vector<double>& concentration() const    { return concentration_; }
-        const std::vector<double>& maxconcentration() const { return cmax_; }
-
-    private:
-        std::vector<double> concentration_;
-        std::vector<double> cmax_;
+        PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases);
     };
 
 } // namespace Opm
diff --git a/opm/polymer/PolymerState.cpp b/opm/polymer/PolymerState.cpp
new file mode 100644
index 000000000..83c323d56
--- /dev/null
+++ b/opm/polymer/PolymerState.cpp
@@ -0,0 +1,42 @@
+/*
+  Copyright 2012 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/>.
+*/
+
+
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/polymer/PolymerState.hpp>
+
+namespace Opm
+{
+    const std::string PolymerState::CONCENTRATION = "CONCENTRATION";
+    const std::string PolymerState::CMAX = "CMAX";
+
+    PolymerState::PolymerState(int number_of_cells, int number_of_faces, int num_phases) :
+        SimulationDataContainer( number_of_cells , number_of_faces , num_phases )
+    {
+        registerCellData(CONCENTRATION , 1 , 0 );
+        registerCellData(CMAX , 1 , 0 );
+    }
+
+
+} // namespace Opm
+
+
+
+
diff --git a/opm/polymer/PolymerState.hpp b/opm/polymer/PolymerState.hpp
index cf6cd5fd6..87f810a00 100644
--- a/opm/polymer/PolymerState.hpp
+++ b/opm/polymer/PolymerState.hpp
@@ -20,8 +20,8 @@
 #ifndef OPM_POLYMERSTATE_HEADER_INCLUDED
 #define OPM_POLYMERSTATE_HEADER_INCLUDED
 
+#include <opm/common/data/SimulationDataContainer.hpp>
 
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/grid.h>
 #include <vector>
 
@@ -29,34 +29,13 @@ namespace Opm
 {
 
     /// Simulator state for a two-phase simulator with polymer.
-    class PolymerState : public SimulatorState
+    class PolymerState : public SimulationDataContainer
     {
     public:
+        static const std::string CONCENTRATION;
+        static const std::string CMAX;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            SimulatorState::init( number_of_cells , number_of_faces , num_phases );
-            registerCellData("CONCENTRATION" , 1 , 0 );
-            registerCellData("CMAX" , 1 , 0 );
-        }
-
-        const std::vector<double>& concentration() const {
-            return getCellData("CONCENTRATION");
-        }
-
-        const std::vector<double>& maxconcentration() const {
-            return getCellData("CMAX");
-        }
-
-
-        std::vector<double>& concentration()  {
-            return getCellData("CONCENTRATION");
-        }
-
-        std::vector<double>& maxconcentration()  {
-            return getCellData("CMAX");
-        }
-
+        PolymerState(int number_of_cells, int number_of_faces, int num_phases);
     };
 
 } // namespace Opm
diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp
index a483cda58..b20e91492 100644
--- a/opm/polymer/SimulatorCompressiblePolymer.cpp
+++ b/opm/polymer/SimulatorCompressiblePolymer.cpp
@@ -268,7 +268,7 @@ namespace Opm
         total_timer.start();
         double init_surfvol[2] = { 0.0 };
         double inplace_surfvol[2] = { 0.0 };
-        double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
         double polymass_adsorbed = computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_);
         double init_polymass = polymass + polymass_adsorbed;
         double tot_injected[2] = { 0.0 };
@@ -301,7 +301,7 @@ namespace Opm
         if (check_well_controls_) {
             computeFractionalFlow(props_, poly_props_, allcells_,
                                   state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
-                                  state.concentration(), state.maxconcentration(),
+                                  state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ) ,
                                   fractional_flows);
             wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
         }
@@ -397,7 +397,7 @@ namespace Opm
                            state.pressure(), state.temperature(), &initial_porevol[0], &porevol[0],
                            &transport_src[0], &polymer_inflow_c[0], stepsize,
                            state.saturation(), state.surfacevol(),
-                           state.concentration(), state.maxconcentration());
+                           state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             double substep_injected[2] = { 0.0 };
             double substep_produced[2] = { 0.0 };
             double substep_polyinj = 0.0;
@@ -416,7 +416,7 @@ namespace Opm
             if (gravity_ != 0 && use_segregation_split_) {
                 tsolver_.solveGravity(columns_, stepsize,
                                       state.saturation(), state.surfacevol(),
-                                      state.concentration(), state.maxconcentration());
+                                      state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             }
         }
         transport_timer.stop();
@@ -426,7 +426,7 @@ namespace Opm
 
         // Report volume balances.
         Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
-        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
         polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
                                                         state, rock_comp_props_);
         tot_injected[0] += injected[0];
@@ -539,8 +539,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             dm["surfvol"] = &state.surfacevol();
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@@ -556,8 +556,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             dm["surfvol"] = &state.surfacevol();
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp
index 615bbd4eb..973b61e23 100644
--- a/opm/polymer/SimulatorPolymer.cpp
+++ b/opm/polymer/SimulatorPolymer.cpp
@@ -292,8 +292,8 @@ namespace Opm
         total_timer.start();
         double init_satvol[2] = { 0.0 };
         double satvol[2] = { 0.0 };
-        double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
-        double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
+        double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
+        double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
         double init_polymass = polymass + polymass_adsorbed;
         double injected[2] = { 0.0 };
         double produced[2] = { 0.0 };
@@ -330,7 +330,7 @@ namespace Opm
         // Solve pressure.
         if (check_well_controls_) {
             computeFractionalFlow(props_, poly_props_, allcells_,
-                                  state.saturation(), state.concentration(), state.maxconcentration(),
+                                  state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ),
                                   fractional_flows);
             wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
         }
@@ -425,7 +425,7 @@ namespace Opm
         injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
         for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
             tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
-                           state.saturation(), state.concentration(), state.maxconcentration());
+                           state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             Opm::computeInjectedProduced(props_, poly_props_,
                                          state,
                                          transport_src, polymer_inflow_c, stepsize,
@@ -438,7 +438,7 @@ namespace Opm
             polyprod += substep_polyprod;
             if (use_segregation_split_) {
                 tsolver_.solveGravity(columns_, &porevol[0], stepsize,
-                                      state.saturation(), state.concentration(), state.maxconcentration());
+                                      state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             }
         }
         transport_timer.stop();
@@ -448,8 +448,8 @@ namespace Opm
 
         // Report volume balances.
         Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
-        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
-        polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
+        polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
         tot_injected[0] += injected[0];
         tot_injected[1] += injected[1];
         tot_produced[0] += produced[0];
@@ -559,8 +559,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;
@@ -575,8 +575,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;
@@ -611,8 +611,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
index 83dead022..485c48e14 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
@@ -125,8 +125,10 @@ namespace Opm {
                 WellState& well_state)
     {
         Base::prepareStep(dt, reservoir_state, well_state);
+        auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
         // Initial max concentration of this time step from PolymerBlackoilState.
-        cmax_ = Eigen::Map<const V>(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_));
+
+        cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
     }
 
 
@@ -168,9 +170,10 @@ namespace Opm {
 
         // Initial polymer concentration.
         if (has_polymer_) {
-            assert (not x.concentration().empty());
-            const int nc = x.concentration().size();
-            const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
+            const auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert (not concentration.empty());
+            const int nc = concentration.size();
+            const V c = Eigen::Map<const V>(concentration.data() , nc);
             // Concentration belongs after other reservoir vars but before well vars.
             auto concentration_pos = vars0.begin() + fluid_.numPhases();
             assert(concentration_pos == vars0.end() - 2);
@@ -255,12 +258,14 @@ namespace Opm {
     template <class Grid>
     void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
     {
-        const int nc = AutoDiffGrid::numCells(grid_);
-        V tmp = V::Zero(nc);
-        for (int i = 0; i < nc; ++i) {
-            tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
-        }
-        std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
+        auto& max_concentration = state.getCellData( state.CMAX );
+        const auto& concentration = state.getCellData( state.CONCENTRATION );
+        std::transform( max_concentration.begin() ,
+                        max_concentration.end() ,
+                        concentration.begin() ,
+                        max_concentration.begin() ,
+                        [](double c_max , double c) { return std::max( c_max , c ); });
+
     }
 
 
@@ -397,10 +402,13 @@ namespace Opm {
             // Call base version.
             Base::updateState(modified_dx, reservoir_state, well_state);
 
-            // Update concentration.
-            const V c_old = Eigen::Map<const V>(&reservoir_state.concentration()[0], nc, 1);
-            const V c = (c_old - dc).max(zero);
-            std::copy(&c[0], &c[0] + nc, reservoir_state.concentration().begin());
+            {
+                auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
+                // Update concentration.
+                const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
+                const V c = (c_old - dc).max(zero);
+                std::copy(&c[0], &c[0] + nc, concentration.begin());
+            }
         } else {
             // Just forward call to base version.
             Base::updateState(dx, reservoir_state, well_state);
diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
index 1ca4df5dd..6e9e195e7 100644
--- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
+++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
@@ -206,7 +206,7 @@ namespace {
         const std::vector<double>& polymer_inflow = xw.polymerInflow();
 
         // Initial max concentration of this time step from PolymerBlackoilState.
-        cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
+        cmax_ = Eigen::Map<V>(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_));
 
         const SolutionState state = constantState(x, xw);
         computeAccum(state, 0);
@@ -366,9 +366,18 @@ namespace {
         state.saturation[1] = ADB::constant(so, bpat);
 
         // Concentration
-        assert(not x.concentration().empty());
-        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
-        state.concentration = ADB::constant(c);
+        {
+            auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert(concentration.empty());
+            const V c = Eigen::Map<const V>(concentration.data(), nc);
+            // Do not understand:
+            //concentration = ADB::constant(c);
+            // Old code based on concentraton() method had the statement:
+            //
+            //   state.concentration = ADB::constant(c)
+            //
+            // This looks like it was a method assignment - how did it even compile?
+        }
 
         // Well rates.
         assert (not xw.wellRates().empty());
@@ -413,9 +422,12 @@ namespace {
         vars0.push_back(sw0);
 
         // Initial concentration.
-        assert (not x.concentration().empty());
-        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
-        vars0.push_back(c);
+        {
+            auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert (not concentration.empty());
+            const V c = Eigen::Map<const V>(concentration.data() , nc);
+            vars0.push_back(c);
+        }
 
         // Initial well rates.
         assert (not xw.wellRates().empty());
@@ -511,11 +523,14 @@ namespace {
     {
         const int nc = grid_.number_of_cells;
         V tmp = V::Zero(nc);
+        const auto& concentration = state.getCellData( state.CONCENTRATION );
+        auto& cmax = state.getCellData( state.CMAX );
+
         for (int i = 0; i < nc; ++i) {
-            tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
+            tmp[i] = std::max(cmax[i], concentration[i]);
         }
 
-        std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
+        std::copy(&tmp[0], &tmp[0] + nc, cmax.begin());
     }
 
 
@@ -764,11 +779,14 @@ namespace {
         // Concentration updates.
 //        const double dcmax = 0.3 * polymer_props_ad_.cMax();
 //		std::cout << "\n the max concentration: " << dcmax / 0.3 << std::endl;
-        const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc, 1);
-//        const V dc_limited = sign(dc) * dc.abs().min(dcmax);
-//        const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
-        const V c = (c_old - dc).max(zero);
-        std::copy(&c[0], &c[0] + nc, state.concentration().begin());
+        {
+            auto& concentration = state.getCellData( state.CONCENTRATION );
+            const V c_old = Eigen::Map<const V>(concentration.data() , nc, 1);
+            //        const V dc_limited = sign(dc) * dc.abs().min(dcmax);
+            //        const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
+            const V c = (c_old - dc).max(zero);
+            std::copy(&c[0], &c[0] + nc, concentration.begin());
+        }
 
         // Qs update.
         // Since we need to update the wellrates, that are ordered by wells,
diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp
index 7fa369fa3..971066f5f 100644
--- a/opm/polymer/polymerUtilities.cpp
+++ b/opm/polymer/polymerUtilities.cpp
@@ -211,8 +211,8 @@ namespace Opm
             OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
         }
         const std::vector<double>& s = state.saturation();
-        const std::vector<double>& c = state.concentration();
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
         std::fill(injected, injected + np, 0.0);
         std::fill(produced, produced + np, 0.0);
         polyinj = 0.0;
@@ -282,8 +282,8 @@ namespace Opm
         const std::vector<double>& temp = state.temperature();
         const std::vector<double>& s = state.saturation();
         const std::vector<double>& z = state.surfacevol();
-        const std::vector<double>& c = state.concentration();
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
         std::fill(injected, injected + np, 0.0);
         std::fill(produced, produced + np, 0.0);
         polyinj = 0.0;
@@ -406,7 +406,7 @@ namespace Opm
             porosity.assign(props.porosity(), props.porosity() + num_cells);
         }
         double abs_mass = 0.0;
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
 	for (int cell = 0; cell < num_cells; ++cell) {
             double c_ads;
             polyprops.simpleAdsorption(cmax[cell], c_ads);
diff --git a/tests/test_rateconverter.cpp b/tests/test_rateconverter.cpp
index 6c9e28898..15ec67aee 100644
--- a/tests/test_rateconverter.cpp
+++ b/tests/test_rateconverter.cpp
@@ -32,6 +32,7 @@
 
 #include <boost/test/unit_test.hpp>
 
+#include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
 #include <opm/core/utility/Units.hpp>
@@ -107,8 +108,7 @@ BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture<SetupSimple>)
     Region reg{ 0 };
     RCvrt  cvrt(ad_props, reg);
 
-    Opm::BlackoilState x;
-    x.init(*grid.c_grid(), 3);
+    Opm::BlackoilState x( Opm::UgGridHelpers::numCells( *grid.c_grid()) , Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3);
 
     cvrt.defineState(x);
 
diff --git a/tests/test_singlecellsolves.cpp b/tests/test_singlecellsolves.cpp
index 25f298d85..1a9ba9800 100644
--- a/tests/test_singlecellsolves.cpp
+++ b/tests/test_singlecellsolves.cpp
@@ -71,7 +71,7 @@ try
 
     boost::scoped_ptr<GridManager> grid;
     boost::scoped_ptr<IncompPropertiesInterface> props;
-    PolymerState state;
+    std::unique_ptr<PolymerState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -81,23 +81,35 @@ try
     // Grid init.
     grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0));
     // Rock and fluid init.
-    props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-    // Init state variables (saturation and pressure).
-    initStateBasic(*grid->c_grid(), *props, param, 0.0, state);
-    // Init Polymer state
-    if (param.has("poly_init")) {
-        double poly_init = param.getDefault("poly_init", 0.0);
-        for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-            double smin[2], smax[2];
-            props->satRange(1, &cell, smin, smax);
-            if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                state.concentration()[cell] = poly_init;
-                state.maxconcentration()[cell] = poly_init;
-            } else {
-                state.saturation()[2*cell + 0] = 0.;
-                state.saturation()[2*cell + 1] = 1.;
-                state.concentration()[cell] = 0.;
-                state.maxconcentration()[cell] = 0.;
+    {
+        const UnstructuredGrid& ug_grid = *(grid->c_grid());
+        state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+
+        props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+        // Init state variables (saturation and pressure).
+        initStateBasic(*grid->c_grid(), *props, param, 0.0, *state);
+        // Init Polymer state
+        if (param.has("poly_init")) {
+            double poly_init = param.getDefault("poly_init", 0.0);
+            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
+                double smin[2], smax[2];
+                props->satRange(1, &cell, smin, smax);
+
+                auto& saturation = state->saturation();
+                auto& concentration = state->getCellData( state->CONCENTRATION );
+                auto& max_concentration = state->getCellData( state->CMAX );
+
+                props->satRange(1, &cell, smin, smax);
+                if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                    concentration[cell] = poly_init;
+                    max_concentration[cell] = poly_init;
+                } else {
+                    saturation[2*cell + 0] = 0.;
+                    saturation[2*cell + 1] = 1.;
+                    concentration[cell] = 0.;
+                    max_concentration[cell] = 0.;
+                }
+
             }
         }
     }
@@ -210,7 +222,7 @@ try
     if (face01 == -1) {
         OPM_THROW(std::runtime_error, "Could not find face adjacent to cells [0 1]");
     }
-    state.faceflux()[face01] = src[0];
+    state->faceflux()[face01] = src[0];
     for (int sats = 0; sats < num_sats; ++sats) {
         const double s = double(sats)/double(num_sats - 1);
         const double ff = s; // Simplified a lot...
@@ -220,20 +232,26 @@ try
             // std::cout << "(s, c) = (" << s << ", " << c << ")\n";
             transport_src[0] = src[0]*ff;
             // Resetting the state for next run.
-            state.saturation()[0] = 0.0;
-            state.saturation()[1] = 0.0;
-            state.concentration()[0] = 0.0;
-            state.concentration()[1] = 0.0;
-            state.maxconcentration()[0] = 0.0;
-            state.maxconcentration()[1] = 0.0;
-            reorder_model.solve(&state.faceflux()[0],
+            auto& saturation = state->saturation();
+            auto& concentration = state->getCellData( state->CONCENTRATION );
+            auto& max_concentration = state->getCellData( state->CMAX );
+
+
+            saturation[0] = 0.0;
+            saturation[1] = 0.0;
+            concentration[0] = 0.0;
+            concentration[1] = 0.0;
+            max_concentration[0] = 0.0;
+            max_concentration[1] = 0.0;
+
+            reorder_model.solve(&state->faceflux()[0],
                                 &porevol[0],
                                 &transport_src[0],
                                 &polymer_inflow_c[0],
                                 dt,
-                                state.saturation(),
-                                state.concentration(),
-                                state.maxconcentration());
+                                saturation,
+                                concentration,
+                                max_concentration);
 
 #ifdef PROFILING
             // Extract residual counts.