diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp
index f57b4d46..ea126dba 100644
--- a/examples/wells_example.cpp
+++ b/examples/wells_example.cpp
@@ -57,8 +57,7 @@ int main(int argc, char** argv) {
     computeTotalMobilityOmega(incomp_properties, all_cells, state.saturation(), totmob, omega);
     
     std::vector<double> wdp;
-    std::vector<double> densities(incomp_properties.density(), incomp_properties.density() + incomp_properties.numPhases());
-    computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), densities, wdp);
+    computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), incomp_properties.density(), wdp, true);
     
     std::vector<double> src;
     Opm::FlowBCManager bcs;
diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp
index bd38b8cc..1eb33938 100644
--- a/opm/core/utility/miscUtilities.cpp
+++ b/opm/core/utility/miscUtilities.cpp
@@ -405,10 +405,12 @@ namespace Opm
     }
 
     void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
-            const std::vector<double>& densities, std::vector<double>& wdp, bool per_grid_cell)
+                    const double* densities, std::vector<double>& wdp, bool per_grid_cell)
     {
-        const size_t np = densities.size();
         const int nw = wells.number_of_wells;
+        const size_t np = per_grid_cell ?
+            saturations.size()/grid.number_of_cells
+            : saturations.size()/wells.well_connpos[nw];
         // Simple for now:
         for (int i = 0; i < nw; i++) {
             double depth_ref = wells.depth_ref[i];
diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp
index 40e0c6be..dece112e 100644
--- a/opm/core/utility/miscUtilities.hpp
+++ b/opm/core/utility/miscUtilities.hpp
@@ -188,7 +188,7 @@ namespace Opm
     /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per 
     ///                          well cell.
     void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
-                    const std::vector<double>& densities, std::vector<double>& wdp, bool per_grid_cell = true);
+                    const double* densities, std::vector<double>& wdp, bool per_grid_cell);
     
     /// Computes (sums) the flow rate for each well. 
     /// \param[in] wells                The wells for which the flow rate should be computed.