diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index c34eb3543..da069b8cd 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -26,7 +26,8 @@ Application
 
 Description
     Transient solver for incompressible, laminar flow of Newtonian fluids.
-    Consistent formulation without time-step and relaxation dependence by Jasak
+    Consistent formulation without time-step and relaxation dependence by
+    Jasak and Tukovic.
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved
@@ -66,24 +67,23 @@ int main(int argc, char *argv[])
           - fvm::laplacian(nu, U)
         );
 
-        // Time derivative matrix
-        fvVectorMatrix ddtUEqn(fvm::ddt(U));
-
         if (piso.momentumPredictor())
         {
-            solve(ddtUEqn + HUEqn == -fvc::grad(p));
+            solve(fvm::ddt(U) + HUEqn == -fvc::grad(p));
         }
 
-        // Prepare clean Ap without time derivative contribution
-        // HJ, 26/Oct/2015
-        volScalarField aU = HUEqn.A();
+        // Prepare clean 1/a_p without time derivative contribution
+        volScalarField rAU = 1.0/HUEqn.A();
 
         // --- PISO loop
 
         while (piso.correct())
         {
-            U = HUEqn.H()/aU;
-            phi = (fvc::interpolate(U) & mesh.Sf());
+            // Calculate U from convection-diffusion matrix
+            U = rAU*HUEqn.H();
+            
+            // Consistently calculate flux and face velocity
+            phi = piso.timeConsistentFlux(U, rAU);
 
             adjustPhi(phi, U, p);
 
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(1/aU, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -108,12 +108,8 @@ int main(int argc, char *argv[])
 
 #           include "continuityErrs.H"
 
-            // Note: cannot call H(U) here because the velocity is not complete
-            // HJ, 22/Jan/2016
-            U = 1.0/(aU + ddtUEqn.A())*
-                (
-                    U*aU - fvc::grad(p) + ddtUEqn.H()
-                );
+            // Consistently reconstruct velocity after pressure equation
+            U = piso.timeConsistentVelocity(U, rAU, p);
             U.correctBoundaryConditions();
         }