diff --git a/opm/core/InjectionSpecification.cpp b/opm/core/InjectionSpecification.cpp
index 691830c6c..dca7fc76b 100644
--- a/opm/core/InjectionSpecification.cpp
+++ b/opm/core/InjectionSpecification.cpp
@@ -8,7 +8,8 @@ namespace Opm
       surface_flow_max_rate_(-1e100),
       reservoir_flow_max_rate_(-1e100),
       BHP_limit_(-1e100),
-      reinjection_fraction_target_(-1e100),
+      reinjection_fraction_target_(1),
+      voidage_replacment_fraction_(1),
       guide_rate_(1.0),
       guide_rate_type_(NONE_GRT)
     {
diff --git a/opm/core/InjectionSpecification.hpp b/opm/core/InjectionSpecification.hpp
index 62c9e6678..107296c66 100644
--- a/opm/core/InjectionSpecification.hpp
+++ b/opm/core/InjectionSpecification.hpp
@@ -32,6 +32,7 @@ namespace Opm
         double reservoir_flow_max_rate_;
         double BHP_limit_;
         double reinjection_fraction_target_;
+        double voidage_replacment_fraction_;
         double guide_rate_;
         GuideRateType guide_rate_type_;
     };
diff --git a/opm/core/ProductionSpecification.hpp b/opm/core/ProductionSpecification.hpp
index ff5ace8c8..dd5bc7c09 100644
--- a/opm/core/ProductionSpecification.hpp
+++ b/opm/core/ProductionSpecification.hpp
@@ -11,7 +11,7 @@ namespace Opm
 
         enum ControlMode
         {
-            NONE, ORAT, WRAT, GRAT, LRAT, CRAT, RESV, PRBL, BHP, THP, GRUP, FLD
+            NONE = 0, ORAT = 1, WRAT=2, GRAT=3, LRAT=4, CRAT=5, RESV=6, PRBL=7, BHP=8, THP=9, GRUP=10, FLD=11
         };
 
         enum Procedure
diff --git a/opm/core/WellCollection.cpp b/opm/core/WellCollection.cpp
index f873c259f..585c5eaaf 100644
--- a/opm/core/WellCollection.cpp
+++ b/opm/core/WellCollection.cpp
@@ -122,4 +122,22 @@ namespace Opm
             roots_[i]->applyInjGroupControls();
         }
     }
+
+    /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+    /// \param[in]    well_reservoirrates_phase
+    ///                         A vector containing reservoir rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    /// \param[in]    well_surfacerates_phase
+    ///                         A vector containing surface rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+
+    void WellCollection::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                          const std::vector<double>& well_surfacerates_phase)
+    {
+        for (size_t i = 0; i < roots_.size(); ++i) {
+            roots_[i]->applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
+        }
+    }
 }
diff --git a/opm/core/WellCollection.hpp b/opm/core/WellCollection.hpp
index a5db76d71..3e71ab020 100644
--- a/opm/core/WellCollection.hpp
+++ b/opm/core/WellCollection.hpp
@@ -88,6 +88,18 @@ namespace Opm
         /// Applies all group controls (injection and production)
         void applyGroupControls();
         
+        /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+        /// \param[in]    well_reservoirrates_phase
+        ///                         A vector containing reservoir rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[in]    well_surfacerates_phase
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                              const std::vector<double>& well_surfacerates_phase);
+        
     private:
         // To account for the possibility of a forest
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;
diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp
index 931087a55..006b7610f 100644
--- a/opm/core/WellsGroup.cpp
+++ b/opm/core/WellsGroup.cpp
@@ -453,7 +453,7 @@ namespace Opm
                 // as that would check if we're under group control, something we're not.
                 const double children_guide_rate = children_[i]->productionGuideRate(true);
                 children_[i]->applyProdGroupControl(prod_mode, 
-                                                    (children_guide_rate / my_guide_rate) * getTarget(prod_mode), 
+                                                   (children_guide_rate / my_guide_rate) * getTarget(prod_mode), 
                                                     false);
             }
             break;
@@ -489,8 +489,9 @@ namespace Opm
             }
             return;
         }
+        case InjectionSpecification::VREP:
         case InjectionSpecification::REIN:
-            std::cout << "WARNING: Ignoring control type REIN" << std::endl;
+            std::cout << "Replacement keywords found, remember to call applyExplicitReinjectionControls." << std::endl;
             return;
         case InjectionSpecification::FLD:
         case InjectionSpecification::NONE:
@@ -528,7 +529,90 @@ namespace Opm
         return sum;
     }
 
+    /// Gets the total production flow of the given phase. 
+    /// \param[in] phase_flows      A vector containing rates by phase for each well.
+    ///                             Is assumed to be ordered the same way as the related Wells-struct,
+    ///                             with all phase rates of a single well adjacent in the array.
+    /// \param[in] phase            The phase for which to sum up.
 
+    double WellsGroup::getTotalProductionFlow(const std::vector<double>& phase_flows,
+                                             const BlackoilPhases::PhaseIndex phase)
+    {
+        double sum = 0.0;
+        for (size_t i = 0; i < children_.size(); ++i) {
+            sum += children_[i]->getTotalProductionFlow(phase_flows, phase);
+        }
+        return sum;
+    }
+
+    /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+    /// \param[in]    well_reservoirrates_phase
+    ///                         A vector containing reservoir rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    /// \param[in]    well_surfacerates_phase
+    ///                         A vector containing surface rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    void WellsGroup::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                      const std::vector<double>& well_surfacerates_phase)
+    {
+        if (injSpec().control_mode_ == InjectionSpecification::REIN) {
+            BlackoilPhases::PhaseIndex phase;
+            switch (injSpec().injector_type_) {
+            case InjectionSpecification::WATER:
+                phase = BlackoilPhases::Aqua;
+                break;
+            case InjectionSpecification::GAS:
+                phase = BlackoilPhases::Vapour;
+                break;
+            case InjectionSpecification::OIL:
+                phase = BlackoilPhases::Liquid;
+                break;
+            }
+            const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
+            const double my_guide_rate = injectionGuideRate(true);
+            for (size_t i = 0; i < children_.size(); ++i) {
+                // Apply for all children. 
+                // Note, we do _not_ want to call the applyProdGroupControl in this object,
+                // as that would check if we're under group control, something we're not.
+                const double children_guide_rate = children_[i]->injectionGuideRate(true);
+#ifdef DIRTY_WELLCTRL_HACK
+                children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
+                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
+                        false);
+#else
+                children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
+                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
+                        false);
+#endif
+            }
+        }
+        else if (injSpec().control_mode_ == InjectionSpecification::VREP) {
+            double total_produced = 0.0;
+            if (phaseUsage().phase_used[BlackoilPhases::Aqua]) {
+                total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Aqua);
+            }
+            if (phaseUsage().phase_used[BlackoilPhases::Liquid]) {
+                total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Liquid);
+            }
+            if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
+                total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
+            }
+            
+                  const double my_guide_rate = injectionGuideRate(true);
+            for (size_t i = 0; i < children_.size(); ++i) {
+                // Apply for all children. 
+                // Note, we do _not_ want to call the applyProdGroupControl in this object,
+                // as that would check if we're under group control, something we're not.
+                const double children_guide_rate = children_[i]->injectionGuideRate(true);
+                children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
+                        (children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
+                        false);
+            }
+            
+        }
+    }
 
     // ==============    WellNode members   ============
 
@@ -691,7 +775,7 @@ namespace Opm
     {
         // Not changing if we're not forced to change
         if (!forced 
-             && (injSpec().control_mode_ != InjectionSpecification::GRUP || injSpec().control_mode_ != InjectionSpecification::NONE)) {
+             && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) {
             return;
         }
         if (!wells_->type[self_index_] == INJECTOR) {
@@ -728,13 +812,47 @@ namespace Opm
     }
 
 
+    /// Gets the total production flow of the given phase. 
+    /// \param[in] phase_flows      A vector containing rates by phase for each well.
+    ///                             Is assumed to be ordered the same way as the related Wells-struct,
+    ///                             with all phase rates of a single well adjacent in the array.
+    /// \param[in] phase            The phase for which to sum up.
+
+    double WellNode::getTotalProductionFlow(const std::vector<double>& phase_flows,
+                                            const BlackoilPhases::PhaseIndex phase)
+    {
+        if (type() == INJECTOR) {
+            return 0.0;
+        }
+        return phase_flows[self_index_*phaseUsage().num_phases + phaseUsage().phase_pos[phase]];
+    }
+    
+    WellType WellNode::type() const {
+        return wells_->type[self_index_];
+    }
+
+    /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+    /// \param[in]    well_reservoirrates_phase
+    ///                         A vector containing reservoir rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    /// \param[in]    well_surfacerates_phase
+    ///                         A vector containing surface rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    void WellNode::applyExplicitReinjectionControls(const std::vector<double>&,
+                                                    const std::vector<double>&)
+    {
+        // Do nothing at well level.
+    }
     void WellNode::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
                                          const double target,
                                          const bool forced)
     {
         // Not changing if we're not forced to change
         if (!forced && (prodSpec().control_mode_ != ProductionSpecification::GRUP
-                        || prodSpec().control_mode_ != ProductionSpecification::NONE)) {
+                        && prodSpec().control_mode_ != ProductionSpecification::NONE)) {
+            std::cout << "Returning" << std::endl;
             return;
         }
         if (!wells_->type[self_index_] == PRODUCER) {
@@ -771,6 +889,7 @@ namespace Opm
             distr[phase_pos[BlackoilPhases::Vapour]] = 1.0;
             break;
         case ProductionSpecification::LRAT:
+            std::cout << "applying rate" << std::endl;  
             wct = SURFACE_RATE;
             if (!phase_used[BlackoilPhases::Liquid]) {
                 THROW("Oil phase not active and LRAT control specified.");
@@ -981,6 +1100,8 @@ namespace Opm
                         injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
                         injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
                         injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_;
+                        injection_specification.reinjection_fraction_target_ = line.reinjection_fraction_target_;
+                        injection_specification.voidage_replacment_fraction_ = line.voidage_replacement_fraction_;
                     }
                 }
             }
diff --git a/opm/core/WellsGroup.hpp b/opm/core/WellsGroup.hpp
index 6fe5508bc..3c14950dd 100644
--- a/opm/core/WellsGroup.hpp
+++ b/opm/core/WellsGroup.hpp
@@ -179,6 +179,27 @@ namespace Opm
         ///                       wells under group control
         virtual double injectionGuideRate(bool only_group) = 0;
         
+        /// Gets the total production flow of the given phase. 
+        /// \param[in] phase_flows      A vector containing rates by phase for each well.
+        ///                             Is assumed to be ordered the same way as the related Wells-struct,
+        ///                             with all phase rates of a single well adjacent in the array.
+        /// \param[in] phase            The phase for which to sum up.
+        virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
+                                              const BlackoilPhases::PhaseIndex phase) = 0;
+        
+        /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+        /// \param[in]    well_reservoirrates_phase
+        ///                         A vector containing reservoir rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[in]    well_surfacerates_phase
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                      const std::vector<double>& well_surfacerates_phase) = 0;
+        
+        
     protected:
         /// Calculates the correct rate for the given ProductionSpecification::ControlMode
         double rateByMode(const double* res_rates, 
@@ -258,6 +279,26 @@ namespace Opm
         /// \param[in] only_group If true, will only accumelate guide rates for 
         ///                       wells under group control
         virtual double injectionGuideRate(bool only_group);
+        
+        /// Gets the total production flow of the given phase. 
+        /// \param[in] phase_flows      A vector containing rates by phase for each well.
+        ///                             Is assumed to be ordered the same way as the related Wells-struct,
+        ///                             with all phase rates of a single well adjacent in the array.
+        /// \param[in] phase            The phase for which to sum up.
+        virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
+                                              const BlackoilPhases::PhaseIndex phase);
+        
+        /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+        /// \param[in]    well_reservoirrates_phase
+        ///                         A vector containing reservoir rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[in]    well_surfacerates_phase
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                      const std::vector<double>& well_surfacerates_phase);
 
     private:
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > children_;
@@ -327,6 +368,29 @@ namespace Opm
         /// \param[in] only_group If true, will only accumelate guide rates for 
         ///                       wells under group control
         virtual double injectionGuideRate(bool only_group);
+        
+        /// Gets the total production flow of the given phase. 
+        /// \param[in] phase_flows      A vector containing rates by phase for each well.
+        ///                             Is assumed to be ordered the same way as the related Wells-struct,
+        ///                             with all phase rates of a single well adjacent in the array.
+        /// \param[in] phase            The phase for which to sum up.
+        virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
+                                              const BlackoilPhases::PhaseIndex phase);
+
+        /// Returns the type of the well.
+        WellType type() const;
+        
+        /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+        /// \param[in]    well_reservoirrates_phase
+        ///                         A vector containing reservoir rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[in]    well_surfacerates_phase
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                      const std::vector<double>& well_surfacerates_phase);
 
     private:
         Wells* wells_;
diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp
index 8c04f9135..8b041be0b 100644
--- a/opm/core/WellsManager.cpp
+++ b/opm/core/WellsManager.cpp
@@ -702,4 +702,20 @@ namespace Opm
                                               well_surfacerates_phase);
     }
 
+    /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+    /// \param[in]    well_reservoirrates_phase
+    ///                         A vector containing reservoir rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+    /// \param[in]    well_surfacerates_phase
+    ///                         A vector containing surface rates by phase for each well.
+    ///                         Is assumed to be ordered the same way as the related Wells-struct,
+    ///                         with all phase rates of a single well adjacent in the array.
+
+    void WellsManager::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                                        const std::vector<double>& well_surfacerates_phase)
+    {
+        well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
+    }
+
 } // namespace Opm
diff --git a/opm/core/WellsManager.hpp b/opm/core/WellsManager.hpp
index 32b04e034..3fbdbd1a1 100644
--- a/opm/core/WellsManager.hpp
+++ b/opm/core/WellsManager.hpp
@@ -88,6 +88,18 @@ namespace Opm
         bool conditionsMet(const std::vector<double>& well_bhp,
                            const std::vector<double>& well_reservoirrates_phase,
                            const std::vector<double>& well_surfacerates_phase);
+        
+        /// Applies explicit reinjection controls. This must be called at each timestep to be correct.
+        /// \param[in]    well_reservoirrates_phase
+        ///                         A vector containing reservoir rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[in]    well_surfacerates_phase
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
+                                              const std::vector<double>& well_surfacerates_phase);
 
     private:
 	// Disable copying and assignment.
diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp
index cd0a07b6b..c7222b947 100644
--- a/opm/core/pressure/IncompTpfa.cpp
+++ b/opm/core/pressure/IncompTpfa.cpp
@@ -26,6 +26,7 @@
 #include <opm/core/linalg/sparse_sys.h>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/newwells.h>
+#include <string.h>
 
 namespace Opm
 {
@@ -210,13 +211,16 @@ namespace Opm
         F.bc = bcs;
         F.totmob = &totmob[0];
         F.wdp = &wdp[0];
-        
+        int ok = true;
         if (rock_comp.empty()) {
-            ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
+            ok = ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
         } else {
-            ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
+            ok = ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
                                        &porevol[0], &rock_comp[0], dt, &pressure[0], h_);
         }
+        if (!ok) {
+            THROW("Failed assembling pressure system.");
+        }
 
 	linsolver_.solve(h_->A, h_->b, h_->x);
 
@@ -233,9 +237,108 @@ namespace Opm
             soln.well_flux = &well_rate[0];
             soln.well_press = &well_bhp[0];
         }
+
         ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
     }
 
+    
+    void IncompTpfa::solveIncrement(const std::vector<double>& totmob,
+				    const std::vector<double>& omega,
+				    const std::vector<double>& src,
+				    const std::vector<double>& wdp,
+				    const FlowBoundaryConditions* bcs,
+				    const std::vector<double>& porevol,
+				    const std::vector<double>& rock_comp,
+				    const std::vector<double>& prev_pressure,
+				    const std::vector<double>& initial_porevol,
+				    const double dt,
+				    std::vector<double>& pressure_increment)
 
+    {
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+	if (!omega.empty()) {
+	    if (gpress_.empty()) {
+		THROW("Nozero omega argument given, but gravity was null in constructor.");
+	    }
+	    mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+				  &omega[0],
+				  &gpress_[0], &gpress_omegaweighted_[0]);
+	} else {
+	    if (!gpress_.empty()) {
+		THROW("Empty omega argument given, but gravity was non-null in constructor.");
+	    }
+	}
+
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
+        if (! src.empty()) { F.src = &src[0]; }
+        F.bc = bcs;
+        F.totmob = &totmob[0];
+        F.wdp = &wdp[0];
+        
+        if (rock_comp.empty()) {
+            ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
+        } else {
+            ifs_tpfa_assemble_comprock_increment(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
+                                       &porevol[0], &rock_comp[0], dt, &prev_pressure[0], 
+				       &initial_porevol[0], h_);
+        }
+
+	linsolver_.solve(h_->A, h_->b, &pressure_increment[0]);
+
+    }
+
+    void IncompTpfa::computeFaceFlux(const std::vector<double>& totmob,
+				     const std::vector<double>& omega,
+				     const std::vector<double>& src,
+				     const std::vector<double>& wdp,
+				     const FlowBoundaryConditions* bcs,
+				     std::vector<double>& pressure,
+				     std::vector<double>& faceflux,
+				     std::vector<double>& well_bhp,
+				     std::vector<double>& well_rate) {
+
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+	if (!omega.empty()) {
+	    if (gpress_.empty()) {
+		THROW("Nozero omega argument given, but gravity was null in constructor.");
+	    }
+	    mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+				  &omega[0],
+				  &gpress_[0], &gpress_omegaweighted_[0]);
+	} else {
+	    if (!gpress_.empty()) {
+		THROW("Empty omega argument given, but gravity was non-null in constructor.");
+	    }
+	}
+
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
+        if (! src.empty()) { F.src = &src[0]; }
+        F.bc = bcs;
+        F.totmob = &totmob[0];
+        F.wdp = &wdp[0];
+
+	ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
+	
+	faceflux.resize(grid_.number_of_faces);
+
+        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
+        soln.cell_press = &pressure[0];
+        soln.face_flux  = &faceflux[0];
+
+        if(wells_ != NULL) {
+            well_bhp.resize(wells_->number_of_wells);
+            well_rate.resize(wells_->well_connpos[ wells_->number_of_wells ]);
+            soln.well_flux = &well_rate[0];
+            soln.well_press = &well_bhp[0];
+        }
+
+	memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
+
+        ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
+    }
 
 } // namespace Opm
diff --git a/opm/core/pressure/IncompTpfa.hpp b/opm/core/pressure/IncompTpfa.hpp
index f40d3e232..f14d9108f 100644
--- a/opm/core/pressure/IncompTpfa.hpp
+++ b/opm/core/pressure/IncompTpfa.hpp
@@ -121,6 +121,28 @@ namespace Opm
                    std::vector<double>& well_bhp,
                    std::vector<double>& well_rate);
 
+        void solveIncrement(const std::vector<double>& totmob,
+			    const std::vector<double>& omega,
+			    const std::vector<double>& src,
+			    const std::vector<double>& wdp,
+			    const FlowBoundaryConditions* bcs,
+			    const std::vector<double>& porevol,
+			    const std::vector<double>& rock_comp,
+			    const std::vector<double>& prev_pressure,
+			    const std::vector<double>& initial_porevol,
+			    const double dt,
+			    std::vector<double>& pressure_increment);
+
+	void computeFaceFlux(const std::vector<double>& totmob,
+					 const std::vector<double>& omega,
+					 const std::vector<double>& src,
+					 const std::vector<double>& wdp,
+					 const FlowBoundaryConditions* bcs,
+					 std::vector<double>& pressure,
+					 std::vector<double>& faceflux,
+					 std::vector<double>& well_bhp,
+					 std::vector<double>& well_rate);
+
         /// Expose read-only reference to internal half-transmissibility.
         const ::std::vector<double>& getHalfTrans() const { return htrans_; }
 
diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c
index c3566c931..030bae422 100644
--- a/opm/core/pressure/tpfa/ifs_tpfa.c
+++ b/opm/core/pressure/tpfa/ifs_tpfa.c
@@ -11,6 +11,16 @@
 #include <opm/core/pressure/tpfa/ifs_tpfa.h>
 
 
+void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v) {
+    size_t i,j;
+    for (j = 0; j < A->m; ++j) {
+        v[j] = 0;
+        for (i = (size_t) (A->ia[j]); i < (size_t) (A->ia[j+1]); ++i) {
+            v[j] += A->sa[i]*u[A->ja[i]];
+        }
+    }
+}
+
 struct ifs_tpfa_impl {
     double *fgrav;              /* Accumulated grav contrib/face */
 
@@ -287,7 +297,7 @@ assemble_shut_well(int nc, int w,
     size_t jw;
     double trans;
 
-    /* The equation added for a shut well w is 
+    /* The equation added for a shut well w is
      *   \sum_{c~w} T_{w,c} * mt_c p_w = 0.0
      * where c~w indicates the cell perforated by w, T are production
      * indices, mt is total mobility and p_w is the bottom-hole
@@ -371,7 +381,7 @@ assemble_well_contrib(int                   nc ,
                 /* We cannot handle this case, since we do
                  * not have access to formation volume factors,
                  * needed to convert between reservoir and
-                 * surface rates 
+                 * surface rates
                  */
                 fprintf(stderr, "ifs_tpfa cannot handle SURFACE_RATE well controls.\n");
                 *ok = 0;
@@ -725,6 +735,49 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid      *G        ,
     return ok;
 }
 
+/* ---------------------------------------------------------------------- */
+int
+ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
+                                     const struct ifs_tpfa_forces *F        ,
+                                     const double                 *trans    ,
+                                     const double                 *gpress   ,
+                                     const double                 *porevol  ,
+                                     const double                 *rock_comp,
+                                     const double                  dt       ,
+                                     const double                 *prev_pressure,
+                                     const double                 *initial_porevolume,
+                                     struct ifs_tpfa_data         *h        )
+/* ---------------------------------------------------------------------- */
+{
+    int c;
+    size_t j;
+    int system_singular, ok;
+    double *v;
+
+    assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
+
+    v = malloc(h->A->m * sizeof *v);
+
+    mult_csr_matrix(h->A, prev_pressure, v);
+
+    /* We want to solve a Newton step for the residual
+     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible 
+     *
+     */
+
+    if (ok) {
+        for (c = 0; c < G->number_of_cells; ++c) {
+            j = csrmatrix_elm_index(c, c, h->A);
+            h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
+            h->b[c]     += -(porevol[c] - initial_porevolume[c])/dt - v[c];
+        }
+    }
+
+    free(v);
+
+    return ok;
+}
+
 
 /* ---------------------------------------------------------------------- */
 void
diff --git a/opm/core/pressure/tpfa/ifs_tpfa.h b/opm/core/pressure/tpfa/ifs_tpfa.h
index 8f8aed757..c06c2f660 100644
--- a/opm/core/pressure/tpfa/ifs_tpfa.h
+++ b/opm/core/pressure/tpfa/ifs_tpfa.h
@@ -61,6 +61,10 @@ struct ifs_tpfa_data *
 ifs_tpfa_construct(struct UnstructuredGrid *G,
                    struct Wells            *W);
 
+
+void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v);
+
+
 int
 ifs_tpfa_assemble(struct UnstructuredGrid      *G     ,
                   const struct ifs_tpfa_forces *F     ,
@@ -78,6 +82,19 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid      *G        ,
                            const double                  dt       ,
                            const double                 *pressure ,
                            struct ifs_tpfa_data         *h        );
+int
+ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
+				     const struct ifs_tpfa_forces *F        ,
+				     const double                 *trans    ,
+				     const double                 *gpress   ,
+				     const double                 *porevol  ,
+				     const double                 *rock_comp,
+				     const double                  dt       ,
+				     const double                 *prev_pressure ,
+				     const double                 *initial_porevolume,
+				     struct ifs_tpfa_data         *h        );
+
+
 void
 ifs_tpfa_press_flux(struct UnstructuredGrid      *G    ,
                     const struct ifs_tpfa_forces *F    ,