diff --git a/applications/solvers/basic/laplacianFoam/laplacianFoam.C b/applications/solvers/basic/laplacianFoam/laplacianFoam.C
index 830a266a40..5ed2849be9 100644
--- a/applications/solvers/basic/laplacianFoam/laplacianFoam.C
+++ b/applications/solvers/basic/laplacianFoam/laplacianFoam.C
@@ -52,7 +52,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+        while (simple.correctNonOrthogonal())
         {
             solve
             (
diff --git a/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C b/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
index 7dbbc97dac..8eade6a545 100644
--- a/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
+++ b/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+        while (simple.correctNonOrthogonal())
         {
             solve
             (
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index 40a48cfbdb..bf70e83c5d 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -114,12 +114,12 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "bEqn.H"
                 #include "ftEqn.H"
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
index 557c8aa1fb..0c5082df04 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@@ -167,13 +167,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "bEqn.H"
                 #include "ftEqn.H"
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 3c6a6341f9..e2a2a471e2 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -24,12 +24,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*invA, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -44,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -53,12 +50,9 @@ else
           - fvm::laplacian(rho*invA, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index 71b03064fa..069e5b8caa 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
@@ -105,8 +105,8 @@ int main(int argc, char *argv[])
                 hu == h;
             }
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index bef433d388..784e9ca6b7 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -24,12 +24,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -44,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -53,12 +50,9 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index 7a0f3acd09..6ef3a599c1 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -74,12 +74,12 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "hEqn.H"
                 #include "pEqn.H"
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 6fd1981be9..3c60e9c9ea 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
@@ -83,14 +83,14 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/dieselEngineFoam/pEqn.H b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
index 42228cc7b0..0493a1ff6c 100644
--- a/applications/solvers/combustion/dieselEngineFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
@@ -12,7 +12,7 @@ if (pimple.transonic())
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -23,12 +23,9 @@ if (pimple.transonic())
             Sevap
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -39,7 +36,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -50,12 +47,9 @@ else
             Sevap
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index 70808e8dbe..5bcb6b3273 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C
@@ -80,14 +80,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/dieselFoam/pEqn.H b/applications/solvers/combustion/dieselFoam/pEqn.H
index 5c4bd5ddf6..8ef13cde9d 100644
--- a/applications/solvers/combustion/dieselFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,12 +26,9 @@ if (pimple.transonic())
             Sevap
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -46,7 +43,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -57,14 +54,11 @@ else
             Sevap
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi == pEqn.flux();
         }
     }
 }
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 5b2256ba45..84472420c5 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
@@ -111,8 +111,8 @@ int main(int argc, char *argv[])
                 hu == h;
             }
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 1b75a93398..580be2da56 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -12,7 +12,7 @@ if (pimple.transonic())
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -21,12 +21,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -37,7 +34,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -46,12 +43,9 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index e95a4318c7..1d58470cbf 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -91,13 +91,13 @@ int main(int argc, char *argv[])
             #include "rhoEqn.H"
 
             // --- PIMPLE loop
-            for (pimple.start(); pimple.loop(); pimple++)
+            while (pimple.loop())
             {
                 #include "UEqn.H"
                 #include "YhsEqn.H"
 
-                // --- PISO loop
-                for (int corr=1; corr<=pimple.nCorr(); corr++)
+                // --- Pressure corrector loop
+                while (pimple.correct())
                 {
                     #include "pEqn.H"
                 }
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 5c7a4ae33e..65259391da 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -15,7 +15,7 @@ surfaceScalarField phiU
 
 phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
@@ -28,12 +28,9 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
       + surfaceFilm.Srho()
     );
 
-    p_rghEqn.solve
-    (
-        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-    );
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    if (pimple.finalNonOrthogonalIter())
     {
         phi += p_rghEqn.flux();
     }
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index bef433d388..784e9ca6b7 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -24,12 +24,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -44,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -53,12 +50,9 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 7f3347f780..5344eea3be 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -66,14 +66,14 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index 724f45e189..b22fc89079 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -30,7 +30,7 @@
           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
         );
 
-        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -38,12 +38,9 @@
               - fvm::laplacian(rho*rAU, p)
             );
 
-            pEqn.solve
-            (
-                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            );
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-            if (nonOrth == pimple.nNonOrthCorr())
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi += pEqn.flux();
             }
@@ -64,7 +61,7 @@
           + fvc::div(phi)
         );
 
-        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -72,12 +69,9 @@
               - fvm::laplacian(rho*rAU, p)
             );
 
-            pEqn.solve
-            (
-                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            );
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-            if (nonOrth == pimple.nNonOrthCorr())
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi += pEqn.flux();
             }
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index 73b252e579..aa8f1a5501 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -68,14 +68,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index eea0ea9bff..1361c6441f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 U = rAU*UEqn().H();
 
-if (pimple.nCorr() <= 1)
+if (pimple.nCorrPIMPLE() <= 1)
 {
     UEqn.clear();
 }
@@ -22,7 +22,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -31,12 +31,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -51,7 +48,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -61,12 +58,9 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 0336ceb74d..26d11a14aa 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -69,13 +69,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "hEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index 66ffb08275..de32f1f5c6 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -77,15 +77,15 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             turbulence->correct();
 
             #include "UEqn.H"
             #include "hEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index 5e54f970cd..cd10ed78e2 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -6,7 +6,7 @@ rho.relax();
 volScalarField rAU(1.0/UEqn().A());
 U = rAU*UEqn().H();
 
-if (pimple.nCorr() <= 1)
+if (pimple.nCorrPIMPLE() <= 1)
 {
     UEqn.clear();
 }
@@ -24,7 +24,7 @@ if (pimple.transonic())
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -33,12 +33,9 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -54,7 +51,7 @@ else
         );
     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -64,12 +61,9 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 8de8a412c3..068de89952 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -72,13 +72,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "hEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index df120f6dea..4d3f1ac76e 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -17,7 +17,7 @@ if (simple.transonic())
         fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -32,7 +32,7 @@ if (simple.transonic())
 
         pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -43,7 +43,7 @@ else
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -54,7 +54,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index 4fc258870b..b1f83304cd 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -20,7 +20,7 @@ if (simple.transonic())
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -37,7 +37,7 @@ if (simple.transonic())
 
         tpEqn().solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi == tpEqn().flux();
         }
@@ -50,7 +50,7 @@ else
 
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -67,7 +67,7 @@ else
 
         tpEqn().solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= tpEqn().flux();
         }
diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
index 9eb783feca..ccc7a1b21e 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
@@ -14,7 +14,7 @@ bool closedVolume = false;
 
 if (simple.transonic())
 {
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         surfaceScalarField phid
         (
@@ -44,17 +44,9 @@ if (simple.transonic())
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // Retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
+        pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi == phic + pEqn.flux();
         }
@@ -62,7 +54,7 @@ if (simple.transonic())
 }
 else
 {
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         phi = fvc::interpolate(rho*U) & mesh.Sf();
         closedVolume = adjustPhi(phi, U, p);
@@ -77,18 +69,9 @@ else
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // Retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
+        pEqn.solve();
 
-
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C b/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
index c0311807d4..37bea3c176 100644
--- a/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
+++ b/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
 
     runTime++;
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         solve(fvm::laplacian(murf, psi) + fvc::div(murf*Mrf));
     }
diff --git a/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H
index 17bf590f29..d2793d7192 100644
--- a/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H
@@ -14,7 +14,7 @@
     surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -24,7 +24,7 @@
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index c00c1fe8d8..461578e864 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -79,13 +79,13 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "TEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 8ac13dc93e..14e87ea844 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -10,7 +10,7 @@
     surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -19,12 +19,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 5d643837ff..5a9ec66707 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -11,7 +11,7 @@
     surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -22,7 +22,7 @@
 
         p_rghEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 4239771168..3e6b46a4bc 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -72,13 +72,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "hEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index b47d65b842..5c60d78b35 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -25,7 +25,7 @@
       + fvc::div(phi)
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -33,12 +33,9 @@
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi += p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 17bf590f29..d2793d7192 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -14,7 +14,7 @@
     surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -24,7 +24,7 @@
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index 1ba675c4bf..bedc12b3ba 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -120,7 +120,7 @@ int main(int argc, char *argv[])
             adjustPhi(phi, U, p);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+            while (simple.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -130,7 +130,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == simple.nNonOrthCorr())
+                if (simple.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }
@@ -181,7 +181,7 @@ int main(int argc, char *argv[])
             adjustPhi(phia, Ua, pa);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+            while (simple.correctNonOrthogonal())
             {
                 fvScalarMatrix paEqn
                 (
@@ -191,7 +191,7 @@ int main(int argc, char *argv[])
                 paEqn.setReference(paRefCell, paRefValue);
                 paEqn.solve();
 
-                if (nonOrth == simple.nNonOrthCorr())
+                if (simple.finalNonOrthogonalIter())
                 {
                     phia -= paEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
index 30725cd070..ba222b5034 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
@@ -65,12 +65,12 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UrelEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index e4f9190cc5..a8ade8c1ca 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -1,7 +1,7 @@
 volScalarField rAUrel(1.0/UrelEqn().A());
 Urel = rAUrel*UrelEqn().H();
 
-if (pimple.nCorr() <= 1)
+if (pimple.nCorrPIMPLE() <= 1)
 {
     UrelEqn.clear();
 }
@@ -12,7 +12,7 @@ phi = (fvc::interpolate(Urel) & mesh.Sf())
 adjustPhi(phi, Urel, p);
 
 // Non-orthogonal pressure corrector loop
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     // Pressure corrector
     fvScalarMatrix pEqn
@@ -22,12 +22,9 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve
-    (
-        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-    );
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    if (pimple.finalNonOrthogonalIter())
     {
         phi -= pEqn.flux();
     }
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index da9b7581bf..2d6b37db67 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
 U = rAU*UEqn().H();
 
-if (pimple.nCorr() <= 1)
+if (pimple.nCorrPIMPLE() <= 1)
 {
     UEqn.clear();
 }
@@ -11,7 +11,7 @@ phi = (fvc::interpolate(U) & mesh.Sf())
 adjustPhi(phi, U, p);
 
 // Non-orthogonal pressure corrector loop
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     // Pressure corrector
     fvScalarMatrix pEqn
@@ -21,12 +21,9 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve
-    (
-        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-    );
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    if (pimple.finalNonOrthogonalIter())
     {
         phi -= pEqn.flux();
     }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
index e112db4621..aed4884c4c 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
@@ -50,7 +50,7 @@
         pcorrTypes
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -60,7 +60,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 3e8a587768..5986a92adf 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
 U = rAU*UEqn().H();
 
-if (pimple.nCorr() <= 1)
+if (pimple.nCorrPIMPLE() <= 1)
 {
     UEqn.clear();
 }
@@ -19,7 +19,7 @@ if (p.needReference())
     fvc::makeAbsolute(phi, U);
 }
 
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
@@ -28,12 +28,9 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve
-    (
-        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-    );
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    if (pimple.finalNonOrthogonalIter())
     {
         phi -= pEqn.flux();
     }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
index 55ee07dbc7..b2b0c48c13 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
@@ -86,12 +86,12 @@ int main(int argc, char *argv[])
         }
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index e8b58bbea5..da5e2974c9 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -64,12 +64,12 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index 9fbe918a9f..b491552647 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
         #include "CourantNo.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
 
@@ -89,8 +89,8 @@ int main(int argc, char *argv[])
                 }
             }
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 volScalarField rAU(1.0/hUEqn.A());
                 surfaceScalarField ghrAUf(magg*fvc::interpolate(h*rAU));
@@ -110,7 +110,7 @@ int main(int argc, char *argv[])
                     + fvc::ddtPhiCorr(rAU, h, hU, phi)
                     - phih0;
 
-                for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+                while (pimple.correctNonOrthogonal())
                 {
                     fvScalarMatrix hEqn
                     (
@@ -119,15 +119,9 @@ int main(int argc, char *argv[])
                       - fvm::laplacian(ghrAUf, h)
                     );
 
-                    hEqn.solve
-                    (
-                        mesh.solver
-                        (
-                            h.select(pimple.finalInnerIter(corr, nonOrth))
-                        )
-                    );
+                    hEqn.solve(mesh.solver(h.select(pimple.finalInnerIter())));
 
-                    if (nonOrth == pimple.nNonOrthCorr())
+                    if (pimple.finalNonOrthogonalIter())
                     {
                         phi += hEqn.flux();
                     }
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
index 81f5c20480..2d2540d303 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
@@ -10,7 +10,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -20,7 +20,7 @@
         pEqn.setReference(pRefCell, pRefValue);
         pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index beb73bd341..54eab9cba7 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -9,7 +9,7 @@
     adjustPhi(phi, Urel, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -20,7 +20,7 @@
 
         pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 699cdcb3cf..966dedb616 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -9,7 +9,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -20,7 +20,7 @@
 
         pEqn.solve();
 
-        if (nonOrth == simple.nNonOrthCorr())
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
index 5908396d18..fbe81daee5 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
@@ -11,7 +11,7 @@ UEqn.clear();
 phi = fvc::interpolate(U) & mesh.Sf();
 adjustPhi(phi, U, p);
 
-for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+while (simple.correctNonOrthogonal())
 {
     tmp<fvScalarMatrix> tpEqn;
 
@@ -25,17 +25,10 @@ for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     }
 
     tpEqn().setReference(pRefCell, pRefValue);
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
-    {
-        tpEqn().solve();
-    }
-    else
-    {
-        tpEqn().solve();
-    }
 
-    if (nonOrth == simple.nNonOrthCorr())
+    tpEqn().solve();
+
+    if (simple.finalNonOrthogonalIter())
     {
         phi -= tpEqn().flux();
     }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index dc642e7c72..23c9e03229 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             turbulence->correct();
 
@@ -95,8 +95,8 @@ int main(int argc, char *argv[])
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index fcbfad3efb..d3aaed21ef 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -32,7 +32,7 @@
       + massSource.SuTot()
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -40,12 +40,9 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 7b5c3e73c7..897f43e3c3 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -88,14 +88,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index de0c04577a..73e56f07c5 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,12 +26,9 @@ if (pimple.transonic())
             coalParcels.Srho()
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -46,7 +43,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -57,12 +54,9 @@ else
             coalParcels.Srho()
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 63d746fe79..b6e69acead 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -32,7 +32,7 @@
       + massSource.SuTot()
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -40,12 +40,9 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index 283773f4f6..19d7e96d13 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
@@ -90,14 +90,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 5c7a4ae33e..65259391da 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -15,7 +15,7 @@ surfaceScalarField phiU
 
 phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
@@ -28,12 +28,9 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
       + surfaceFilm.Srho()
     );
 
-    p_rghEqn.solve
-    (
-        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-    );
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    if (pimple.finalNonOrthogonalIter())
     {
         phi += p_rghEqn.flux();
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 8dd3e3c9e6..06ef5ddac3 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -83,14 +83,14 @@ int main(int argc, char *argv[])
             #include "rhoEqn.H"
 
             // --- PIMPLE loop
-            for (pimple.start(); pimple.loop(); pimple++)
+            while (pimple.loop())
             {
                 #include "UEqn.H"
                 #include "YEqn.H"
                 #include "hsEqn.H"
 
-                // --- PISO loop
-                for (int corr=1; corr<=pimple.nCorr(); corr++)
+                // --- Pressure corrector loop
+                while (pimple.correct())
                 {
                     #include "pEqn.H"
                 }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 3bf304cecc..8031af20d5 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -15,7 +15,7 @@ if (pimple.transonic())
         )
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,12 +26,9 @@ if (pimple.transonic())
             parcels.Srho()
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -46,7 +43,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -57,12 +54,9 @@ else
             parcels.Srho()
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 94e551ca72..a46cb30f80 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -76,14 +76,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
index 96ad7dee52..50e24a70df 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
@@ -76,14 +76,14 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
index 0d095d3648..4bbaf792af 100644
--- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
+++ b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
@@ -68,14 +68,14 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "alphaEqn.H"
             #include "liftDragCoeffs.H"
             #include "UEqns.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
 
diff --git a/applications/solvers/multiphase/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H
index 453575b19b..268c2b4f0d 100644
--- a/applications/solvers/multiphase/bubbleFoam/pEqn.H
+++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H
@@ -42,7 +42,7 @@
         alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -51,12 +51,9 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index 159eb0f0bd..256a58d953 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -68,13 +68,15 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        for (pimple.start(); pimple.loop(); pimple++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
             #include "rhoEqn.H"
             #include "gammaPsi.H"
             #include "UEqn.H"
 
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index c604539cde..fb8893be2c 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -1,5 +1,5 @@
 {
-    if (pimple.nOuterCorr() == 1)
+    if (pimple.nCorrPIMPLE() == 1)
     {
         p =
         (
@@ -26,7 +26,7 @@
 
     #include "resetPhivPatches.H"
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -37,12 +37,9 @@
           - fvm::laplacian(rAUf, p)
         );
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phiv += (phiGradp + pEqn.flux())/rhof;
         }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
index 9161c80563..2d70a5bab5 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
@@ -30,7 +30,8 @@
         #include "alphaEqns.H"
     }
 
-    if (pimple.corr() == 0)
+    // correct interface on first PIMPLE corrector
+    if (pimple.corr() == 1)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
index 532d2bcc05..ade8af0081 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
@@ -32,7 +32,8 @@
         #include "alphaEqns.H"
     }
 
-    if (pimple.corr() == 0)
+    // correct interface on first PIMPLE corrector
+    if (pimple.corr() == 1)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 4b3b5127d4..1ac1596c4d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -116,8 +116,8 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        // --- Outer-corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
             #include "alphaEqnsSubCycle.H"
 
@@ -125,8 +125,8 @@ int main(int argc, char *argv[])
 
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
index c3b3726d12..1d7b9ca624 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
@@ -42,7 +42,7 @@
 
     adjustPhi(phi, U, pcorr);
 
-    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -51,7 +51,7 @@
 
         pcorrEqn.solve();
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index cdc541f0cd..26666c4120 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -39,7 +39,7 @@
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqnIncomp
         (
@@ -55,10 +55,10 @@
             )
            *p_rghEqnComp()
           + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
index 38bacbf4ea..1a718434a8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
@@ -4,7 +4,7 @@
 
     label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
-    if (nAlphaSubCycles > 1 && pimple.nOuterCorr() != 1)
+    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
     {
         FatalErrorIn(args.executable())
             << "Sub-cycling alpha is only allowed for PISO, "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 02d6764ae2..4ad1b3d01d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -74,8 +74,8 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        // --- Outer-corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
             #include "alphaEqnsSubCycle.H"
 
@@ -83,8 +83,8 @@ int main(int argc, char *argv[])
 
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 28c42abe79..035e8e237d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -39,7 +39,7 @@
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqnIncomp
         (
@@ -55,10 +55,10 @@
             )
            *p_rghEqnComp()
           + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
diff --git a/applications/solvers/multiphase/compressibleInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
index 954fcb27f2..87a055d641 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
@@ -1,19 +1,13 @@
-   #include "readTimeControls.H"
+    #include "readTimeControls.H"
 
-    label nAlphaCorr
-    (
-        readLabel(pimple.dict().lookup("nAlphaCorr"))
-    );
+    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
 
-    label nAlphaSubCycles
-    (
-        readLabel(pimple.dict().lookup("nAlphaSubCycles"))
-    );
+    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
-    if (nAlphaSubCycles > 1 && pimple.nOuterCorr() != 1)
+    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
     {
         FatalErrorIn(args.executable())
-            << "Sub-cycling alpha is only allowed for PISO, "
+            << "Sub-cycling alpha is only allowed for PISO operation, "
                "i.e. when the number of outer-correctors = 1"
             << exit(FatalError);
     }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
index c764472cb7..d29da482ba 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "alphaEqn.H"
             #include "kEpsilon.H"
@@ -85,8 +85,8 @@ int main(int argc, char *argv[])
             #include "TEqns.H"
             #include "UEqns.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index ea3e825448..d6af50988f 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -69,7 +69,7 @@
       + alpha2f*rAlphaAU2f/fvc::interpolate(rho2))
     );
 
-    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
         (
@@ -82,12 +82,11 @@
             (
                 (alpha1/rho1)*pEqnComp1()
               + (alpha2/rho2)*pEqnComp2()
-            ) +
-            pEqnIncomp,
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+            )
+          + pEqnIncomp,
+            mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             surfaceScalarField mSfGradp = pEqnIncomp.flux()/Dp;
             phi1 += rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 938ccbcc02..06811c2191 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -81,12 +81,12 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index b4f87264b9..042a17e1f3 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -84,12 +84,12 @@ int main(int argc, char *argv[])
         #include "zonePhaseVolumes.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index d281e177c8..6fc51c7cee 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -19,7 +19,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -28,12 +28,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/correctPhi.H
index cb3503e7b8..a2a61cef21 100644
--- a/applications/solvers/multiphase/interFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/correctPhi.H
@@ -34,7 +34,7 @@
 
     adjustPhi(phi, U, pcorr);
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -44,7 +44,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 1face6ddd5..a0c4827633 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -116,12 +116,12 @@ int main(int argc, char *argv[])
         #include "alphaEqnSubCycle.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 4b90d99499..d94e075424 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -21,7 +21,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -30,12 +30,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index d3895c142c..559916f2ce 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -84,12 +84,12 @@ int main(int argc, char *argv[])
         #include "alphaEqnSubCycle.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index d75a76766a..a88b5627e2 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -18,7 +18,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -27,12 +27,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 251a36c997..891736c77f 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -85,12 +85,12 @@ int main(int argc, char *argv[])
         #include "alphaEqnSubCycle.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
index 73cbee3ae4..d4b1194ba1 100644
--- a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
@@ -78,12 +78,12 @@ int main(int argc, char *argv[])
         #define twoPhaseProperties threePhaseProperties
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
index b68160e62d..6e02524e26 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
@@ -7,7 +7,7 @@ surfaceScalarField rhoPhi
         mesh
     ),
     mesh,
-    dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0)
+    dimensionedScalar("0", dimMass/dimTime, 0)
 );
 
 {
@@ -40,7 +40,7 @@ surfaceScalarField rhoPhi
         #include "alphaEqn.H"
     }
 
-    if (pimple.nOuterCorr() == 1)
+    if (pimple.nCorrPIMPLE() == 1)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index 70fd424e62..5ed807bbad 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -85,12 +85,12 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 6cfa59cc7c..c9b65eb9d5 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -13,8 +13,9 @@
 
     adjustPhi(phiU, U, p_rgh);
 
-    phi = phiU +
-        (
+    phi =
+        phiU
+      + (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf();
@@ -23,7 +24,7 @@
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -33,12 +34,9 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
index 4bfc58978e..a85e136eaa 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
@@ -39,6 +39,15 @@ if (mesh.nInternalFaces())
         fvc::surfaceSum(mag(phi))().internalField()
     );
 
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        sumPhi = max
+        (
+            sumPhi,
+            fvc::surfaceSum(mag(iter().phi()))().internalField()
+        );
+    }
+
     CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 
     meanCoNum =
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H
deleted file mode 100644
index 9bd18a9fcc..0000000000
--- a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H
+++ /dev/null
@@ -1,12 +0,0 @@
-#   include "CourantNo.H"
-
-// {
-//     scalar UrCoNum = 0.5*gMax
-//     (
-//         fvc::surfaceSum(mag(phi1 - phi2))().internalField()/mesh.V().field()
-//     )*runTime.deltaTValue();
-
-//     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
-
-//     CoNum = max(CoNum, UrCoNum);
-// }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
index 92c8676225..78f93384b0 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
@@ -34,7 +34,7 @@
 
     adjustPhi(phi, U, pcorr);
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -44,7 +44,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
index 6b859f068f..484ceae4ac 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
@@ -42,7 +42,7 @@
         dimensionedScalar("phi", dimArea*dimVelocity, 0)
     );
 
-    multiphaseSystem fluid(mesh, phi);
+    multiphaseSystem fluid(U, phi);
 
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -77,9 +77,8 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
 
-    singlePhaseTransportModel laminarTransport(U, phi);
 
     autoPtr<incompressible::LESModel> sgsModel
     (
-        incompressible::LESModel::New(U, phi, laminarTransport)
+        incompressible::LESModel::New(U, phi, fluid)
     );
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files
index 43869fde40..ff81afc634 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files
@@ -8,6 +8,7 @@ dragModels/Gibilaro/Gibilaro.C
 dragModels/WenYu/WenYu.C
 dragModels/SyamlalOBrien/SyamlalOBrien.C
 dragModels/blended/blended.C
+dragModels/interface/interface.C
 
 heatTransferModels/heatTransferModel/heatTransferModel.C
 heatTransferModels/heatTransferModel/newHeatTransferModel.C
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C
new file mode 100644
index 0000000000..0ef07da220
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interface.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+    defineTypeNameAndDebug(interface, 0);
+
+    addToRunTimeSelectionTable
+    (
+        dragModel,
+        interface,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dragModels::interface::interface
+(
+    const dictionary& interfaceDict,
+    const phaseModel& phase1,
+    const phaseModel& phase2
+)
+:
+    dragModel(interfaceDict, phase1, phase2)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dragModels::interface::~interface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::dragModels::interface::K
+(
+    const volScalarField& Ur
+) const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "K",
+                Ur.mesh().time().timeName(),
+                Ur.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            Ur.mesh(),
+            dimensionedScalar("K", dimDensity/dimTime, 0)
+        )
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H
new file mode 100644
index 0000000000..bc7c8223cf
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::dragModels::interface
+
+Description
+    Drag between phase separated by a VoF resolved interface.
+
+SourceFiles
+    interface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interface_H
+#define interface_H
+
+#include "dragModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class interface Declaration
+\*---------------------------------------------------------------------------*/
+
+class interface
+:
+    public dragModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("interface");
+
+
+    // Constructors
+
+        //- Construct from components
+        interface
+        (
+            const dictionary& interfaceDict,
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        );
+
+
+    //- Destructor
+    virtual ~interface();
+
+
+    // Member Functions
+
+        tmp<volScalarField> K(const volScalarField& Ur) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace dragModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index b4cc579835..4ca01f1a3a 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "correctPhi.H"
-    #include "CourantNos.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,25 +68,26 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "CourantNos.H"
+        #include "CourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             sgsModel->correct();
             fluid.solve();
             rho = fluid.rho();
+            #include "zonePhaseVolumes.H"
 
             //#include "interfacialCoeffs.H"
             //#include "TEqns.H"
             #include "UEqns.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
new file mode 100644
index 0000000000..635b3d21f8
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiphaseFixedFluxPressureFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    phi0Name_("phi0"),
+    phiName_("phi"),
+    rhoName_("rho")
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
+    phi0Name_(ptf.phi0Name_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_)
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    phi0Name_(dict.lookupOrDefault<word>("phi0", "phi0")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+{
+    if (dict.found("gradient"))
+    {
+        gradient() = scalarField("gradient", dict, p.size());
+        fixedGradientFvPatchScalarField::updateCoeffs();
+        fixedGradientFvPatchScalarField::evaluate();
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(patchInternalField());
+        gradient() = 0.0;
+    }
+}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf
+)
+:
+    fixedGradientFvPatchScalarField(wbppsf),
+    phi0Name_(wbppsf.phi0Name_),
+    phiName_(wbppsf.phiName_),
+    rhoName_(wbppsf.rhoName_)
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(wbppsf, iF),
+    phi0Name_(wbppsf.phi0Name_),
+    phiName_(wbppsf.phiName_),
+    rhoName_(wbppsf.rhoName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const surfaceScalarField& phi0 =
+        db().lookupObject<surfaceScalarField>(phi0Name_);
+
+    const surfaceScalarField& phi =
+        db().lookupObject<surfaceScalarField>(phiName_);
+
+    fvsPatchField<scalar> phi0p =
+        patch().patchField<surfaceScalarField, scalar>(phi0);
+
+    fvsPatchField<scalar> phip =
+        patch().patchField<surfaceScalarField, scalar>(phi);
+
+    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+    {
+        const fvPatchField<scalar>& rhop =
+            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+        phip /= rhop;
+    }
+
+    const fvsPatchField<scalar>& Dpp =
+        patch().lookupPatchField<surfaceScalarField, scalar>("Dp");
+
+    gradient() = (phi0p - phip)/patch().magSf()/Dpp;
+
+    fixedGradientFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>(os, "phi0", "phi0", phi0Name_);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    gradient().writeEntry("gradient", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        multiphaseFixedFluxPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
new file mode 100644
index 0000000000..f1699309b1
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multiphaseFixedFluxPressureFvPatchScalarField
+
+Description
+    Foam::multiphaseFixedFluxPressureFvPatchScalarField
+
+SourceFiles
+    multiphaseFixedFluxPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H
+#define multiphaseFixedFluxPressureFvPatchScalarFields_H
+
+#include "fvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class multiphaseFixedFluxPressureFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseFixedFluxPressureFvPatchScalarField
+:
+    public fixedGradientFvPatchScalarField
+{
+    // Private data
+
+        //- Name of the predicted flux transporting the field
+        word phi0Name_;
+
+        //- Name of the flux transporting the field
+        word phiName_;
+
+        //- Name of the density field used to normalise the mass flux
+        //  if neccessary
+        word rhoName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseFixedFluxPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  multiphaseFixedFluxPressureFvPatchScalarField onto a new patch
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new multiphaseFixedFluxPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new multiphaseFixedFluxPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 3382504153..2403597d9c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -441,24 +441,15 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 
 Foam::multiphaseSystem::multiphaseSystem
 (
-    const fvMesh& mesh,
+    const volVectorField& U,
     const surfaceScalarField& phi
 )
 :
-    IOdictionary
-    (
-        IOobject
-        (
-            "transportProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-    phases_(lookup("phases"), phaseModel::iNew(mesh)),
+    transportModel(U, phi),
 
-    mesh_(mesh),
+    phases_(lookup("phases"), phaseModel::iNew(U.mesh())),
+
+    mesh_(U.mesh()),
     phi_(phi),
 
     alphas_
@@ -514,7 +505,6 @@ Foam::multiphaseSystem::multiphaseSystem
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 {
     PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
@@ -530,6 +520,21 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
+{
+    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+
+    tmp<volScalarField> tnu = iter()*iter().nu();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        tnu() += iter()*iter().nu();
+    }
+
+    return tnu;
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
 (
     const phaseModel& phase
@@ -770,19 +775,81 @@ void Foam::multiphaseSystem::solve()
 
     label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
-    volScalarField& alpha = phases_.first();
-
     if (nAlphaSubCycles > 1)
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();
 
+        PtrList<volScalarField> alpha0s(phases_.size());
+        PtrList<surfaceScalarField> phiSums(phases_.size());
+
+        int phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            alpha0s.set
+            (
+                phasei,
+                new volScalarField(alpha.oldTime())
+            );
+
+            phiSums.set
+            (
+                phasei,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "phiSum" + alpha.name(),
+                        runTime.timeName(),
+                        mesh_
+                    ),
+                    mesh_,
+                    dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
+                )
+            );
+
+            phasei++;
+        }
+
         for
         (
-            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            subCycleTime alphaSubCycle
+            (
+                const_cast<Time&>(runTime),
+                nAlphaSubCycles
+            );
             !(++alphaSubCycle).end();
         )
         {
             solveAlphas();
+
+            int phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+            {
+                phiSums[phasei] += (runTime.deltaT()/totalDeltaT)*iter().phi();
+                phasei++;
+            }
+        }
+
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            phase.phi() = phiSums[phasei];
+
+            // Correct the time index of the field
+            // to correspond to the global time
+            alpha.timeIndex() = runTime.timeIndex();
+
+            // Reset the old-time field value
+            alpha.oldTime() = alpha0s[phasei];
+            alpha.oldTime().timeIndex() = runTime.timeIndex();
+
+            phasei++;
         }
     }
     else
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index 408fdcdafa..4c3e52e0ea 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -61,7 +61,7 @@ namespace Foam
 
 class multiphaseSystem
 :
-    public IOdictionary
+    public transportModel
 {
 
 public:
@@ -247,7 +247,7 @@ public:
         //- Construct from components
         multiphaseSystem
         (
-            const fvMesh& mesh,
+            const volVectorField& U,
             const surfaceScalarField& phi
         );
 
@@ -274,6 +274,9 @@ public:
         //- Return the mixture density
         tmp<volScalarField> rho() const;
 
+        //- Return the mixture laminar viscosity
+        tmp<volScalarField> nu() const;
+
         //- Return the virtual-mass coefficient for the given phase
         tmp<volScalarField> Cvm(const phaseModel& phase) const;
 
@@ -305,6 +308,10 @@ public:
         //- Solve for the mixture phase-fractions
         void solve();
 
+        //- Dummy correct
+        void correct()
+        {}
+
         //- Read base transportProperties dictionary
         bool read();
 };
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 6cbfbf2d16..98fa7b139e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -41,6 +41,20 @@
         phasei++;
     }
 
+    surfaceScalarField phi0
+    (
+        IOobject
+        (
+            "phi0",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("phi0", dimArea*dimVelocity, 0)
+    );
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -60,7 +74,8 @@
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
         );
         mrfZones.relativeFlux(phase.phi());
-        phase.phi() += rAlphaAUfs[phasei]*(g & mesh.Sf());
+        surfaceScalarField pphi0("pphi0", phase.phi());
+        pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf());
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();
@@ -99,13 +114,16 @@
                     phasej++;
                 }
 
-                phase.phi() +=
+                pphi0 +=
                     fvc::interpolate
                     ((1.0/phase.rho())*rAUs[phasei]*(*dcIter()))
                    *phi0s[phasej];
             }
         }
 
+        phi0 += alphafs[phasei]*pphi0;
+        phase.phi() = pphi0;
+
         phasei++;
     }
 
@@ -132,13 +150,13 @@
         phasei++;
     }
     Dp = mag(Dp);
-    adjustPhi(phi, U, p);
+    adjustPhi(phi0, U, p);
 
-    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
         (
-            fvc::div(phi)
+            fvc::div(phi0)
           - fvm::laplacian(Dp, p)
         );
 
@@ -149,10 +167,10 @@
             //   + (alpha2/rho2)*pEqnComp2()
             // ) +
             pEqnIncomp,
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+            mesh.solver(p.select(pimple.finalInnerIter()))
         );
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             surfaceScalarField mSfGradp = pEqnIncomp.flux()/Dp;
 
@@ -169,6 +187,7 @@
             }
 
             // dgdt =
+
             // (
             //     pos(alpha2)*(pEqnComp2 & p)/rho2
             //   - pos(alpha1)*(pEqnComp1 & p)/rho1
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H
new file mode 100644
index 0000000000..5f8f3fb19a
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H
@@ -0,0 +1,26 @@
+{
+    const scalarField& V = mesh.V();
+
+    forAll(mesh.cellZones(), czi)
+    {
+        const labelList& cellLabels = mesh.cellZones()[czi];
+
+        forAllConstIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            const volScalarField& alpha = iter();
+            scalar phaseVolume = 0;
+
+            forAll(cellLabels, cli)
+            {
+                label celli = cellLabels[cli];
+                phaseVolume += alpha[celli]*V[celli];
+            }
+
+            reduce(phaseVolume, sumOp<scalar>());
+
+            Info<< alpha.name()
+                << " phase volume in zone " << mesh.cellZones()[czi].name()
+                << " = " << phaseVolume*1e6 << " ml " << endl;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index 6769a47474..0179b98e68 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -76,12 +76,12 @@ int main(int argc, char *argv[])
         #include "zonePhaseVolumes.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index 9975ebac83..efb181f3f2 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -14,13 +14,14 @@
 
     adjustPhi(phiU, U, p_rgh);
 
-    phi = phiU +
-    (
-        mixture.surfaceTensionForce()
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+    phi =
+        phiU
+      + (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf();
 
-    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -29,12 +30,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index 3575c25cda..9f3fec616b 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -73,12 +73,12 @@ int main(int argc, char *argv[])
         rho = mixture.rho();
 
          // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index 78bfc89bb4..e7e5bc14cf 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -19,7 +19,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -28,12 +28,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 7b62411d2d..d89fd3ef4e 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -17,7 +17,7 @@ phi =
 surfaceScalarField phiU("phiU", phi);
 phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
 
-for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
@@ -25,9 +25,10 @@ for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     );
 
     p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
-    p_rghEqn.solve();
 
-    if (nonOrth == pimple.nNonOrthCorr())
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
     {
         phi -= p_rghEqn.flux();
     }
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index 76f141bbea..3d9fe516d8 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -65,7 +65,8 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        for (pimple.start(); pimple.loop(); pimple++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
             #include "calcVdj.H"
 
@@ -75,8 +76,8 @@ int main(int argc, char *argv[])
 
             #include "correctViscosity.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index 01b6442e9d..ac7fc35f68 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -14,7 +14,7 @@
 
     phi = phiU - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -23,12 +23,9 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve
-        (
-            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 465020a700..8a2b20e153 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -67,7 +67,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             twoPhaseProperties.correct();
 
@@ -75,8 +75,8 @@ int main(int argc, char *argv[])
 
             #include "UEqn.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 53d11f9e06..1788dd2808 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -58,7 +58,7 @@
         alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
     );
 
-    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -67,12 +67,9 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        pEqn.solve
-        (
-            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        );
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (nonOrth == pimple.nNonOrthCorr())
+        if (pimple.finalNonOrthogonalIter())
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 43f353cd1c..19810a15f1 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -76,14 +76,14 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (pimple.start(); pimple.loop(); pimple++)
+        while (pimple.loop())
         {
             #include "alphaEqn.H"
             #include "liftDragCoeffs.H"
             #include "UEqns.H"
 
-            // --- PISO loop
-            for (int corr=0; corr<pimple.nCorr(); corr++)
+            // --- Pressure corrector loop
+            while (pimple.correct())
             {
                 #include "pEqn.H"
 
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
index 9f256bd689..c8e3ccba56 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
@@ -12,18 +12,80 @@
         )
     );
 
-    dimensionedScalar rho(mechanicalProperties.lookup("rho"));
-    dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
-    dimensionedScalar nu(mechanicalProperties.lookup("nu"));
+    const dictionary& rhoDict(mechanicalProperties.subDict("rho"));
+    word rhoType(rhoDict.lookup("rho"));
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimMass/dimVolume, 0.0)
+    );
+
+    if (rhoType == "rhoInf")
+    {
+        rho = rhoDict.lookup("rhoInf");
+    }
+
+    volScalarField rhoE
+    (
+        IOobject
+        (
+            "E",
+            runTime.timeName(0),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimMass/dimLength/sqr(dimTime), 0.0)
+    );
+
+    const dictionary& EDict(mechanicalProperties.subDict("E"));
+    word EType(EDict.lookup("E"));
+    if (EType == "EInf")
+    {
+        rhoE = EDict.lookup("EInf");
+    }
+
+
+    volScalarField nu
+    (
+        IOobject
+        (
+            "nu",
+            runTime.timeName(0),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless, 0.0)
+    );
+
+    const dictionary& nuDict(mechanicalProperties.subDict("nu"));
+    word nuType(nuDict.lookup("nu"));
+
+    if (nuType == "nuInf")
+    {
+        nu = nuDict.lookup("nuInf");
+    }
 
     Info<< "Normalising E : E/rho\n" << endl;
-    dimensionedScalar E = rhoE/rho;
+    volScalarField E = rhoE/rho;
 
     Info<< "Calculating Lame's coefficients\n" << endl;
 
-    dimensionedScalar mu = E/(2.0*(1.0 + nu));
-    dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
-    dimensionedScalar threeK = E/(1.0 - 2.0*nu);
+    volScalarField mu = E/(2.0*(1.0 + nu));
+    volScalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
+    volScalarField threeK = E/(1.0 - 2.0*nu);
 
     Switch planeStress(mechanicalProperties.lookup("planeStress"));
 
@@ -38,7 +100,3 @@
     {
         Info<< "Plane Strain\n" << endl;
     }
-
-    Info<< "mu = " << mu.value() << " Pa/rho\n";
-    Info<< "lambda = " << lambda.value() << " Pa/rho\n";
-    Info<< "threeK = " << threeK.value() << " Pa/rho\n";
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
index 12e10a607b..4d7bb43fc9 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
@@ -1,46 +1,122 @@
-    Info<< "Reading thermal properties\n" << endl;
+Info<< "Reading thermal properties\n" << endl;
 
-    IOdictionary thermalProperties
+IOdictionary thermalProperties
+(
+    IOobject
+    (
+        "thermalProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch thermalStress(thermalProperties.lookup("thermalStress"));
+
+volScalarField threeKalpha
+(
+    IOobject
+    (
+        "threeKalpha",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0)
+);
+
+
+volScalarField DT
+(
+    IOobject
+    (
+        "DT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimensionSet(0, 2, -1 , 0, 0), 0.0)
+);
+
+
+if (thermalStress)
+{
+    volScalarField C
     (
         IOobject
         (
-            "thermalProperties",
-            runTime.constant(),
+            "C",
+            runTime.timeName(0),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE
-        )
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0)
     );
 
-    Switch thermalStress(thermalProperties.lookup("thermalStress"));
-
-    dimensionedScalar threeKalpha
-    (
-        "threeKalpha",
-        dimensionSet(0, 2, -2 , -1, 0),
-        0
-    );
-
-    dimensionedScalar DT
-    (
-        "DT",
-        dimensionSet(0, 2, -1 , 0, 0),
-        0
-    );
-
-    if (thermalStress)
+    const dictionary& CDict(thermalProperties.subDict("C"));
+    word CType(CDict.lookup("C"));
+    if (CType == "CInf")
     {
-        dimensionedScalar C(thermalProperties.lookup("C"));
-        dimensionedScalar rhoK(thermalProperties.lookup("k"));
-        dimensionedScalar alpha(thermalProperties.lookup("alpha"));
-
-        Info<< "Normalising k : k/rho\n" << endl;
-        dimensionedScalar k = rhoK/rho;
-
-        Info<< "Calculating thermal coefficients\n" << endl;
-
-        threeKalpha = threeK*alpha;
-        DT.value() = (k/C).value();
-
-        Info<< "threeKalpha = " << threeKalpha.value() << " Pa/rho\n";
+        C = CDict.lookup("CInf");
     }
+
+
+    volScalarField rhoK
+    (
+        IOobject
+        (
+            "k",
+            runTime.timeName(0),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, 1, -3 , -1, 0), 0.0)
+    );
+
+    const dictionary& kDict(thermalProperties.subDict("k"));
+    word kType(kDict.lookup("k"));
+    if (kType == "kInf")
+    {
+        rhoK = kDict.lookup("kInf");
+    }
+
+    volScalarField alpha
+    (
+        IOobject
+        (
+            "alpha",
+            runTime.timeName(0),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(0, 0, 0 , -1, 0), 0.0)
+    );
+
+    const dictionary& alphaDict(thermalProperties.subDict("alpha"));
+    word alphaType(alphaDict.lookup("alpha"));
+
+    if (alphaType == "alphaInf")
+    {
+        alpha = alphaDict.lookup("alphaInf");
+    }
+
+    Info<< "Normalising k : k/rho\n" << endl;
+    volScalarField k = rhoK/rho;
+
+    Info<< "Calculating thermal coefficients\n" << endl;
+
+    threeKalpha = threeK*alpha;
+    DT = k/C;
+
+}
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C b/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
index 85f3c67b38..b24239c0f8 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
@@ -149,14 +149,20 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs()
     const dictionary& thermalProperties =
         db().lookupObject<IOdictionary>("thermalProperties");
 
-    dimensionedScalar rho(mechanicalProperties.lookup("rho"));
-    dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
-    dimensionedScalar nu(mechanicalProperties.lookup("nu"));
 
-    dimensionedScalar E = rhoE/rho;
-    dimensionedScalar mu = E/(2.0*(1.0 + nu));
-    dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
-    dimensionedScalar threeK = E/(1.0 - 2.0*nu);
+    const fvPatchField<scalar>& rho =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const fvPatchField<scalar>& rhoE =
+        patch().lookupPatchField<volScalarField, scalar>("E");
+
+    const fvPatchField<scalar>& nu =
+        patch().lookupPatchField<volScalarField, scalar>("nu");
+
+    scalarField E = rhoE/rho;
+    scalarField mu = E/(2.0*(1.0 + nu));
+    scalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
+    scalarField threeK = E/(1.0 - 2.0*nu);
 
     Switch planeStress(mechanicalProperties.lookup("planeStress"));
 
@@ -166,7 +172,7 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs()
         threeK = E/(1.0 - nu);
     }
 
-    scalar twoMuLambda = (2*mu + lambda).value();
+    scalarField twoMuLambda = (2*mu + lambda);
 
     vectorField n(patch().nf());
 
@@ -175,7 +181,7 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs()
 
     gradient() =
     (
-        (traction_ + pressure_*n)/rho.value()
+        (traction_ + pressure_*n)/rho
       + twoMuLambda*fvPatchField<vector>::snGrad() - (n & sigmaD)
     )/twoMuLambda;
 
@@ -183,13 +189,13 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs()
 
     if (thermalStress)
     {
-        dimensionedScalar alpha(thermalProperties.lookup("alpha"));
-        dimensionedScalar threeKalpha = threeK*alpha;
+        const fvPatchField<scalar>&  threeKalpha=
+            patch().lookupPatchField<volScalarField, scalar>("threeKalpha");
 
         const fvPatchField<scalar>& T =
             patch().lookupPatchField<volScalarField, scalar>("T");
 
-        gradient() += n*threeKalpha.value()*T/twoMuLambda;
+        gradient() += n*threeKalpha*T/twoMuLambda;
     }
 
     fixedGradientFvPatchVectorField::updateCoeffs();
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
index 9c65f5cd82..fa95c852c1 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
@@ -34,7 +34,7 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+#include "argList.H"
 #include "conformalVoronoiMesh.H"
 
 using namespace Foam;
diff --git a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
index 1b1d317138..2f56ee07de 100644
--- a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
+++ b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
@@ -44,7 +44,9 @@ int main(int argc, char *argv[])
     argList::validArgs.append("FOAMThermodynamicsFile");
     argList args(argc, argv);
 
-    chemkinReader cr(args[1], args[2]);
+    speciesTable species;
+
+    chemkinReader cr(args[1], species, args[2]);
 
     OFstream reactionsFile(args[3]);
     reactionsFile
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 995c75cdc0..3211df1514 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -655,6 +655,27 @@ Foam::scalar Foam::face::sweptVol
     const pointField& newPoints
 ) const
 {
+    if (size() == 3)
+    {
+        return
+        (
+            triPointRef
+            (
+                oldPoints[operator[](0)],
+                oldPoints[operator[](1)],
+                oldPoints[operator[](2)]
+            ).sweptVol
+            (
+                triPointRef
+                (
+                    newPoints[operator[](0)],
+                    newPoints[operator[](1)],
+                    newPoints[operator[](2)]
+                )
+            )
+        );
+    }
+
     scalar sv = 0;
 
     // Calculate the swept volume by breaking the face into triangles and
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index cc2f3e4225..d30f23b9bd 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -280,6 +280,7 @@ $(limitedSchemes)/limitWith/limitWith.C
 multivariateSchemes = $(surfaceInterpolation)/multivariateSchemes
 $(multivariateSchemes)/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationSchemes.C
 $(multivariateSchemes)/multivariateSelectionScheme/multivariateSelectionSchemes.C
+$(multivariateSchemes)/multivariateIndependentScheme/multivariateIndependentSchemes.C
 $(multivariateSchemes)/upwind/multivariateUpwind.C
 $(multivariateSchemes)/Gamma/multivariateGamma.C
 $(multivariateSchemes)/vanLeer/multivariateVanLeer.C
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index 87c3f13090..98a742e16b 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -42,8 +42,8 @@ void Foam::pimpleControl::read()
 
     // Read solution controls
     const dictionary& pimpleDict = dict();
-    nOuterCorr_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
-    nCorr_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
+    nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
     turbOnFinalIterOnly_ =
         pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
 }
@@ -51,12 +51,14 @@ void Foam::pimpleControl::read()
 
 bool Foam::pimpleControl::criteriaSatisfied()
 {
-    if ((corr_ == 0) || residualControl_.empty() || finalIter())
+    // no checks on first iteration - nothing has been calculated yet
+    if ((corr_ == 1) || residualControl_.empty() || finalIter())
     {
         return false;
     }
 
-    bool firstIter = corr_ == 1;
+
+    bool storeIni = this->storeInitialResiduals();
 
     bool achieved = true;
     bool checked = false;    // safety that some checks were indeed performed
@@ -73,7 +75,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
 
             checked = true;
 
-            if (firstIter)
+            if (storeIni)
             {
                 residualControl_[fieldI].initialResidual =
                     sp.first().initialResidual();
@@ -83,7 +85,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
             bool relCheck = false;
 
             scalar relative = 0.0;
-            if (!firstIter)
+            if (!storeIni)
             {
                 const scalar iniRes =
                     residualControl_[fieldI].initialResidual
@@ -97,9 +99,10 @@ bool Foam::pimpleControl::criteriaSatisfied()
 
             if (debug)
             {
-                Info<< algorithmName_ << "loop statistics:" << endl;
+                Info<< algorithmName_ << " loop:" << endl;
 
-                Info<< "    " << variableName << " iter " << corr_
+                Info<< "    " << variableName
+                    << " PIMPLE iter " << corr_
                     << ": ini res = "
                     << residualControl_[fieldI].initialResidual
                     << ", abs tol = " << residual
@@ -120,25 +123,25 @@ bool Foam::pimpleControl::criteriaSatisfied()
 Foam::pimpleControl::pimpleControl(fvMesh& mesh)
 :
     solutionControl(mesh, "PIMPLE"),
-    nOuterCorr_(0),
-    nCorr_(0),
-    corr_(0),
+    nCorrPIMPLE_(0),
+    nCorrPISO_(0),
+    corrPISO_(0),
     turbOnFinalIterOnly_(true)
 {
     read();
 
-    if (nOuterCorr_ > 1)
+    if (nCorrPIMPLE_ > 1)
     {
         Info<< nl;
         if (residualControl_.empty())
         {
             Info<< algorithmName_ << ": no residual control data found. "
-                << "Calculations will employ " << nOuterCorr_
+                << "Calculations will employ " << nCorrPIMPLE_
                 << " corrector loops" << nl << endl;
         }
         else
         {
-            Info<< algorithmName_ << ": max iterations = " << nOuterCorr_
+            Info<< algorithmName_ << ": max iterations = " << nCorrPIMPLE_
                 << endl;
             forAll(residualControl_, i)
             {
@@ -164,4 +167,60 @@ Foam::pimpleControl::~pimpleControl()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::pimpleControl::loop()
+{
+    read();
+
+    corr_++;
+
+    if (debug)
+    {
+        Info<< algorithmName_ << " loop: corr = " << corr_ << endl;
+    }
+
+    if (corr_ == nCorrPIMPLE_ + 1)
+    {
+        if ((!residualControl_.empty()) && (nCorrPIMPLE_ != 1))
+        {
+            Info<< algorithmName_ << ": not converged within "
+                << nCorrPIMPLE_ << " iterations" << endl;
+        }
+
+        corr_ = 0;
+        mesh_.data::remove("finalIteration");
+        return false;
+    }
+
+    bool completed = false;
+    if (criteriaSatisfied())
+    {
+        Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
+            << endl;
+        completed = true;
+    }
+    else
+    {
+        if (finalIter())
+        {
+            mesh_.data::add("finalIteration", true);
+        }
+
+        if (corr_ <= nCorrPIMPLE_)
+        {
+            if (nCorrPIMPLE_ != 1)
+            {
+                Info<< algorithmName_ << ": iteration " << corr_ << endl;
+                storePrevIterFields();
+            }
+
+            completed = false;
+        }
+    }
+
+    return !completed;
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index 3b891a5aa9..9a5912e4fb 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -55,13 +55,13 @@ protected:
         // Solution controls
 
             //- Maximum number of PIMPLE correctors
-            label nOuterCorr_;
+            label nCorrPIMPLE_;
 
             //- Maximum number of PISO correctors
-            label nCorr_;
+            label nCorrPISO_;
 
-            //- Current PIMPLE corrector
-            label corr_;
+            //- Current PISO corrector
+            label corrPISO_;
 
             //- Flag to indicate whether to only solve turbulence on final iter
             bool turbOnFinalIterOnly_;
@@ -105,41 +105,35 @@ public:
 
         // Access
 
-            //- Current corrector index
-            inline label corr() const;
-
             //- Maximum number of PIMPLE correctors
-            inline label nOuterCorr() const;
+            inline label nCorrPIMPLE() const;
 
             //- Maximum number of PISO correctors
-            inline label nCorr() const;
+            inline label nCorrPISO() const;
+
+            //- Current PISO corrector index
+            inline label corrPISO() const;
 
 
         // Solution control
 
-            //- Loop start
-            inline bool start();
+            //- PIMPLE loop
+            virtual bool loop();
 
-            //- Loop loop
-            inline bool loop();
+            //- Pressure corrector loop
+            inline bool correct();
+
+            //- Helper function to identify when to store the intial residuals
+            inline bool storeInitialResiduals() const;
 
             //- Helper function to identify final PIMPLE (outer) iteration
             inline bool finalIter() const;
 
             //- Helper function to identify final inner iteration
-            inline bool finalInnerIter
-            (
-                const label corr,
-                const label nonOrth
-            ) const;
+            inline bool finalInnerIter() const;
 
             //- Helper function to identify whether to solve for turbulence
             inline bool turbCorr() const;
-
-
-    // Member Operators
-
-        void operator++(int);
 };
 
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index 59e4e14204..1fe6098731 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -25,89 +25,64 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::pimpleControl::corr() const
+inline Foam::label Foam::pimpleControl::nCorrPIMPLE() const
 {
-    return corr_;
+    return nCorrPIMPLE_;
 }
 
 
-inline Foam::label Foam::pimpleControl::nOuterCorr() const
+inline Foam::label Foam::pimpleControl::nCorrPISO() const
 {
-    return nOuterCorr_;
+    return nCorrPISO_;
 }
 
 
-inline Foam::label Foam::pimpleControl::nCorr() const
+inline Foam::label Foam::pimpleControl::corrPISO() const
 {
-    return nCorr_;
+    return corrPISO_;
 }
 
 
-inline bool Foam::pimpleControl::start()
+inline bool Foam::pimpleControl::correct()
 {
-    corr_ = 0;
+    corrPISO_++;
 
-    return true;
-}
-
-
-inline bool Foam::pimpleControl::loop()
-{
-    read();
-
-    if (criteriaSatisfied())
+    if (debug)
     {
-        Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
-            << endl;
-        return false;
+        Info<< algorithmName_ << " correct: corrPISO = " << corrPISO_ << endl;
+    }
+
+    if (corrPISO_ <= nCorrPISO_)
+    {
+        return true;
     }
     else
     {
-        if (finalIter())
-        {
-            mesh_.data::add("finalIteration", true);
-        }
-
-        if (corr_ < nOuterCorr_)
-        {
-            if (nOuterCorr_ != 1)
-            {
-                Info<< algorithmName_ << ": iteration " << corr_ + 1 << endl;
-                storePrevIterFields();
-            }
-
-            return true;
-        }
-        else
-        {
-            if ((!residualControl_.empty()) && (nOuterCorr_ != 1))
-            {
-                Info<< algorithmName_ << ": not converged within "
-                    << nOuterCorr_ << " iterations" << endl;
-            }
-
-            return false;
-        }
+        corrPISO_ = 0;
+        return false;
     }
 }
 
 
+inline bool Foam::pimpleControl::storeInitialResiduals() const
+{
+    // start from second PIMPLE iteration
+    return (corr_ == 2) && (corrPISO_ == 0) && (corrNonOrtho_ == 0);
+}
+
+
 inline bool Foam::pimpleControl::finalIter() const
 {
-    return corr_ == nOuterCorr_-1;
+    return corr_ == nCorrPIMPLE_;
 }
 
 
-inline bool Foam::pimpleControl::finalInnerIter
-(
-    const label corr,
-    const label nonOrth
-) const
+inline bool Foam::pimpleControl::finalInnerIter() const
 {
     return
-        corr_ == nOuterCorr_-1
-     && corr == nCorr_-1
-     && nonOrth == nNonOrthCorr_;
+       corr_ == nCorrPIMPLE_
+    && corrPISO_ == nCorrPISO_
+    && corrNonOrtho_ == nNonOrthCorr_ + 1;
 }
 
 
@@ -117,15 +92,4 @@ inline bool Foam::pimpleControl::turbCorr() const
 }
 
 
-inline void Foam::pimpleControl::operator++(int)
-{
-    if (finalIter())
-    {
-        mesh_.data::remove("finalIteration");
-    }
-
-    corr_++;
-}
-
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
index d2b603a5e0..b944f38b8f 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
@@ -119,4 +119,38 @@ Foam::simpleControl::~simpleControl()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::simpleControl::loop()
+{
+    read();
+
+    Time& time = const_cast<Time&>(mesh_.time());
+
+    if (initialised_)
+    {
+        if (criteriaSatisfied())
+        {
+            Info<< nl << algorithmName_ << " solution converged in "
+                << time.timeName() << " iterations" << nl << endl;
+
+            // Set to finalise calculation
+            time.writeAndEnd();
+        }
+        else
+        {
+            storePrevIterFields();
+        }
+    }
+    else
+    {
+        initialised_ = true;
+        storePrevIterFields();
+    }
+
+
+    return time.loop();
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H
index ee718ce424..71f8d0e04c 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H
@@ -96,7 +96,7 @@ public:
         // Solution control
 
             //- Loop loop
-            inline bool loop();
+            virtual bool loop();
 };
 
 
@@ -106,10 +106,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "simpleControlI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index a238d94e0a..3de1587e7e 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -170,7 +170,9 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
     algorithmName_(algorithmName),
     nNonOrthCorr_(0),
     momentumPredictor_(true),
-    transonic_(false)
+    transonic_(false),
+    corr_(0),
+    corrNonOrtho_(0)
 {}
 
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
index 94a8b41305..45970fc325 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
@@ -82,6 +82,15 @@ protected:
             bool transonic_;
 
 
+        // Evolution
+
+            //- Current corrector loop index
+            label corr_;
+
+            //- Current non-orthogonal corrector loop index
+            label corrNonOrtho_;
+
+
     // Protected Member Functions
 
         //- Read controls from fvSolution dictionary
@@ -137,17 +146,35 @@ public:
             //- Return the solution dictionary
             inline const dictionary& dict() const;
 
+            //- Current corrector loop index
+            inline label corr() const;
+
+            //- Current non-orthogonal corrector index
+            inline label corrNonOrtho() const;
+
 
         // Solution control
 
             //- Maximum number of non-orthogonal correctors
             inline label nNonOrthCorr() const;
 
+            //- Helper function to identify final non-orthogonal iteration
+            inline bool finalNonOrthogonalIter() const;
+
             //- Flag to indicate to solve for momentum
             inline bool momentumPredictor() const;
 
             //- Flag to indicate to solve using transonic algorithm
             inline bool transonic() const;
+
+
+        // Evolution
+
+            //- Main control loop
+            virtual bool loop() = 0;
+
+            //- Non-orthogonal corrector loop
+            inline bool correctNonOrthogonal();
 };
 
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
index 438c8569d3..b9d975cfba 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
@@ -31,12 +31,30 @@ inline const Foam::dictionary& Foam::solutionControl::dict() const
 }
 
 
+inline Foam::label Foam::solutionControl::corr() const
+{
+    return corr_;
+}
+
+
+inline Foam::label Foam::solutionControl::corrNonOrtho() const
+{
+    return corrNonOrtho_;
+}
+
+
 inline Foam::label Foam::solutionControl::nNonOrthCorr() const
 {
     return nNonOrthCorr_;
 }
 
 
+inline bool Foam::solutionControl::finalNonOrthogonalIter() const
+{
+    return corrNonOrtho_ == nNonOrthCorr_ + 1;
+}
+
+
 inline bool Foam::solutionControl::momentumPredictor() const
 {
     return momentumPredictor_;
@@ -49,4 +67,26 @@ inline bool Foam::solutionControl::transonic() const
 }
 
 
+inline bool Foam::solutionControl::correctNonOrthogonal()
+{
+    corrNonOrtho_++;
+
+    if (debug)
+    {
+        Info<< algorithmName_ << " correctNonOrthogonal: corrNonOrtho = "
+            << corrNonOrtho_ << endl;
+    }
+
+    if (corrNonOrtho_ <= nNonOrthCorr_ + 1)
+    {
+        return true;
+    }
+    else
+    {
+        corrNonOrtho_ = 0;
+        return false;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.C
new file mode 100644
index 0000000000..9539fee912
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.C
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multivariateIndependentScheme.H"
+#include "limitedSurfaceInterpolationScheme.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "upwind.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::multivariateIndependentScheme<Type>::multivariateIndependentScheme
+(
+    const fvMesh& mesh,
+    const typename multivariateSurfaceInterpolationScheme<Type>::
+        fieldTable& fields,
+    const surfaceScalarField& faceFlux,
+    Istream& schemeData
+)
+:
+    multivariateSurfaceInterpolationScheme<Type>
+    (
+        mesh,
+        fields,
+        faceFlux,
+        schemeData
+    ),
+    schemes_(schemeData),
+    faceFlux_(faceFlux)
+{}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.H b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.H
new file mode 100644
index 0000000000..61a91dcc48
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentScheme.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multivariateIndependentScheme
+
+Description
+    Generic multi-variate discretisation scheme class for which any of the
+    NVD, CNVD or NVDV schemes may be selected for each variable and applied
+    independently.
+
+    This is equivalent to using separate "div" terms and schemes for each
+    variable/equation.
+
+SourceFiles
+    multivariateIndependentScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multivariateIndependentScheme_H
+#define multivariateIndependentScheme_H
+
+#include "multivariateSurfaceInterpolationScheme.H"
+#include "limitedSurfaceInterpolationScheme.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class multivariateIndependentScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class multivariateIndependentScheme
+:
+    public multivariateSurfaceInterpolationScheme<Type>
+{
+    // Private data
+
+        dictionary schemes_;
+        const surfaceScalarField& faceFlux_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multivariateIndependentScheme(const multivariateIndependentScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const multivariateIndependentScheme&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multivariateIndependent");
+
+
+    // Constructors
+
+        //- Construct for field, faceFlux and Istream
+        multivariateIndependentScheme
+        (
+            const fvMesh& mesh,
+            const typename multivariateSurfaceInterpolationScheme<Type>::
+                fieldTable& fields,
+            const surfaceScalarField& faceFlux,
+            Istream& schemeData
+        );
+
+
+    // Member Operators
+
+        tmp<surfaceInterpolationScheme<Type> > operator()
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field
+        ) const
+        {
+            return surfaceInterpolationScheme<Type>::New
+            (
+                faceFlux_.mesh(),
+                faceFlux_,
+                schemes_.lookup(field.name())
+            );
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "multivariateIndependentScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentSchemes.C
similarity index 65%
rename from src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H
rename to src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentSchemes.C
index 9080c9b245..1d6f297ada 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateIndependentScheme/multivariateIndependentSchemes.C
@@ -23,40 +23,26 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Time.H"
+#include "multivariateIndependentScheme.H"
+#include "volFields.H"
+#include "surfaceFields.H"
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inline bool Foam::simpleControl::loop()
+namespace Foam
 {
-    read();
 
-    Time& time = const_cast<Time&>(mesh_.time());
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    if (initialised_)
-    {
-        if (criteriaSatisfied())
-        {
-            Info<< nl << algorithmName_ << " solution converged in "
-                << time.timeName() << " iterations" << nl << endl;
+defineNamedTemplateTypeNameAndDebug(multivariateIndependentScheme<scalar>, 0);
 
-            // Set to finalise calculation
-            time.writeAndEnd();
-        }
-        else
-        {
-            storePrevIterFields();
-        }
-    }
-    else
-    {
-        initialised_ = true;
-        storePrevIterFields();
-    }
+multivariateSurfaceInterpolationScheme<scalar>::addIstreamConstructorToTable
+    <multivariateIndependentScheme<scalar> >
+    addMultivariateIndependentSchemeScalarConstructorToTable_;
 
 
-    return time.loop();
-}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
index 855db45c88..dfef9f4835 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VoidFraction.H"
-#
+
 // * * * * * * * * * * * * * Protectd Member Functions * * * * * * * * * * * //
 
 template<class CloudType>
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.C
index d6dd42a045..75a21a1690 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.C
@@ -31,7 +31,8 @@ template<class ThermoType>
 Foam::autoPtr<Foam::chemistryReader<ThermoType> >
 Foam::chemistryReader<ThermoType>::New
 (
-    const dictionary& thermoDict
+    const dictionary& thermoDict,
+    speciesTable& species
 )
 {
     // Let the chemistry reader type default to CHEMKIN
@@ -50,7 +51,7 @@ Foam::chemistryReader<ThermoType>::New
     {
         FatalErrorIn
         (
-            "chemistryReader::New(const dictionary& thermoDict)"
+            "chemistryReader::New(const dictionary&, speciesTable&)"
         )   << "Unknown chemistryReader type "
             << chemistryReaderTypeName << nl << nl
             << "Valid chemistryReader types are:" << nl
@@ -58,7 +59,10 @@ Foam::chemistryReader<ThermoType>::New
             << exit(FatalError);
     }
 
-    return autoPtr<chemistryReader<ThermoType> >(cstrIter()(thermoDict));
+    return autoPtr<chemistryReader<ThermoType> >
+    (
+        cstrIter()(thermoDict, species)
+    );
 }
 
 
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.H b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.H
index 839100b7a3..13313b9772 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.H
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.H
@@ -85,16 +85,21 @@ public:
             chemistryReader,
             dictionary,
             (
-                const dictionary& thermoDict
+                const dictionary& thermoDict,
+                speciesTable& species
             ),
-            (thermoDict)
+            (thermoDict, species)
         );
 
 
     // Selectors
 
         //- Select constructed from dictionary
-        static autoPtr<chemistryReader> New(const dictionary& thermoDict);
+        static autoPtr<chemistryReader> New
+        (
+            const dictionary& thermoDict,
+            speciesTable& species
+        );
 
 
     //- Destructor
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
index a95f9439c9..26fb1b1478 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
@@ -854,29 +854,31 @@ void Foam::chemkinReader::read
 Foam::chemkinReader::chemkinReader
 (
     const fileName& CHEMKINFileName,
+    speciesTable& species,
     const fileName& thermoFileName
 )
 :
     lineNo_(1),
     specieNames_(10),
-    speciesTable_(),
+    speciesTable_(species),
     reactions_(speciesTable_, speciesThermo_)
 {
     read(CHEMKINFileName, thermoFileName);
 }
 
 
-Foam::chemkinReader::chemkinReader(const dictionary& thermoDict)
+Foam::chemkinReader::chemkinReader
+(
+    const dictionary& thermoDict,
+    speciesTable& species
+)
 :
     lineNo_(1),
     specieNames_(10),
-    speciesTable_(),
+    speciesTable_(species),
     reactions_(speciesTable_, speciesThermo_)
 {
-    fileName chemkinFile
-    (
-        fileName(thermoDict.lookup("CHEMKINFile")).expand()
-    );
+    fileName chemkinFile(fileName(thermoDict.lookup("CHEMKINFile")).expand());
 
     fileName thermoFile = fileName::null;
 
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H
index f00889d056..6e0d4dc5c1 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H
@@ -193,7 +193,7 @@ private:
         HashTable<label> specieIndices_;
 
         //- Table of species
-        speciesTable speciesTable_;
+        speciesTable& speciesTable_;
 
         //- Specie phase
         HashTable<phase> speciePhase_;
@@ -318,11 +318,12 @@ public:
         chemkinReader
         (
             const fileName& chemkinFile,
+            speciesTable& species,
             const fileName& thermoFileName = fileName::null
         );
 
         //- Construct by getting the CHEMKIN III file name from dictionary
-        chemkinReader(const dictionary& thermoDict);
+        chemkinReader(const dictionary& thermoDict, speciesTable& species);
 
 
     //- Destructor
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.C
index 3024473270..0efe8e2b8a 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.C
@@ -27,12 +27,28 @@ License
 #include "IFstream.H"
 #include "addToRunTimeSelectionTable.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ThermoType>
+Foam::speciesTable& Foam::foamChemistryReader<ThermoType>::setSpecies
+(
+    const dictionary& dict,
+    speciesTable& species
+)
+{
+    wordList s(dict.lookup("species"));
+    species.transfer(s);
+    return species;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 template<class ThermoType>
 Foam::foamChemistryReader<ThermoType>::foamChemistryReader
 (
     const fileName& reactionsFileName,
+    speciesTable& species,
     const fileName& thermoFileName
 )
 :
@@ -51,8 +67,8 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
             fileName(thermoFileName).expand()
         )()
     ),
+    speciesTable_(setSpecies(chemDict_, species)),
     speciesThermo_(thermoDict_),
-    speciesTable_(chemDict_.lookup("species")),
     reactions_(speciesTable_, speciesThermo_, chemDict_)
 {}
 
@@ -60,7 +76,8 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
 template<class ThermoType>
 Foam::foamChemistryReader<ThermoType>::foamChemistryReader
 (
-    const dictionary& thermoDict
+    const dictionary& thermoDict,
+    speciesTable& species
 )
 :
     chemistryReader<ThermoType>(),
@@ -79,7 +96,7 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
         )()
     ),
     speciesThermo_(thermoDict_),
-    speciesTable_(chemDict_.lookup("species")),
+    speciesTable_(setSpecies(chemDict_, species)),
     reactions_(speciesTable_, speciesThermo_, chemDict_)
 {}
 
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.H b/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.H
index 84f18d2221..925692dc7a 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.H
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.H
@@ -67,7 +67,7 @@ class foamChemistryReader
         HashPtrTable<ThermoType> speciesThermo_;
 
         //- Table of species
-        speciesTable speciesTable_;
+        speciesTable& speciesTable_;
 
         //- List of the reactions
         ReactionList<ThermoType> reactions_;
@@ -75,6 +75,9 @@ class foamChemistryReader
 
     // Private Member Functions
 
+        //- Set the species list
+        speciesTable& setSpecies(const dictionary& dict, speciesTable& species);
+
         //- Disallow default bitwise copy construct
         foamChemistryReader(const foamChemistryReader&);
 
@@ -94,12 +97,17 @@ public:
         foamChemistryReader
         (
             const fileName& reactionsFileName,
+            speciesTable& species,
             const fileName& thermoFileName
         );
 
         //- Construct by getting the foamChemistry and thermodynamics file names
         //  from dictionary
-        foamChemistryReader(const dictionary& thermoDict);
+        foamChemistryReader
+        (
+            const dictionary& thermoDict,
+            speciesTable& species
+        );
 
 
     //- Destructor
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
index d5707fe8fd..3d78443860 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
@@ -49,7 +49,7 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
 template<class ThermoType>
 void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
 {
-        // It changes Yt patches to "calculated"
+    // It changes Yt patches to "calculated"
     volScalarField Yt("Yt", 1.0*Y_[0]);
 
     for (label n=1; n<Y_.size(); n++)
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.C
index 4303a0e7fb..6c8cd1f50a 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.C
@@ -35,14 +35,15 @@ Foam::reactingMixture<ThermoType>::reactingMixture
     const fvMesh& mesh
 )
 :
+    speciesTable(),
     autoPtr<chemistryReader<ThermoType> >
     (
-        chemistryReader<ThermoType>::New(thermoDict)
+        chemistryReader<ThermoType>::New(thermoDict, *this)
     ),
     multiComponentMixture<ThermoType>
     (
         thermoDict,
-        autoPtr<chemistryReader<ThermoType> >::operator()().species(),
+        *this,
         autoPtr<chemistryReader<ThermoType> >::operator()().speciesThermo(),
         mesh
     ),
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
index 5c004eefeb..b29c0758ad 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
@@ -35,6 +35,7 @@ SourceFiles
 #ifndef reactingMixture_H
 #define reactingMixture_H
 
+#include "speciesTable.H"
 #include "chemistryReader.H"
 #include "multiComponentMixture.H"
 
@@ -50,6 +51,7 @@ namespace Foam
 template<class ThermoType>
 class reactingMixture
 :
+    public speciesTable,
     public autoPtr<chemistryReader<ThermoType> >,
     public multiComponentMixture<ThermoType>,
     public PtrList<Reaction<ThermoType> >
@@ -84,6 +86,16 @@ public:
 
         //- Read dictionary
         void read(const dictionary&);
+
+        label size() const
+        {
+            return PtrList<Reaction<ThermoType> >::size();
+        }
+
+        Reaction<ThermoType>& operator [] (const label i)
+        {
+            return PtrList<Reaction<ThermoType> >::operator[](i);
+        }
 };
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
index c4a66e7159..6c90e7c082 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
@@ -68,8 +68,8 @@ NonEquilibriumReversibleReaction
 )
 :
     Reaction<ReactionThermo>(species, thermoDatabase, dict),
-    fk_(species, dict),
-    rk_(species, dict)
+    fk_(species, dict.subDict("forward")),
+    rk_(species, dict.subDict("reverse"))
 {}
 
 
@@ -136,8 +136,20 @@ void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::write
 ) const
 {
     Reaction<ReactionThermo>::write(os);
+
+    os  << indent << "forward" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     fk_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "reverse" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     rk_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
index 7a477bcd4b..8a1c8e7e4d 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
@@ -67,10 +67,10 @@ FallOffReactionRate
     const dictionary& dict
 )
 :
-    k0_(species, dict),
-    kInf_(species, dict),
-    F_(dict),
-    thirdBodyEfficiencies_(species, dict)
+    k0_(species, dict.subDict("k0")),
+    kInf_(species, dict.subDict("kInf")),
+    F_(dict.subDict("F")),
+    thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
 {}
 
 
@@ -100,10 +100,33 @@ inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::write
     Ostream& os
 ) const
 {
+    os  << indent << "k0" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     k0_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "kInf" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     kInf_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "F" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     F_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "thirdBodyEfficiencies" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     thirdBodyEfficiencies_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
 }
 
 
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
index 9154126ada..bd8bdb4573 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
@@ -135,6 +135,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
index f031aa63f6..31bd392b5d 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
@@ -89,6 +89,16 @@ Foam::eConstThermo<EquationOfState>::New(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class EquationOfState>
 inline Foam::scalar Foam::eConstThermo<EquationOfState>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
index 014de5eef2..9a8352a064 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
@@ -133,6 +133,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
index 76ec41b620..0a4a35358c 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
@@ -89,6 +89,16 @@ Foam::hConstThermo<equationOfState>::New(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class equationOfState>
 inline Foam::scalar Foam::hConstThermo<equationOfState>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
index 7a4dcb3867..31954a4a56 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
@@ -151,6 +151,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index 620019e6d1..a164392c61 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -82,6 +82,16 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
 template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
index 9c414afaac..b552fad094 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
@@ -119,9 +119,6 @@ private:
         //- Check that input data is valid
         void checkInputData() const;
 
-        //- Check given temperature is within the range of the fitted coeffs
-        inline void checkT(const scalar T) const;
-
         //- Return the coefficients corresponding to the given temperature
         inline const coeffArray& coeffs(const scalar T) const;
 
@@ -153,6 +150,9 @@ public:
 
     // Member Functions
 
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
         // Fundamental properties
 
             //- Heat capacity at constant pressure [J/(kmol K)]
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index 05f0c0c0cc..de54035ea4 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -52,22 +52,6 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
 }
 
 
-template<class EquationOfState>
-inline void Foam::janafThermo<EquationOfState>::checkT(const scalar T) const
-{
-    if (T < Tlow_ || T > Thigh_)
-    {
-        FatalErrorIn
-        (
-            "janafThermo<EquationOfState>::checkT(const scalar T) const"
-        )   << "attempt to use janafThermo<EquationOfState>"
-               " out of temperature range "
-            << Tlow_ << " -> " << Thigh_ << ";  T = " << T
-            << abort(FatalError);
-    }
-}
-
-
 template<class EquationOfState>
 inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
 Foam::janafThermo<EquationOfState>::coeffs
@@ -75,8 +59,6 @@ Foam::janafThermo<EquationOfState>::coeffs
     const scalar T
 ) const
 {
-    checkT(T);
-
     if (T < Tcommon_)
     {
         return lowCpCoeffs_;
@@ -112,6 +94,31 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    if (T < Tlow_ || T > Thigh_)
+    {
+        WarningIn
+        (
+            "janafThermo<EquationOfState>::limit(const scalar T) const"
+        )   << "attempt to use janafThermo<EquationOfState>"
+               " out of temperature range "
+            << Tlow_ << " -> " << Thigh_ << ";  T = " << T
+            << endl;
+
+        return min(max(T, Tlow_), Thigh_);
+    }
+    else
+    {
+        return T;
+    }
+}
+
+
 template<class EquationOfState>
 inline Foam::scalar Foam::janafThermo<EquationOfState>::cp
 (
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
index c6f0f7a316..93268b9e6b 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
@@ -103,14 +103,15 @@ class specieThermo
 
     // Private Member Functions
 
-        //- return the temperature corresponding to the value of the
+        //- Return the temperature corresponding to the value of the
         //  thermodynamic property f, given the function f = F(T) and dF(T)/dT
         inline scalar T
         (
             scalar f,
             scalar T0,
             scalar (specieThermo::*F)(const scalar) const,
-            scalar (specieThermo::*dFdT)(const scalar) const
+            scalar (specieThermo::*dFdT)(const scalar) const,
+            scalar (specieThermo::*limit)(const scalar) const
         ) const;
 
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
index f8632376c2..b2d2edd5c1 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
@@ -43,7 +43,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
     scalar f,
     scalar T0,
     scalar (specieThermo<Thermo>::*F)(const scalar) const,
-    scalar (specieThermo<Thermo>::*dFdT)(const scalar) const
+    scalar (specieThermo<Thermo>::*dFdT)(const scalar) const,
+    scalar (specieThermo<Thermo>::*limit)(const scalar) const
 ) const
 {
     scalar Test = T0;
@@ -54,7 +55,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
     do
     {
         Test = Tnew;
-        Tnew = Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test);
+        Tnew =
+            (this->*limit)(Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test));
 
         if (iter++ > maxIter_)
         {
@@ -276,7 +278,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TH
     const scalar T0
 ) const
 {
-    return T(h, T0, &specieThermo<Thermo>::H, &specieThermo<Thermo>::Cp);
+    return T
+    (
+        h,
+        T0,
+        &specieThermo<Thermo>::H,
+        &specieThermo<Thermo>::Cp,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
@@ -287,7 +296,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::THs
     const scalar T0
 ) const
 {
-    return T(hs, T0, &specieThermo<Thermo>::Hs, &specieThermo<Thermo>::Cp);
+    return T
+    (
+        hs,
+        T0,
+        &specieThermo<Thermo>::Hs,
+        &specieThermo<Thermo>::Cp,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
@@ -298,7 +314,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TE
     const scalar T0
 ) const
 {
-    return T(e, T0, &specieThermo<Thermo>::E, &specieThermo<Thermo>::Cv);
+    return T
+    (
+        e,
+        T0,
+        &specieThermo<Thermo>::E,
+        &specieThermo<Thermo>::Cv,
+        &specieThermo<Thermo>::limit
+    );
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/Make/files b/src/turbulenceModels/incompressible/RAS/Make/files
index f18b124e06..e28a685903 100644
--- a/src/turbulenceModels/incompressible/RAS/Make/files
+++ b/src/turbulenceModels/incompressible/RAS/Make/files
@@ -16,6 +16,7 @@ LienCubicKELowRe/LienCubicKELowRe.C
 NonlinearKEShih/NonlinearKEShih.C
 LienLeschzinerLowRe/LienLeschzinerLowRe.C
 LamBremhorstKE/LamBremhorstKE.C
+kkLOmega/kkLOmega.C
 
 /* Wall functions */
 wallFunctions = derivedFvPatchFields/wallFunctions
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
new file mode 100644
index 0000000000..ad864bb3c7
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
@@ -0,0 +1,779 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kkLOmega.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "backwardsCompatibilityWallFunctions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(kkLOmega, 0);
+addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
+{
+    return(1.0 - exp(-sqrt(Ret)/Av_));
+}
+
+
+tmp<volScalarField> kkLOmega::fINT() const
+{
+    return
+    (
+        min
+        (
+            kl_/(Cint_*(kl_ + kt_ + kMin_)),
+            dimensionedScalar("1.0", dimless, 1.0)
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::fSS(const volScalarField& omega) const
+{
+    return(exp(-sqr(Css_*nu()*omega/(kt_ + kMin_))));
+}
+
+
+tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
+{
+    return(1.0/(A0_ + As_*(S/(omega_ + omegaMin_))));
+}
+
+
+tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& Rew) const
+{
+    return(scalar(1.0) - exp(-sqr(max(Rew - CtsCrit_, 0.0))/Ats_));
+}
+
+
+tmp<volScalarField> kkLOmega::fTaul
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& ktL
+) const
+{
+    return
+    (
+        scalar(1.0)
+      - exp
+        (
+            -CtauL_*ktL
+          /
+            (
+                sqr
+                (
+                    lambdaEff*omega_
+                  + dimensionedScalar
+                    (
+                        "ROTVSMALL",
+                        dimLength*inv(dimTime),
+                        ROOTVSMALL
+                    )
+                )
+            )
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::alphaT
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& fv,
+    const volScalarField& ktS
+) const
+{
+    return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
+}
+
+
+tmp<volScalarField> kkLOmega::fOmega
+(
+    const volScalarField& lambdaEff,
+    const volScalarField& lambdaT
+) const
+{
+    return
+    (
+        scalar(1.0)
+      - exp
+        (
+            -0.41
+           * pow4
+            (
+                lambdaEff
+              / (
+                    lambdaT
+                  + dimensionedScalar
+                    (
+                        "ROTVSMALL",
+                        lambdaT.dimensions(),
+                        ROOTVSMALL
+                    )
+                )
+            )
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
+{
+    return
+    (
+        max
+        (
+            kt_/nu()
+          /
+            (
+                omega
+              + dimensionedScalar("ROTVSMALL", omega.dimensions(), ROOTVSMALL)
+            )
+          -
+            CbpCrit_
+        ,
+            0.0
+        )
+    );
+}
+
+
+tmp<volScalarField> kkLOmega::gammaNAT
+(
+    const volScalarField& ReOmega,
+    const volScalarField& fNatCrit
+) const
+{
+    return
+    (
+        max
+        (
+            ReOmega
+          - CnatCrit_
+          / (
+                fNatCrit + dimensionedScalar("ROTVSMALL", dimless, ROOTVSMALL)
+            )
+          ,
+            0.0
+        )
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+kkLOmega::kkLOmega
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    transportModel& transport,
+    const word& turbulenceModelName,
+    const word& modelName
+)
+:
+    RASModel(modelName, U, phi, transport, turbulenceModelName),
+
+    A0_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "A0",
+            coeffDict_,
+            4.04
+        )
+    ),
+    As_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "As",
+            coeffDict_,
+            2.12
+        )
+    ),
+    Av_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Av",
+            coeffDict_,
+            6.75
+        )
+    ),
+    Abp_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Abp",
+            coeffDict_,
+            0.6
+        )
+    ),
+    Anat_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Anat",
+            coeffDict_,
+            200
+        )
+    ),
+    Ats_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ats",
+            coeffDict_,
+            200
+        )
+    ),
+    CbpCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CbpCrit",
+            coeffDict_,
+            1.2
+        )
+    ),
+    Cnc_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cnc",
+            coeffDict_,
+            0.1
+        )
+    ),
+    CnatCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CnatCrit",
+            coeffDict_,
+            1250
+        )
+    ),
+    Cint_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cint",
+            coeffDict_,
+            0.75
+        )
+    ),
+    CtsCrit_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CtsCrit",
+            coeffDict_,
+            1000
+        )
+    ),
+    CrNat_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CrNat",
+            coeffDict_,
+            0.02
+        )
+    ),
+    C11_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C11",
+            coeffDict_,
+            3.4e-6
+        )
+    ),
+    C12_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C12",
+            coeffDict_,
+            1.0e-10
+        )
+    ),
+    CR_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CR",
+            coeffDict_,
+            0.12
+        )
+    ),
+    CalphaTheta_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CalphaTheta",
+            coeffDict_,
+            0.035
+        )
+    ),
+    Css_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Css",
+            coeffDict_,
+            1.5
+        )
+    ),
+    CtauL_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CtauL",
+            coeffDict_,
+            4360
+        )
+    ),
+    Cw1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw1",
+            coeffDict_,
+            0.44
+        )
+    ),
+    Cw2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw2",
+            coeffDict_,
+            0.92
+        )
+    ),
+    Cw3_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw3",
+            coeffDict_,
+            0.3
+        )
+    ),
+    CwR_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CwR",
+            coeffDict_,
+            1.5
+        )
+    ),
+    Clambda_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Clambda",
+            coeffDict_,
+            2.495
+        )
+    ),
+    CmuStd_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CmuStd",
+            coeffDict_,
+            0.09
+        )
+    ),
+    Prtheta_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Prtheta",
+            coeffDict_,
+            0.85
+        )
+    ),
+    Sigmak_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Sigmak",
+            coeffDict_,
+            1
+        )
+    ),
+    Sigmaw_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Sigmaw",
+            coeffDict_,
+            1.17
+        )
+    ),
+    kt_
+    (
+        IOobject
+        (
+            "kt",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateK("kt", mesh_)
+    ),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateOmega("omega", mesh_)
+    ),
+    kl_
+    (
+        IOobject
+        (
+            "kl",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateK("kl", mesh_)
+    ),
+    nut_
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateNut("nut", mesh_)
+    ),
+    y_(mesh_)
+{
+    bound(kt_, kMin_);
+    bound(kl_, kMin_);
+    bound(omega_, omegaMin_);
+
+    nut_ = kt_/(omega_ + omegaMin_);
+    nut_.correctBoundaryConditions();
+
+    printCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volSymmTensorField> kkLOmega::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "R",
+                runTime_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            ((2.0/3.0)*I)*(kt_) - nut_*twoSymm(fvc::grad(U_)),
+            kt_.boundaryField().types()
+        )
+    );
+}
+
+
+tmp<volSymmTensorField> kkLOmega::devReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "devRhoReff",
+                runTime_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+           -nuEff()*dev(twoSymm(fvc::grad(U_)))
+        )
+    );
+}
+
+
+tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
+{
+    return
+    (
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+bool kkLOmega::read()
+{
+    if (RASModel::read())
+    {
+        A0_.readIfPresent(coeffDict());
+        As_.readIfPresent(coeffDict());
+        Av_.readIfPresent(coeffDict());
+        Abp_.readIfPresent(coeffDict());
+        Anat_.readIfPresent(coeffDict());
+        Abp_.readIfPresent(coeffDict());
+        Ats_.readIfPresent(coeffDict());
+        CbpCrit_.readIfPresent(coeffDict());
+        Cnc_.readIfPresent(coeffDict());
+        CnatCrit_.readIfPresent(coeffDict());
+        Cint_.readIfPresent(coeffDict());
+        CtsCrit_.readIfPresent(coeffDict());
+        CrNat_.readIfPresent(coeffDict());
+        C11_.readIfPresent(coeffDict());
+        C12_.readIfPresent(coeffDict());
+        CR_.readIfPresent(coeffDict());
+        CalphaTheta_.readIfPresent(coeffDict());
+        Css_.readIfPresent(coeffDict());
+        CtauL_.readIfPresent(coeffDict());
+        Cw1_.readIfPresent(coeffDict());
+        Cw2_.readIfPresent(coeffDict());
+        Cw3_.readIfPresent(coeffDict());
+        CwR_.readIfPresent(coeffDict());
+        Clambda_.readIfPresent(coeffDict());
+        CmuStd_.readIfPresent(coeffDict());
+        Prtheta_.readIfPresent(coeffDict());
+        Sigmak_.readIfPresent(coeffDict());
+        Sigmaw_.readIfPresent(coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void kkLOmega::correct()
+{
+    RASModel::correct();
+
+    if (!turbulence_)
+    {
+        return;
+    }
+
+    if (mesh_.changing())
+    {
+        y_.correct();
+        y_.boundaryField() = max(y_.boundaryField(), VSMALL);
+    }
+
+
+    const volScalarField kT(kt_ + kl_);
+
+    const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));
+
+    const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
+
+    const volScalarField fw
+    (
+        lambdaEff/(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL))
+    );
+
+    const volTensorField gradU(fvc::grad(U_));
+
+    const volScalarField omega(sqrt(2.0)*mag(skew(gradU)));
+
+    const volScalarField S2(2.0*magSqr(symm(gradU)));
+
+    const volScalarField ktS(fSS(omega)*fw*kt_);
+
+    const volScalarField nuts
+    (
+        fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
+       *fINT()
+       *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
+    );
+    const volScalarField Pkt(nuts*S2);
+
+    const volScalarField ktL(kt_ - ktS);
+    const volScalarField ReOmega(sqr(y_)*omega/nu());
+    const volScalarField nutl
+    (
+        min
+        (
+            C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
+          * sqrt(ktL)*lambdaEff/nu()
+          + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
+        ,
+            0.5*(kl_ + ktL)/sqrt(S2)
+        )
+    );
+
+    const volScalarField Pkl(nutl*S2);
+
+    const volScalarField alphaTEff
+    (
+        alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
+    );
+
+    // By pass s0urce term divided by kl_
+
+    const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
+
+    const volScalarField Rbp
+    (
+        CR_*(1.0 - exp(-gammaBP(omega)()/Abp_))*omega_
+      / (fw + fwMin)
+    );
+
+    const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
+    // Natural source term divided by kl_
+    const volScalarField Rnat
+    (
+        CrNat_*(1.0 - exp(-gammaNAT(ReOmega, fNatCrit)/Anat_))*omega
+    );
+
+    const volScalarField Dt(nu()*magSqr(fvc::grad(sqrt(kt_))));
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> ktEqn
+    (
+        fvm::ddt(kt_)
+      + fvm::div(phi_, kt_)
+      - fvm::Sp(fvc::div(phi_), kt_)
+      - fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
+     ==
+        Pkt
+      + (Rbp + Rnat)*kl_
+      - Dt
+      - fvm::Sp(omega_, kt_)
+    );
+
+    ktEqn().relax();
+    ktEqn().boundaryManipulate(kt_.boundaryField());
+
+    solve(ktEqn);
+    bound(kt_, kMin_);
+
+
+    const volScalarField Dl(nu()*magSqr(fvc::grad(sqrt(kl_))));
+
+    // Laminar kinetic energy equation
+    tmp<fvScalarMatrix> klEqn
+    (
+        fvm::ddt(kl_)
+      + fvm::div(phi_, kl_)
+      - fvm::Sp(fvc::div(phi_), kl_)
+      - fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
+     ==
+        Pkl
+      - fvm::Sp(Rbp, kl_)
+      - fvm::Sp(Rnat, kl_)
+      - Dl
+    );
+
+    klEqn().relax();
+    klEqn().boundaryManipulate(kl_.boundaryField());
+
+    solve(klEqn);
+    bound(kl_, kMin_);
+
+
+    omega_.boundaryField().updateCoeffs();
+    // Turbulence specific dissipation rate equation
+    tmp<fvScalarMatrix> omegaEqn
+    (
+        fvm::ddt(omega_)
+      + fvm::div(phi_, omega_)
+      - fvm::Sp(fvc::div(phi_), omega_)
+      - fvm::laplacian
+        (
+            DomegaEff(alphaTEff),
+            omega_,
+            "laplacian(alphaTEff,omega)"
+        )
+     ==
+        Cw1_*Pkt*omega_/(kt_ + kMin_)
+      + fvm::SuSp
+        (
+            (CwR_/(fw + fwMin) - 1.0)*kl_*(Rbp + Rnat)/(kt_ + kMin_)
+          , omega_
+        )
+      - fvm::Sp(Cw2_*omega_, omega_)
+      + Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)/pow3(y_)
+    );
+
+
+    omegaEqn().relax();
+    omegaEqn().boundaryManipulate(omega_.boundaryField());
+
+    solve(omegaEqn);
+    bound(omega_, omegaMin_);
+
+    // Re-calculate viscosity
+    nut_ = nuts + nutl;
+    nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
new file mode 100644
index 0000000000..b7d451b315
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::kkLOmega
+
+Description
+    Low Reynolds-number k-kl-omega turbulence model for
+    incompressible flows.
+
+    Turbulence model described in:
+    \verbatim
+        D. Keith Walters, Davor Cokljat
+        "A Three-Equation Eddy-Viscosity Model for Reynold-Averaged
+        Navier-Stokes Simulations of Transitional Flow"
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        kkLOmegaCoeffs
+        {
+            A0             4.04
+            As             2.12
+            Av             6.75
+            Abp            0.6
+            Anat           200
+            Ats            200
+            CbpCrit        1.2
+            Cnc            0.1
+            CnatCrit       1250
+            Cint           0.75
+            CtsCrit        1000
+            CrNat          0.02
+            C11            3.4e-6
+            C12            1.0e-10
+            CR             0.12
+            CalphaTheta    0.035
+            Css            1.5
+            CtauL          4360
+            Cw1            0.44
+            Cw2            0.92
+            Cw3            0.3
+            CwR            1.5
+            Clambda        2.495
+            CmuStd         0.09
+            Prtheta        0.85
+            Sigmak         1
+            Sigmaw         1.17
+        }
+    \endverbatim
+
+SourceFiles
+    kkLOmega.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kkLOmega_H
+#define kkLOmega_H
+
+#include "RASModel.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kkLOmega Declaration
+\*---------------------------------------------------------------------------*/
+
+class kkLOmega
+:
+    public RASModel
+{
+    // Private memmber functions
+
+        tmp<volScalarField> fv(const volScalarField& Ret) const;
+
+        tmp<volScalarField> fINT() const;
+
+        tmp<volScalarField> fSS(const volScalarField& omega) const;
+
+        tmp<volScalarField> Cmu(const volScalarField& S) const;
+
+        tmp<volScalarField> BetaTS(const volScalarField& Rew) const;
+
+        tmp<volScalarField> fTaul
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& ktL
+        ) const;
+
+        tmp<volScalarField> alphaT
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& fv,
+            const volScalarField& ktS
+        ) const;
+
+        tmp<volScalarField> fOmega
+        (
+            const volScalarField& lambdaEff,
+            const volScalarField& lambdaT
+        ) const;
+
+        tmp<volScalarField> gammaBP(const volScalarField& omega) const;
+
+        tmp<volScalarField> gammaNAT
+        (
+            const volScalarField& ReOmega,
+            const volScalarField& fNatCrit
+        ) const;
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar A0_;
+            dimensionedScalar As_;
+            dimensionedScalar Av_;
+            dimensionedScalar Abp_;
+            dimensionedScalar Anat_;
+            dimensionedScalar Ats_;
+            dimensionedScalar CbpCrit_;
+            dimensionedScalar Cnc_;
+            dimensionedScalar CnatCrit_;
+            dimensionedScalar Cint_;
+            dimensionedScalar CtsCrit_;
+            dimensionedScalar CrNat_;
+            dimensionedScalar C11_;
+            dimensionedScalar C12_;
+            dimensionedScalar CR_;
+            dimensionedScalar CalphaTheta_;
+            dimensionedScalar Css_;
+            dimensionedScalar CtauL_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar CwR_;
+            dimensionedScalar Clambda_;
+            dimensionedScalar CmuStd_;
+            dimensionedScalar Prtheta_;
+            dimensionedScalar Sigmak_;
+            dimensionedScalar Sigmaw_;
+
+
+        // Fields
+
+            volScalarField kt_;
+            volScalarField omega_;
+            volScalarField kl_;
+            volScalarField nut_;
+            wallDist y_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kkLOmega");
+
+    // Constructors
+
+        //- Construct from components
+        kkLOmega
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            transportModel& transport,
+            const word& turbulenceModelName = turbulenceModel::typeName,
+            const word& modelName = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kkLOmega()
+    {}
+
+
+    // Member Functions
+
+        //- Return the turbulence viscosity
+        virtual tmp<volScalarField> nut() const
+        {
+            return nut_;
+        }
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff(const volScalarField& alphaT) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DkEff", alphaT/Sigmak_ + nu())
+            );
+        }
+
+        //- Return the effective diffusivity for omega
+        tmp<volScalarField> DomegaEff(const volScalarField& alphaT) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DomegaEff", alphaT/Sigmaw_ + nu())
+            );
+        }
+
+        //- Return the laminar kinetic energy
+        virtual tmp<volScalarField> kl() const
+        {
+            return kl_;
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return kt_;
+        }
+
+        //- Return the turbulence specific dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "epsilon",
+                        mesh_.time().timeName(),
+                        mesh_
+                    ),
+                    kt_*omega_,
+                    omega_.boundaryField().types()
+                )
+            );
+        }
+
+        //- Return the Reynolds stress tensor
+        virtual tmp<volSymmTensorField> R() const;
+
+        //- Return the effective stress tensor including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+
+        //- Read RASProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
index 9210b03c39..caad132d7b 100755
--- a/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
+++ b/tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
@@ -1,7 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # run from this directory
 
-
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
index 9ae921391c..c92d02a0f7 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
@@ -15,11 +15,23 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-rho             rho [ 1 -3 0 0 0 0 0 ] 7854;
+rho
+{
+    rho         rhoInf;
+    rhoInf      rhoInf [ 1 -3 0 0 0 0 0 ] 7854;
+}
 
-nu              nu [ 0 0 0 0 0 0 0 ] 0.3;
+nu
+{
+    nu          nuInf;
+    nuInf       nuInf [ 0 0 0 0 0 0 0 ] 0.3;
+}
 
-E               E [ 1 -1 -2 0 0 0 0 ] 2e+11;
+E
+{
+    E           EInf;
+    EInf        EInf [ 1 -1 -2 0 0 0 0 ] 2e+11;
+}
 
 planeStress     yes;
 
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
index 9d09fa7901..f68549dbba 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
@@ -15,11 +15,23 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-C               C [ 0 2 -2 -1 0 0 0 ] 434;
+C
+{
+    C           CInf;
+    CInf        CInf [ 0 2 -2 -1 0 0 0 ] 434;
+}
 
-k               k [ 1 1 -3 -1 0 0 0 ] 60.5;
+k
+{
+    k           kInf;
+    kInf        kInf [ 1 1 -3 -1 0 0 0 ] 60.5;
+}
 
-alpha           alpha [ 0 0 0 -1 0 0 0 ] 1.1e-05;
+alpha
+{
+    alpha       alphaInf;
+    alphaInf    alphaInf [ 0 0 0 -1 0 0 0 ] 1.1e-05;
+}
 
 thermalStress   no;
 
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes
index 0f667a4132..a7c9b4c8d2 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes
@@ -20,6 +20,11 @@ d2dt2Schemes
     default         steadyState;
 }
 
+ddtSchemes
+{
+    default         Euler;
+}
+
 gradSchemes
 {
     default         leastSquares;