diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index 164f5f6e9..5073f4e61 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -26,6 +26,10 @@ Application
 
 Description
     Transient solver for incompressible, laminar flow of Newtonian fluids.
+    Consistent formulation without time-step and relaxation dependence by Jasak
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
 
 \*---------------------------------------------------------------------------*/
 
@@ -53,32 +57,36 @@ int main(int argc, char *argv[])
 #       include "readPISOControls.H"
 #       include "CourantNo.H"
 
-        fvVectorMatrix UEqn
+        // Convection-diffusion matrix
+        fvVectorMatrix HUEqn
         (
-            fvm::ddt(U)
-          + fvm::div(phi, U)
+            fvm::div(phi, U)
           - fvm::laplacian(nu, U)
         );
 
-        solve(UEqn == -fvc::grad(p));
+        // Time derivative matrix
+        fvVectorMatrix ddtUEqn(fvm::ddt(U));
+
+        solve(ddtUEqn + HUEqn == -fvc::grad(p));
+
+        // Prepare clean Ap without time derivative contribution
+        // HJ, 26/Oct/2015
+        volScalarField aU = HUEqn.A();
 
         // --- PISO loop
 
         for (int corr = 0; corr < nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
-
-            U = rUA*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
+            U = HUEqn.H()/aU;
+            phi = (fvc::interpolate(U) & mesh.Sf());
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
+                    fvm::laplacian(1/aU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -92,7 +100,12 @@ int main(int argc, char *argv[])
 
 #           include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            // 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()
+                );
             U.correctBoundaryConditions();
         }
 
diff --git a/tutorials/incompressible/icoFoam/elbow/system/fvSchemes b/tutorials/incompressible/icoFoam/elbow/system/fvSchemes
index 0cd903c76..7c96ea2c4 100644
--- a/tutorials/incompressible/icoFoam/elbow/system/fvSchemes
+++ b/tutorials/incompressible/icoFoam/elbow/system/fvSchemes
@@ -28,7 +28,7 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,U)      Gauss linearUpwind leastSquares;
 }
 
 laplacianSchemes
diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U
index c9b0cd956..e5afe4b5d 100644
--- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U
+++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U
@@ -54,12 +54,6 @@ boundaryField
         type            mixingPlane;
     }
 
-    top
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U
index c9b0cd956..e5afe4b5d 100644
--- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U
+++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U
@@ -54,12 +54,6 @@ boundaryField
         type            mixingPlane;
     }
 
-    top
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U
index c9b0cd956..e5afe4b5d 100644
--- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U
+++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U
@@ -54,12 +54,6 @@ boundaryField
         type            mixingPlane;
     }
 
-    top
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
     frontAndBack
     {
         type            empty;