diff --git a/applications/solvers/incompressible/icoDyMFoam/correctPhi.H b/applications/solvers/incompressible/icoDyMFoam/correctPhi.H
index 6c78a0960..d29b9e726 100644
--- a/applications/solvers/incompressible/icoDyMFoam/correctPhi.H
+++ b/applications/solvers/incompressible/icoDyMFoam/correctPhi.H
@@ -36,7 +36,7 @@
 
     mesh.setFluxRequired(pcorr.name());
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -46,7 +46,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/incompressible/icoDyMFoam/createFields.H b/applications/solvers/incompressible/icoDyMFoam/createFields.H
index 9c918b565..d4d011004 100644
--- a/applications/solvers/incompressible/icoDyMFoam/createFields.H
+++ b/applications/solvers/incompressible/icoDyMFoam/createFields.H
@@ -71,16 +71,9 @@
         zeroGradientFvPatchScalarField::typeName
     );
 
+    bool momentumPredictor = piso.momentumPredictor();
 
     {
-        const dictionary& piso = mesh.solutionDict().subDict("PISO");
-
-        bool momentumPredictor = true;
-        if (piso.found("momentumPredictor"))
-        {
-            momentumPredictor = Switch(piso.lookup("momentumPredictor"));
-        }
-
         bool momentumPredictorSave = momentumPredictor;
         momentumPredictor = false;
 #       include "UEqn.H"
diff --git a/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C b/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
index 84259e4eb..2bd63eb2a 100644
--- a/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
+++ b/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
@@ -35,6 +35,7 @@ Author
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,6 +44,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "initContinuityErrs.H"
 #   include "initTotalVolume.H"
 #   include "createFields.H"
@@ -90,7 +94,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             rAU = 1.0/UEqn.A();
 
@@ -100,7 +104,8 @@ int main(int argc, char *argv[])
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -109,7 +114,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
+                if (piso.finalInnerIter())
                 {
                     pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
                 }
@@ -118,7 +123,7 @@ int main(int argc, char *argv[])
                     pEqn.solve(mesh.solutionDict().solver(p.name()));
                 }
 
-                if (nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/icoDyMFoam/readControls.H b/applications/solvers/incompressible/icoDyMFoam/readControls.H
index 9336616a1..a6b4cce58 100644
--- a/applications/solvers/incompressible/icoDyMFoam/readControls.H
+++ b/applications/solvers/incompressible/icoDyMFoam/readControls.H
@@ -1,14 +1,11 @@
 #   include "readTimeControls.H"
-#   include "readPISOControls.H"
 
-    bool correctPhi = false;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
+    bool correctPhi
+    (
+        piso.dict().lookupOrDefault("correctPhi", false)
+    );
 
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+    bool checkMeshCourantNo
+    (
+        piso.dict().lookupOrDefault("checkMeshCourantNo", false)
+    );
diff --git a/applications/solvers/incompressible/icoDyMSimpleFoam/correctPhi.H b/applications/solvers/incompressible/icoDyMSimpleFoam/correctPhi.H
index 3f9a190d4..eb9730b95 100644
--- a/applications/solvers/incompressible/icoDyMSimpleFoam/correctPhi.H
+++ b/applications/solvers/incompressible/icoDyMSimpleFoam/correctPhi.H
@@ -26,7 +26,7 @@
 
     mesh.setFluxRequired(pcorr.name());
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -36,7 +36,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/incompressible/icoDyMSimpleFoam/createFields.H b/applications/solvers/incompressible/icoDyMSimpleFoam/createFields.H
index 9c918b565..d4d011004 100644
--- a/applications/solvers/incompressible/icoDyMSimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/icoDyMSimpleFoam/createFields.H
@@ -71,16 +71,9 @@
         zeroGradientFvPatchScalarField::typeName
     );
 
+    bool momentumPredictor = piso.momentumPredictor();
 
     {
-        const dictionary& piso = mesh.solutionDict().subDict("PISO");
-
-        bool momentumPredictor = true;
-        if (piso.found("momentumPredictor"))
-        {
-            momentumPredictor = Switch(piso.lookup("momentumPredictor"));
-        }
-
         bool momentumPredictorSave = momentumPredictor;
         momentumPredictor = false;
 #       include "UEqn.H"
diff --git a/applications/solvers/incompressible/icoDyMSimpleFoam/icoDyMSimpleFoam.C b/applications/solvers/incompressible/icoDyMSimpleFoam/icoDyMSimpleFoam.C
index b29c58ad9..bd86b5208 100644
--- a/applications/solvers/incompressible/icoDyMSimpleFoam/icoDyMSimpleFoam.C
+++ b/applications/solvers/incompressible/icoDyMSimpleFoam/icoDyMSimpleFoam.C
@@ -36,6 +36,7 @@ Author
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "initContinuityErrs.H"
 #   include "initTotalVolume.H"
 #   include "createFields.H"
@@ -89,7 +93,7 @@ int main(int argc, char *argv[])
 
         // --- SIMPLE loop
 
-        for (int ocorr = 0; ocorr < nOuterCorr; ocorr++)
+        while (piso.correct())
         {
 #           include "CourantNo.H"
 
@@ -105,7 +109,8 @@ int main(int argc, char *argv[])
 
             p.storePrevIter();
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -114,11 +119,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                if
-                (
-                    ocorr == nOuterCorr - 1
-                 && nonOrth == nNonOrthCorr
-                )
+                if(piso.finalInnerIter())
                 {
                     pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
                 }
@@ -127,7 +128,7 @@ int main(int argc, char *argv[])
                     pEqn.solve(mesh.solutionDict().solver(p.name()));
                 }
 
-                if (nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/icoDyMSimpleFoam/readControls.H b/applications/solvers/incompressible/icoDyMSimpleFoam/readControls.H
index 9336616a1..a6b4cce58 100644
--- a/applications/solvers/incompressible/icoDyMSimpleFoam/readControls.H
+++ b/applications/solvers/incompressible/icoDyMSimpleFoam/readControls.H
@@ -1,14 +1,11 @@
 #   include "readTimeControls.H"
-#   include "readPISOControls.H"
 
-    bool correctPhi = false;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
+    bool correctPhi
+    (
+        piso.dict().lookupOrDefault("correctPhi", false)
+    );
 
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+    bool checkMeshCourantNo
+    (
+        piso.dict().lookupOrDefault("checkMeshCourantNo", false)
+    );
diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index 164f5f6e9..c7efa80cf 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,6 +40,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -50,7 +54,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
 #       include "CourantNo.H"
 
         fvVectorMatrix UEqn
@@ -63,8 +66,7 @@ int main(int argc, char *argv[])
         solve(UEqn == -fvc::grad(p));
 
         // --- PISO loop
-
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             volScalarField rUA = 1.0/UEqn.A();
 
@@ -74,7 +76,8 @@ int main(int argc, char *argv[])
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -82,9 +85,17 @@ int main(int argc, char *argv[])
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
+                    
+                if (piso.finalInnerIter())
+                {
+                    pEqn.solve(mesh.solutionDict().solver("pFinal"));
+                }
+                else
+                {
+                    pEqn.solve();
+                }
 
-                if (nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }