diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp
index daa68b65..da4270e1 100644
--- a/examples/spu_2p.cpp
+++ b/examples/spu_2p.cpp
@@ -324,7 +324,11 @@ main(int argc, char** argv)
         // Gravity.
         gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
         // Init state variables (saturation and pressure).
-        initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+	if (param.has("init_saturation")) {
+	    initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state);
+	} else {
+	    initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+	}
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
@@ -518,10 +522,11 @@ main(int argc, char** argv)
     std::vector<double> well_perfrates;
     std::vector<double> fractional_flows;
     std::vector<double> well_resflows_phase;
+    int num_wells = 0;
     if (wells->c_wells()) {
-        const int nw = wells->c_wells()->number_of_wells;
-        well_bhp.resize(nw, 0.0);
-        well_perfrates.resize(wells->c_wells()->well_connpos[nw], 0.0);
+        num_wells = wells->c_wells()->number_of_wells;
+        well_bhp.resize(num_wells, 0.0);
+        well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0);
         well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0);
         wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates);
     }
@@ -557,15 +562,17 @@ main(int argc, char** argv)
                 std::vector<double> initial_pressure = state.pressure();
 		std::vector<double> initial_porevolume(num_cells);
 		computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
-		std::vector<double> pressure_increment(num_cells);
-                std::vector<double> prev_pressure;
+		std::vector<double> pressure_increment(num_cells + num_wells);
+                std::vector<double> prev_pressure(num_cells + num_wells);
 		for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
 		
 		    for (int cell = 0; cell < num_cells; ++cell) {
 			rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
 		    }
 		    computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
-		    prev_pressure = state.pressure();
+		    std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
+		    std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells);
+		    // prev_pressure = state.pressure();
 		
 		    // compute pressure increment
 		    psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, 
@@ -577,6 +584,10 @@ main(int argc, char** argv)
 			state.pressure()[cell] += pressure_increment[cell];
 			max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
 		    }
+		    for (int well = 0; well < num_wells; ++well) {
+			well_bhp[well] += pressure_increment[num_cells + well];
+			max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
+		    }
 		
 		    std::cout << "Pressure iter " << iter << "   max change = " << max_change << std::endl;
 		    if (max_change < nl_pressure_tolerance) {
diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp
index c7222b94..62618651 100644
--- a/opm/core/pressure/IncompTpfa.cpp
+++ b/opm/core/pressure/IncompTpfa.cpp
@@ -336,7 +336,10 @@ namespace Opm
             soln.well_press = &well_bhp[0];
         }
 
-	memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
+	// memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
+	ASSERT(int(pressure.size()) == grid_.number_of_cells);
+	std::copy(pressure.begin(), pressure.end(), h_->x);
+	std::copy(well_bhp.begin(), well_bhp.end(), h_->x + grid_.number_of_cells);
 
         ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
     }
diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c
index 3bc33be2..eb84dfdd 100644
--- a/opm/core/pressure/tpfa/ifs_tpfa.c
+++ b/opm/core/pressure/tpfa/ifs_tpfa.c
@@ -760,7 +760,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
                                      struct ifs_tpfa_data         *h        )
 /* ---------------------------------------------------------------------- */
 {
-    int     c, system_singular, ok;
+    int     c, w, system_singular, ok;
     size_t  j;
     double *v;
 
@@ -780,6 +780,11 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
             h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
             h->b[c]     += -(porevol[c] - initial_porevolume[c])/dt - v[c];
         }
+	if (F->W != NULL) {
+	    for (w = 0; w < F->W->number_of_wells; ++w) {
+		h->b[c] += -v[c];
+	    }
+	}
     }
 
     return ok;