diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index 78fa851e67..80e13faa8d 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,15 +30,18 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
-
     #include "createTime.H"
     #include "createMesh.H"
+
+    pisoControl piso(mesh);
+
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -50,9 +53,10 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
+        // Momentum predictor
+
         fvVectorMatrix UEqn
         (
             fvm::ddt(U)
@@ -60,11 +64,13 @@ int main(int argc, char *argv[])
           - fvm::laplacian(nu, U)
         );
 
-        solve(UEqn == -fvc::grad(p));
+        if (piso.momentumPredictor())
+        {
+            solve(UEqn == -fvc::grad(p));
+        }
 
         // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        while (piso.correct())
         {
             volScalarField rAU(1.0/UEqn.A());
 
@@ -79,17 +85,21 @@ int main(int argc, char *argv[])
 
             adjustPhi(phiHbyA, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
+                // Pressure corrector
+
                 fvScalarMatrix pEqn
                 (
                     fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi = phiHbyA - pEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
index 30a73b70a6..9d3b567d3b 100644
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,15 +31,18 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
-
     #include "createTime.H"
     #include "createMeshNoClear.H"
+
+    pisoControl piso(mesh);
+
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -51,11 +54,12 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         fluid.correct();
 
+        // Momentum predictor
+
         fvVectorMatrix UEqn
         (
             fvm::ddt(U)
@@ -64,11 +68,13 @@ int main(int argc, char *argv[])
           - (fvc::grad(U) & fvc::grad(fluid.nu()))
         );
 
-        solve(UEqn == -fvc::grad(p));
+        if (piso.momentumPredictor())
+        {
+            solve(UEqn == -fvc::grad(p));
+        }
 
         // --- PISO loop
-
-        for (int corr=0; corr<nCorr; corr++)
+        while (piso.correct())
         {
             volScalarField rAU(1.0/UEqn.A());
 
@@ -83,17 +89,21 @@ int main(int argc, char *argv[])
 
             adjustPhi(phiHbyA, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
+                // Pressure corrector
+
                 fvScalarMatrix pEqn
                 (
                     fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi = phiHbyA - pEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/pisoFoam/pisoFoam.C b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
index e068306470..7bf7abb3ac 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -41,7 +41,6 @@ Description
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
-
     #include "createTime.H"
     #include "createMesh.H"
 
@@ -63,7 +62,6 @@ int main(int argc, char *argv[])
         // Pressure-velocity PISO corrector
         {
             // Momentum predictor
-
             fvVectorMatrix UEqn
             (
                 fvm::ddt(U)
diff --git a/src/finiteVolume/cfdTools/general/include/readPISOControls.H b/src/finiteVolume/cfdTools/general/include/readPISOControls.H
deleted file mode 100644
index 961acfc67e..0000000000
--- a/src/finiteVolume/cfdTools/general/include/readPISOControls.H
+++ /dev/null
@@ -1,17 +0,0 @@
-    const dictionary& pisoDict = mesh.solutionDict().subDict("PISO");
-
-    const int nOuterCorr =
-        pisoDict.lookupOrDefault<int>("nOuterCorrectors", 1);
-
-    const int nCorr =
-        pisoDict.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        pisoDict.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        pisoDict.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        pisoDict.lookupOrDefault("transonic", false);
-