diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp
new file mode 100644
index 000000000..17ad7b85b
--- /dev/null
+++ b/examples/compute_tof_from_files.cpp
@@ -0,0 +1,177 @@
+/*
+  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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/GridManager.hpp>
+#include <opm/core/newwells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/initState.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/fluid/IncompPropertiesBasic.hpp>
+#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/pressure/IncompTpfa.hpp>
+#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
+#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
+
+#include <boost/scoped_ptr.hpp>
+#include <boost/filesystem.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <iterator>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Warning: unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "-------------------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+{
+    using namespace Opm;
+
+    parameter::ParameterGroup param(argc, argv, false);
+
+    // Read grid.
+    GridManager grid_manager(param.get<std::string>("grid_filename"));
+    const UnstructuredGrid& grid = *grid_manager.c_grid();
+
+    // Read porosity, compute pore volume.
+    std::vector<double> porevol;
+    {
+        std::ifstream poro_stream(param.get<std::string>("poro_filename").c_str());
+        std::istream_iterator<double> beg(poro_stream);
+        std::istream_iterator<double> end;
+        porevol.assign(beg, end); // Now contains poro.
+        if (int(porevol.size()) != grid.number_of_cells) {
+            THROW("Size of porosity field differs from number of cells.");
+        }
+        for (int i = 0; i < grid.number_of_cells; ++i) {
+            porevol[i] *= grid.cell_volumes[i];
+        }
+    }
+
+    // Read flux.
+    std::vector<double> flux;
+    {
+        std::ifstream flux_stream(param.get<std::string>("flux_filename").c_str());
+        std::istream_iterator<double> beg(flux_stream);
+        std::istream_iterator<double> end;
+        flux.assign(beg, end);
+        if (int(flux.size()) != grid.number_of_faces) {
+            THROW("Size of flux field differs from number of faces.");
+        }
+    }
+
+    // Read source terms.
+    std::vector<double> src;
+    {
+        std::ifstream src_stream(param.get<std::string>("src_filename").c_str());
+        std::istream_iterator<double> beg(src_stream);
+        std::istream_iterator<double> end;
+        src.assign(beg, end);
+        if (int(src.size()) != grid.number_of_cells) {
+            THROW("Size of source term field differs from number of cells.");
+        }
+    }
+
+    // Choice of tof solver.
+    bool use_dg = param.getDefault("use_dg", false);
+    int dg_degree = -1;
+    bool use_cvi = false;
+    bool use_multidim_upwind = false;
+    if (use_dg) {
+        dg_degree = param.getDefault("dg_degree", 0);
+        use_cvi = param.getDefault("use_cvi", false);
+    } else {
+        use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
+    }
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::ofstream epoch_os;
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            THROW("Creating directories failed: " << fpath);
+        }
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    // Issue a warning if any parameters were unused.
+    warnIfUnusedParams(param);
+
+    // Solve time-of-flight.
+    Opm::time::StopWatch transport_timer;
+    transport_timer.start();
+    std::vector<double> tof;
+    if (use_dg) {
+        Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi);
+        tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
+    } else {
+        Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
+        tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
+    }
+    transport_timer.stop();
+    double tt = transport_timer.secsSinceStart();
+    std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+
+    // Output.
+    if (output) {
+        std::string tof_filename = output_dir + "/tof.txt";
+        std::ofstream tof_stream(tof_filename.c_str());
+        tof_stream.precision(16);
+        std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
+    }
+}
diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp
index 8d93f2480..91b20ef01 100644
--- a/opm/core/transport/reorder/TransportModelTracerTof.cpp
+++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp
@@ -30,8 +30,9 @@ namespace Opm
 
     /// Construct solver.
     /// \param[in] grid      A 2d or 3d grid.
-    TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid)
-        : grid_(grid)
+    TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid,
+                                                     const bool use_multidim_upwind)
+        : grid_(grid), use_multidim_upwind_(use_multidim_upwind)
     {
     }
 
@@ -63,6 +64,10 @@ namespace Opm
         tof.resize(grid_.number_of_cells);
         std::fill(tof.begin(), tof.end(), 0.0);
         tof_ = &tof[0];
+        if (use_multidim_upwind_) {
+            face_tof_.resize(grid_.number_of_faces);
+            std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
+        }
         reorderAndTransport(grid_, darcyflux);
     }
 
@@ -71,6 +76,10 @@ namespace Opm
 
     void TransportModelTracerTof::solveSingleCell(const int cell)
     {
+        if (use_multidim_upwind_) {
+            solveSingleCellMultidimUpwind(cell);
+            return;
+        }
         // Compute flux terms.
         // Sources have zero tof, and therefore do not contribute
         // to upwind_term. Sinks on the other hand, must be added
@@ -94,7 +103,7 @@ namespace Opm
             // Add flux to upwind_term or downwind_flux, if interior.
             if (other != -1) {
                 if (flux < 0.0) {
-                    upwind_term  += flux*tof_[other];
+                    upwind_term += flux*tof_[other];
                 } else {
                     downwind_flux += flux;
                 }
@@ -108,6 +117,60 @@ namespace Opm
 
 
 
+    void TransportModelTracerTof::solveSingleCellMultidimUpwind(const int cell)
+    {
+        // Compute flux terms.
+        // Sources have zero tof, and therefore do not contribute
+        // to upwind_term. Sinks on the other hand, must be added
+        // to the downwind terms (note sign change resulting from
+        // different sign conventions: pos. source is injection,
+        // pos. flux is outflow).
+        double upwind_term = 0.0;
+        double downwind_term_cell_factor = std::max(-source_[cell], 0.0);
+        double downwind_term_face = 0.0;
+        for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
+            int f = grid_.cell_faces[i];
+            double flux;
+            int other;
+            // Compute cell flux
+            if (cell == grid_.face_cells[2*f]) {
+                flux  = darcyflux_[f];
+                other = grid_.face_cells[2*f+1];
+            } else {
+                flux  =-darcyflux_[f];
+                other = grid_.face_cells[2*f];
+            }
+            // Add flux to upwind_term or downwind_flux, if interior.
+            if (other != -1) {
+                if (flux < 0.0) {
+                    upwind_term += flux*face_tof_[f];
+                } else {
+                    double fterm, cterm_factor;
+                    multidimUpwindTerms(f, cell, fterm, cterm_factor);
+                    downwind_term_face += fterm*flux;
+                    downwind_term_cell_factor += cterm_factor*flux;
+                }
+            }
+        }
+
+        // Compute tof for cell.
+        tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor;        // }
+
+        // Compute tof for downwind faces.
+        for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
+            int f = grid_.cell_faces[i];
+            const double outflux_f = (grid_.face_cells[2*f] == cell) ? darcyflux_[f] : -darcyflux_[f];
+            if (outflux_f > 0.0) {
+                double fterm, cterm_factor;
+                multidimUpwindTerms(f, cell, fterm, cterm_factor);
+                face_tof_[f] = fterm + cterm_factor*tof_[cell];
+            }
+        }
+    }
+
+
+
+
     void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
     {
         std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
@@ -119,4 +182,77 @@ namespace Opm
 
 
 
+    // Assumes that face_tof_[f] is known for all upstream faces f of upwind_cell.
+    // Assumes that darcyflux_[face] is != 0.0.
+    // This function returns factors to compute the tof for 'face':
+    //   tof(face) = face_term + cell_term_factor*tof(upwind_cell).
+    // It is not computed here, since these factors are needed to
+    // compute the tof(upwind_cell) itself.
+    void TransportModelTracerTof::multidimUpwindTerms(const int face,
+                                                      const int upwind_cell,
+                                                      double& face_term,
+                                                      double& cell_term_factor) const
+    {
+        // Implements multidim upwind according to
+        // "Multidimensional upstream weighting for multiphase transport on general grids"
+        // by Keilegavlen, Kozdon, Mallison.
+        // However, that article does not give a 3d extension other than noting that using
+        // multidimensional upwinding in the XY-plane and not in the Z-direction may be
+        // a good idea. We have here attempted some generalization, by looking at all
+        // face-neighbours across edges as upwind candidates, and giving them all uniform weight.
+        // This will over-weight the immediate upstream cell value in an extruded 2d grid with
+        // one layer (top and bottom no-flow faces will enter the computation) compared to the
+        // original 2d case. Improvements are welcome.
+
+        // Identify the adjacent faces of the upwind cell.
+        const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face];
+        const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[face + 1];
+        ASSERT(face_nodes_end - face_nodes_beg == 2 || grid_.dimensions != 2);
+        adj_faces_.clear();
+        for (int hf = grid_.cell_facepos[upwind_cell]; hf < grid_.cell_facepos[upwind_cell + 1]; ++hf) {
+            const int f = grid_.cell_faces[hf];
+            if (f != face) {
+                const int* f_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
+                const int* f_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
+                // Find out how many vertices they have in common.
+                // Using simple linear searches since sets are small.
+                int num_common = 0;
+                for (const int* f_iter = f_nodes_beg; f_iter < f_nodes_end; ++f_iter) {
+                    num_common += std::count(face_nodes_beg, face_nodes_end, *f_iter);
+                }
+                if (num_common == grid_.dimensions - 1) {
+                    // Neighbours over an edge (3d) or vertex (2d).
+                    adj_faces_.push_back(f);
+                } else {
+                    ASSERT(num_common == 0);
+                }
+            }
+        }
+
+        // Indentify adjacent faces with inflows, compute omega_star, omega,
+        // add up contributions.
+        const int num_adj = adj_faces_.size();
+        ASSERT(num_adj == face_nodes_end - face_nodes_beg);
+        const double flux_face = std::fabs(darcyflux_[face]);
+        face_term = 0.0;
+        cell_term_factor = 0.0;
+        for (int ii = 0; ii < num_adj; ++ii) {
+            const int f = adj_faces_[ii];
+            const double influx_f = (grid_.face_cells[2*f] == upwind_cell) ? -darcyflux_[f] : darcyflux_[f];
+            const double omega_star = influx_f/flux_face;
+            // SPU
+            // const double omega = 0.0;
+            // TMU
+            // const double omega = omega_star > 0.0 ? std::min(omega_star, 1.0) :  0.0;
+            // SMU
+            const double omega = omega_star > 0.0 ? omega_star/(1.0 + omega_star) : 0.0;
+            face_term += omega*face_tof_[f];
+            cell_term_factor += (1.0 - omega);
+        }
+        face_term /= double(num_adj);
+        cell_term_factor /= double(num_adj);
+    }
+
+
+
 } // namespace Opm
diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp
index 270ac328a..d9540abd2 100644
--- a/opm/core/transport/reorder/TransportModelTracerTof.hpp
+++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp
@@ -43,7 +43,9 @@ namespace Opm
     public:
         /// Construct solver.
         /// \param[in] grid      A 2d or 3d grid.
-        TransportModelTracerTof(const UnstructuredGrid& grid);
+        /// \param[in] use_multidim_upwind  If true, use multidimensional tof upwinding.
+        TransportModelTracerTof(const UnstructuredGrid& grid,
+                                const bool use_multidim_upwind = false);
 
         /// Solve for time-of-flight.
         /// \param[in]  darcyflux         Array of signed face fluxes.
@@ -59,14 +61,21 @@ namespace Opm
 
     private:
         virtual void solveSingleCell(const int cell);
+        void solveSingleCellMultidimUpwind(const int cell);
         virtual void solveMultiCell(const int num_cells, const int* cells);
 
+        void multidimUpwindTerms(const int face, const int upwind_cell,
+                                 double& face_term, double& cell_term_factor) const;
+
     private:
         const UnstructuredGrid& grid_;
         const double* darcyflux_;   // one flux per grid face
         const double* porevolume_;  // one volume per cell
         const double* source_;      // one volumetric source term per cell
         double* tof_;
+        bool use_multidim_upwind_;
+        std::vector<double> face_tof_;       // For multidim upwind face tofs.
+        mutable std::vector<int> adj_faces_; // For multidim upwind logic.
     };
 
 } // namespace Opm