diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp
index 687c2d2b7..38190071d 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp
@@ -1006,10 +1006,16 @@ namespace {
         rq_[ actph ].mob = tr_mult * kr[ phase ] / mu;
 
         const ADB rho   = fluidDensity(phase, state.pressure, state.rs, cells_);
-        const ADB gflux = grav_ * rho;
 
         ADB& head = rq_[ actph ].head;
-        head      = transi*(ops_.ngrad * state.pressure) + gflux;
+
+        // compute gravity potensial using the face average as in eclipse and MRST
+        const ADB rhoavg = ops_.caver * rho;
+
+        const ADB dp = ops_.ngrad * state.pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
+
+        head = transi*dp;
+        //head      = transi*(ops_.ngrad * state.pressure) + gflux;
 
         UpwindSelector<double> upwind(grid_, ops_, head.value());
 
diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp
index fbf646295..112cbc1de 100644
--- a/opm/autodiff/GeoProps.hpp
+++ b/opm/autodiff/GeoProps.hpp
@@ -46,6 +46,7 @@ namespace Opm
             : pvol_ (grid.number_of_cells)
             , trans_(grid.number_of_faces)
             , gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
+            , z_(grid.number_of_cells)
         {
             // Pore volume
             const typename Vector::Index nc = grid.number_of_cells;
@@ -59,6 +60,12 @@ namespace Opm
             tpfa_htrans_compute(ug, props.permeability(), htrans.data());
             tpfa_trans_compute (ug, htrans.data()     , trans_.data());
 
+            // Compute z coordinates
+            for (int c = 0; c<nc; ++c){
+                z_[c] = grid.cell_centroids[c*3 + 2];
+            }
+
+
             // Gravity potential
             std::fill(gravity_, gravity_ + 3, 0.0);
             if (grav != 0) {
@@ -82,15 +89,17 @@ namespace Opm
             }
         }
 
-        const Vector& poreVolume()       const { return pvol_ ; }
-        const Vector& transmissibility() const { return trans_; }
-        const Vector& gravityPotential() const { return gpot_ ; }
-        const double* gravity()          const { return gravity_; }
+        const Vector& poreVolume()       const { return pvol_   ;}
+        const Vector& transmissibility() const { return trans_  ;}
+        const Vector& gravityPotential() const { return gpot_   ;}
+        const Vector& z()                const { return z_      ;}
+        const double* gravity()          const { return gravity_;}
 
     private:
         Vector pvol_ ;
         Vector trans_;
         Vector gpot_ ;
+        Vector z_;
         double gravity_[3]; // Size 3 even if grid is 2-dim.
     };
 }
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
index 5b1b6718f..32ce8457c 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
@@ -199,7 +199,35 @@ namespace Opm
             std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
         }
     }
+    static void outputWellStateMatlab(const Opm::WellState& well_state,
+                                  const int step,
+                                  const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
 
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
 
 #if 0
     static void outputWaterCut(const Opm::Watercut& watercut,
@@ -341,6 +369,8 @@ namespace Opm
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+
             }
 
             SimulatorReport sreport;
@@ -455,6 +485,7 @@ namespace Opm
                 outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
             outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
 #if 0
             outputWaterCut(watercut, output_dir_);
             if (wells_) {