diff --git a/applications/solvers/basic/laplacianFoam/laplacianFoam.C b/applications/solvers/basic/laplacianFoam/laplacianFoam.C
index b498db2d97..e2efee8f45 100644
--- a/applications/solvers/basic/laplacianFoam/laplacianFoam.C
+++ b/applications/solvers/basic/laplacianFoam/laplacianFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,17 +42,17 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createFields.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nCalculating temperature distribution\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
         {
             solve
             (
diff --git a/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C b/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
index ce5580d50b..3c0c0f3b14 100644
--- a/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
+++ b/applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,32 +30,30 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
 
-#   include "setRootCase.H"
+    simpleControl simple(mesh);
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nCalculating scalar transport\n" << endl;
 
-#   include "CourantNo.H"
+    #include "CourantNo.H"
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
         {
             solve
             (
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index d2ef427e48..4bef0f80ea 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -59,11 +59,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -79,7 +77,7 @@ int main(int argc, char *argv[])
             #include "fuhsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 64ebb1a970..6813950e19 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<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     surfaceScalarField rhorAUf(fvc::interpolate(rho*rAU));
 
@@ -28,18 +28,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     p_rghEqn.solve
     (
-        mesh.solver
-        (
-            p_rgh.select
-            (
-                pimple.finalIter()
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-        )
+        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi += p_rghEqn.flux();
     }
diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
index 9697c6e1ed..1e626d75b8 100644
--- a/applications/solvers/combustion/rhoReactingFoam/UEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index 7c227676eb..724f45e189 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -8,7 +8,7 @@
     volScalarField rAU(1.0/UEqn.A());
     U = rAU*UEqn.H();
 
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phiv
         (
@@ -30,7 +30,7 @@
           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
         );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
         {
             fvScalarMatrix pEqn
             (
@@ -40,18 +40,10 @@
 
             pEqn.solve
             (
-                mesh.solver
-                (
-                    p.select
-                    (
-                        pimple.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
+                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
             );
 
-            if (nonOrth == nNonOrthCorr)
+            if (nonOrth == pimple.nNonOrthCorr())
             {
                 phi += pEqn.flux();
             }
@@ -72,7 +64,7 @@
           + fvc::div(phi)
         );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
         {
             fvScalarMatrix pEqn
             (
@@ -82,18 +74,10 @@
 
             pEqn.solve
             (
-                mesh.solver
-                (
-                    p.select
-                    (
-                        pimple.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
+                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
             );
 
-            if (nonOrth == nNonOrthCorr)
+            if (nonOrth == pimple.nNonOrthCorr())
             {
                 phi += pEqn.flux();
             }
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index 36329f1215..539fe9d938 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -62,10 +62,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -81,7 +79,7 @@ int main(int argc, char *argv[])
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 27b06029b8..03bc6e76e7 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -11,7 +11,7 @@ UEqn().relax();
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index b665c8fc62..eea0ea9bff 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -5,12 +5,12 @@ rho.relax();
 
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -22,7 +22,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -33,18 +33,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -59,7 +51,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -71,18 +63,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 6db691f5ed..746fb1224f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -56,10 +56,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -70,7 +68,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -80,7 +78,7 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index 1ee1ed665e..97cc2cc256 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -50,14 +50,14 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPIMPLEControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "setInitialrDeltaT.H"
     #include "createFields.H"
     #include "createZones.H"
     #include "initContinuityErrs.H"
 
-    pimpleControl pimple(mesh);
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -65,10 +65,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -81,7 +79,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -93,7 +91,7 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setInitialrDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setInitialrDeltaT.H
index 54652dd8ed..353abc83a3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setInitialrDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setInitialrDeltaT.H
@@ -1,6 +1,6 @@
 scalar maxDeltaT
 (
-    pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
 );
 
 volScalarField rDeltaT
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setrDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setrDeltaT.H
index 7cad4d2679..d56283d13f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setrDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setrDeltaT.H
@@ -1,4 +1,6 @@
 {
+    const dictionary& pimpleDict = pimple.dict();
+
     scalar maxCo
     (
         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.8)
@@ -29,7 +31,7 @@
        /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
     );
 
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phid
         (
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
index 4efb76bd45..f3b115d0da 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
@@ -15,7 +15,7 @@ pZones.addResistance(UEqn());
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index 32058ffa9e..5e54f970cd 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -6,12 +6,12 @@ rho.relax();
 volScalarField rAU(1.0/UEqn().A());
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -24,7 +24,7 @@ if (transonic)
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -35,18 +35,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -62,7 +54,7 @@ else
         );
     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -74,18 +66,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 0bce64d7a3..5195655a33 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -61,10 +61,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -75,7 +73,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -85,7 +83,7 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4574a096ef..54af64c003 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -9,7 +9,7 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
     surfaceScalarField phid
     (
@@ -17,7 +17,7 @@ if (transonic)
         fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -32,7 +32,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -54,7 +54,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index 44b391a7b5..4fc258870b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -11,7 +11,7 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
     surfaceScalarField phid
     (
@@ -20,7 +20,7 @@ if (transonic)
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -37,7 +37,7 @@ if (transonic)
 
         tpEqn().solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi == tpEqn().flux();
         }
@@ -50,7 +50,7 @@ else
 
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -67,7 +67,7 @@ else
 
         tpEqn().solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= tpEqn().flux();
         }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 83898240ff..e210e971ca 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -36,6 +36,7 @@ Description
 #include "RASModel.H"
 #include "MRFZones.H"
 #include "thermalPorousZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,16 +49,16 @@ int main(int argc, char *argv[])
     #include "createZones.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index a8ddcbeb6c..63b6586cb4 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,16 +45,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
index 36b77f0399..9eb783feca 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
@@ -12,9 +12,9 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         surfaceScalarField phid
         (
@@ -54,7 +54,7 @@ if (transonic)
             pEqn.solve();
         }
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi == phic + pEqn.flux();
         }
@@ -62,7 +62,7 @@ if (transonic)
 }
 else
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         phi = fvc::interpolate(rho*U) & mesh.Sf();
         closedVolume = adjustPhi(phi, U, p);
@@ -88,7 +88,7 @@ else
         }
 
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
@@ -119,7 +119,7 @@ rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 
-if (!transonic)
+if (!simple.transonic())
 {
     rho.relax();
 }
diff --git a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
index c804e4e45d..13b8ebddc4 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
@@ -35,6 +35,7 @@ Description
 #include "RASModel.H"
 #include "mixedFvPatchFields.H"
 #include "bound.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,20 +47,20 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
-        if (!transonic)
+        if (!simple.transonic())
         {
             rho.storePrevIter();
         }
diff --git a/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C b/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
index 151d17a398..1ec8b86022 100644
--- a/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
+++ b/applications/solvers/electromagnetics/magneticFoam/magneticFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ Description
 #include "OSspecific.H"
 #include "magnet.H"
 #include "electromagneticConstants.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +66,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-    #include "readSIMPLEControls.H"
+
+    simpleControl simple(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +75,7 @@ int main(int argc, char *argv[])
 
     runTime++;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         solve(fvm::laplacian(murf, psi) + fvc::div(murf*Mrf));
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
index df6f90ac02..07f46ec998 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index f5cd55daca..adb08c1ee5 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -75,15 +75,13 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p_rgh.storePrevIter();
             }
@@ -92,7 +90,7 @@ int main(int argc, char *argv[])
             #include "TEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 31c4f484e0..8ac13dc93e 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -21,18 +21,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index 9a53619af5..cbe464fc02 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    if (momentumPredictor)
+    if (simple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index 1a5086e130..06d0f2df38 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,16 +61,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 0e8ae4e793..5d643837ff 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -22,7 +22,7 @@
 
         p_rghEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
index c8b9f13180..7dff741e66 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index f2deb2e2d5..ca7724f101 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -62,10 +62,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -76,7 +74,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p_rgh.storePrevIter();
             }
@@ -85,7 +83,7 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 963b4df548..b47d65b842 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -35,18 +35,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi += p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index 5fbfb3f764..ca28910aaf 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    if (momentumPredictor)
+    if (simple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 5891efcbf8..d70b03a45b 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Description
 #include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,16 +46,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
         rho.storePrevIter();
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 6f39aff1c1..17bf590f29 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -24,7 +24,7 @@
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index 51aaaadae3..70fc290088 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "radiationModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,16 +49,16 @@ int main(int argc, char *argv[])
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
         rho.storePrevIter();
 
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index e25669e04a..1c7863c2ff 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 template<class Type>
 void zeroCells
@@ -74,16 +75,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         laminarTransport.lookup("lambda") >> lambda;
@@ -121,7 +122,7 @@ int main(int argc, char *argv[])
             adjustPhi(phi, U, p);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
             {
                 fvScalarMatrix pEqn
                 (
@@ -131,7 +132,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (nonOrth == simple.nNonOrthCorr())
                 {
                     phi -= pEqn.flux();
                 }
@@ -184,7 +185,7 @@ int main(int argc, char *argv[])
             adjustPhi(phia, Ua, pa);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
             {
                 fvScalarMatrix paEqn
                 (
@@ -194,7 +195,7 @@ int main(int argc, char *argv[])
                 paEqn.setReference(paRefCell, paRefValue);
                 paEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (nonOrth == simple.nNonOrthCorr())
                 {
                     phia -= paEqn.flux();
                 }
diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index b897eec4b5..08c06d3f46 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -11,7 +11,7 @@ UEqn().relax();
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 0a4e6b1aaa..da9b7581bf 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 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<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     // Pressure corrector
     fvScalarMatrix pEqn
@@ -23,18 +23,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve
     (
-        mesh.solver
-        (
-            p.select
-            (
-                pimple.finalIter()
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-        )
+        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi -= pEqn.flux();
     }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
index 652773fea7..e43a40b94b 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
@@ -11,7 +11,7 @@ UEqn().relax();
 
 rAU = 1.0/UEqn().A();
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
index f137c0027f..e112db4621 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pcorrEqn
         (
@@ -60,7 +60,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 81a46f227b..7e84ef50fe 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 (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
@@ -14,7 +14,7 @@ if (p.needReference())
     fvc::makeAbsolute(phi, U);
 }
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     fvScalarMatrix pEqn
     (
@@ -25,18 +25,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve
     (
-        mesh.solver
-        (
-            p.select
-            (
-                pimple.finalIter()
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-        )
+        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi -= pEqn.flux();
     }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
index 14c7c6fdc7..6b5901eac8 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
@@ -60,7 +60,6 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
         #include "CourantNo.H"
-        pimple.read();
 
         // Make the fluxes absolute
         fvc::makeAbsolute(phi, U);
@@ -89,7 +88,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -97,7 +96,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H
index 8835bb6b48..8f48f5d7d4 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H
@@ -1,5 +1,6 @@
     #include "readTimeControls.H"
-    #include "readPIMPLEControls.H"
+
+    const dictionary& pimpleDict = pimple.dict();
 
     const bool correctPhi =
         pimpleDict.lookupOrDefault("correctPhi", false);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 6e2683d69c..215947c1e8 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -56,10 +56,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -68,7 +66,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -76,7 +74,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/simpleFoam/createFields.H b/applications/solvers/incompressible/simpleFoam/createFields.H
index b957c72650..947da10c9d 100644
--- a/applications/solvers/incompressible/simpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/createFields.H
@@ -26,14 +26,13 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-
     singlePhaseTransportModel laminarTransport(U, phi);
 
     autoPtr<incompressible::RASModel> turbulence
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 0b1500213a..4bd84f59c4 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -7,7 +7,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -18,7 +18,7 @@
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
index d14adbb493..e4614b7063 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H
@@ -6,13 +6,7 @@
     if (pZones.size())
     {
         // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
-        {
-            nUCorr = readInt
-            (
-                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
-            );
-        }
+        nUCorr = simple.dict().lookupOrDefault<int>("nUCorrectors", 0);
 
         if (nUCorr > 0)
         {
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
index a97ff3983c..5908396d18 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<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
 {
     tmp<fvScalarMatrix> tpEqn;
 
@@ -35,7 +35,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         tpEqn().solve();
     }
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == simple.nNonOrthCorr())
     {
         phi -= tpEqn().flux();
     }
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
index 3507fded85..123fd811ea 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
@@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "porousZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,6 +43,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    simpleControl simple(mesh);
+
     #include "createFields.H"
     #include "createPorousZones.H"
     #include "initContinuityErrs.H"
@@ -50,12 +54,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
diff --git a/applications/solvers/incompressible/simpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
index fdb5fe0a9d..af2a913615 100644
--- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,6 +32,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,16 +44,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index 7b8d04206d..23c30d4561 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -74,11 +74,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
-        pimple.read();
 
         runTime++;
 
@@ -94,7 +92,7 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -106,7 +104,7 @@ int main(int argc, char *argv[])
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
index 35394576f8..e52431f1e8 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
@@ -12,7 +12,7 @@
 
     pZones.addResistance(UEqn);
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index 4accc5dad9..fcbfad3efb 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -42,18 +42,10 @@
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index d58b74ab9c..e7f242fa5f 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -11,7 +11,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index ab9c56dccd..ac0cfb7cf8 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -74,10 +74,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -101,7 +99,7 @@ int main(int argc, char *argv[])
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 2e9099624c..de0c04577a 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -28,18 +28,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -54,7 +46,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -67,18 +59,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
index 61d805e01d..a64e50a2d2 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 8178c40af5..bbfb5be051 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<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     fvScalarMatrix p_rghEqn
     (
@@ -30,18 +30,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     p_rghEqn.solve
     (
-        mesh.solver
-        (
-            p_rgh.select
-            (
-                pimple.finalIter()
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-        )
+        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi += p_rghEqn.flux();
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 5ab5dfb8d8..6d0a023ef0 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -69,11 +69,9 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setMultiRegionDeltaT.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -96,7 +94,7 @@ int main(int argc, char *argv[])
                 #include "hsEqn.H"
 
                 // --- PISO loop
-                for (int corr=1; corr<=nCorr; corr++)
+                for (int corr=1; corr<=pimple.nCorr(); corr++)
                 {
                     #include "pEqn.H"
                 }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
index 2e8f979be4..31e00d8f35 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
@@ -10,7 +10,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 0e7e36f0ef..3bf304cecc 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -28,18 +28,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -54,7 +46,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -67,18 +59,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index fd5599736a..0525585c17 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -67,10 +67,8 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -89,7 +87,7 @@ int main(int argc, char *argv[])
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 01911faaa3..9a5761b59f 100644
--- a/applications/solvers/multiphase/cavitatingFoam/UEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index 2e35e0bfc3..0c6129ebe7 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -36,6 +36,7 @@ Description
 #include "barotropicCompressibilityModel.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +53,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -65,13 +68,13 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
+        for (int outerCorr=0; outerCorr<pimple.nOuterCorr(); outerCorr++)
         {
             #include "rhoEqn.H"
             #include "gammaPsi.H"
             #include "UEqn.H"
 
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 6b6fc73eb9..c604539cde 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -1,5 +1,5 @@
 {
-    if (nOuterCorr == 1)
+    if (pimple.nOuterCorr() == 1)
     {
         p =
         (
@@ -26,7 +26,7 @@
 
     #include "resetPhivPatches.H"
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -37,16 +37,12 @@
           - fvm::laplacian(rAUf, p)
         );
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
-        {
-            pEqn.solve(mesh.solver(p.name() + "Final"));
-        }
-        else
-        {
-            pEqn.solve(mesh.solver(p.name()));
-        }
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phiv += (phiGradp + pEqn.flux())/rhof;
         }
@@ -82,9 +78,9 @@
     U = HbyA - rAU*fvc::grad(p);
 
     // Remove the swirl component of velocity for "wedge" cases
-    if (piso.found("removeSwirl"))
+    if (pimple.dict().found("removeSwirl"))
     {
-        label swirlCmpt(readLabel(piso.lookup("removeSwirl")));
+        label swirlCmpt(readLabel(pimple.dict().lookup("removeSwirl")));
 
         Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
         U.field().replace(swirlCmpt, 0.0);
diff --git a/applications/solvers/multiphase/cavitatingFoam/readControls.H b/applications/solvers/multiphase/cavitatingFoam/readControls.H
index f53e7b9eb1..0dd2b77780 100644
--- a/applications/solvers/multiphase/cavitatingFoam/readControls.H
+++ b/applications/solvers/multiphase/cavitatingFoam/readControls.H
@@ -5,5 +5,3 @@ scalar maxAcousticCo
     readScalar(runTime.controlDict().lookup("maxAcousticCo"))
 );
 
-
-#include "readPISOControls.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 320953ab6f..257f6d48b5 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index a6fe17ee9b..161d30f8a6 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -44,6 +44,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,6 +61,8 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
 
@@ -113,7 +116,7 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         // --- Outer-corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
         {
             #include "alphaEqnsSubCycle.H"
 
@@ -122,7 +125,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<oimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 0c1ad0192d..cdc541f0cd 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -4,7 +4,7 @@
 
     tmp<fvScalarMatrix> p_rghEqnComp;
 
-    if (transonic)
+    if (pimple.transonic())
     {
         p_rghEqnComp =
         (
@@ -39,7 +39,7 @@
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf();
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqnIncomp
         (
@@ -55,18 +55,10 @@
             )
            *p_rghEqnComp()
           + p_rghEqnIncomp,
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    oCorr == nOuterCorr-1
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             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 a2e4ef3747..c93b1c0bfa 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
@@ -1,15 +1,8 @@
-   #include "readPISOControls.H"
-   #include "readTimeControls.H"
+    #include "readTimeControls.H"
 
-    label nAlphaCorr
-    (
-        readLabel(piso.lookup("nAlphaCorr"))
-    );
+    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
 
-    label nAlphaSubCycles
-    (
-        readLabel(piso.lookup("nAlphaSubCycles"))
-    );
+    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
     if (nAlphaSubCycles > 1 && nOuterCorr != 1)
     {
@@ -19,14 +12,8 @@
             << exit(FatalError);
     }
 
-    bool correctPhi = true;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
+    bool correctPhi =
+        pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
 
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+    bool checkMeshCourantNo =
+        pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 4aed5e9c07..fda51d8546 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -41,6 +41,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +57,8 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -71,7 +74,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Outer-corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
         {
             #include "alphaEqnsSubCycle.H"
 
@@ -80,7 +83,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 047f2f3c4c..3378fe4b5d 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqnIncomp
         (
@@ -55,18 +55,10 @@
             )
            *p_rghEqnComp()
           + p_rghEqnIncomp,
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    oCorr == nOuterCorr-1
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             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 7e23354f47..f66da8e5f9 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
@@ -1,14 +1,13 @@
-   #include "readPISOControls.H"
    #include "readTimeControls.H"
 
     label nAlphaCorr
     (
-        readLabel(piso.lookup("nAlphaCorr"))
+        readLabel(pimple.dict().lookup("nAlphaCorr"))
     );
 
     label nAlphaSubCycles
     (
-        readLabel(piso.lookup("nAlphaSubCycles"))
+        readLabel(pimple.dict().lookup("nAlphaSubCycles"))
     );
 
     if (nAlphaSubCycles > 1 && nOuterCorr != 1)
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 7f28e63e94..74cffd7495 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,6 +44,7 @@ Description
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "fvcSmooth.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,9 +53,11 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialrDeltaT.H"
@@ -65,8 +68,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
-
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -76,13 +77,19 @@ int main(int argc, char *argv[])
         twoPhaseProperties.correct();
 
         #include "alphaEqnSubCycle.H"
-        turbulence->correct();
-        #include "UEqn.H"
 
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        turbulence->correct();
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
index 593b3503e1..b59b0ae5fe 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -1,6 +1,8 @@
-label nAlphaCorr(readLabel(piso.lookup("nAlphaCorr")));
+const dictionary& pimpleDict = pimple.dict();
 
-label nAlphaSubCycles(readLabel(piso.lookup("nAlphaSubCycles")));
+label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
+
+label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
index 210c65627d..c7ac175848 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
@@ -1,6 +1,6 @@
 scalar maxDeltaT
 (
-    piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
 );
 
 volScalarField rDeltaT
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
index 208d01b484..8525a8bf67 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@@ -1,52 +1,54 @@
 {
+    const dictionary& pimpleDict = pimple.dict();
+
     scalar maxCo
     (
-        piso.lookupOrDefault<scalar>("maxCo", 0.9)
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
     );
 
     scalar maxAlphaCo
     (
-        piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
     );
 
     scalar rDeltaTSmoothingCoeff
     (
-        piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
     );
 
     label nAlphaSpreadIter
     (
-        piso.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
     );
 
     scalar alphaSpreadDiff
     (
-        piso.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
     );
 
     scalar alphaSpreadMax
     (
-        piso.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
     );
 
     scalar alphaSpreadMin
     (
-        piso.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
     );
 
     label nAlphaSweepIter
     (
-        piso.lookupOrDefault<label>("nAlphaSweepIter", 5)
+        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
     );
 
     scalar rDeltaTDampingCoeff
     (
-        piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
     );
 
     scalar maxDeltaT
     (
-        piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
     );
 
     volScalarField rDeltaT0("rDeltaT0", rDeltaT);
@@ -130,7 +132,7 @@
 
     label nAlphaSubCycles
     (
-        readLabel(piso.lookup("nAlphaSubCycles"))
+        readLabel(pimpleDict.lookup("nAlphaSubCycles"))
     );
 
     rSubDeltaT = rDeltaT*nAlphaSubCycles;
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index 4c0d33eeb2..8eff1637e2 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ Description
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "MRFZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,11 +52,13 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createMRFZones.H"
     #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -66,7 +69,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
@@ -81,12 +83,16 @@ int main(int argc, char *argv[])
         #include "alphaEqnSubCycle.H"
         #include "zonePhaseVolumes.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H
index 80e9d09b3b..c26df02019 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H
@@ -17,7 +17,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index e4f8247bac..d281e177c8 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<=nNonOrthCorr; nonOrth++)
+    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -30,13 +30,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/UEqn.H b/applications/solvers/multiphase/interFoam/UEqn.H
index 320953ab6f..257f6d48b5 100644
--- a/applications/solvers/multiphase/interFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/UEqn.H
@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 073d058bfa..8dbb489c8d 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -1,12 +1,6 @@
-label nAlphaCorr
-(
-    readLabel(piso.lookup("nAlphaCorr"))
-);
+label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
 
-label nAlphaSubCycles
-(
-    readLabel(piso.lookup("nAlphaSubCycles"))
-);
+label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {
@@ -19,7 +13,7 @@ if (nAlphaSubCycles > 1)
         !(++alphaSubCycle).end();
     )
     {
-#       include "alphaEqn.H"
+        #include "alphaEqn.H"
         rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
     }
 
@@ -27,7 +21,7 @@ if (nAlphaSubCycles > 1)
 }
 else
 {
-#       include "alphaEqn.H"
+    #include "alphaEqn.H"
 }
 
 interface.correct();
diff --git a/applications/solvers/multiphase/interFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/correctPhi.H
index e08d3702e8..cb3503e7b8 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pcorrEqn
         (
@@ -44,7 +44,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 174cba3e1a..34dea334e8 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -128,7 +128,7 @@
     (
         p,
         p_rgh,
-        mesh.solutionDict().subDict("PISO"),
+        mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index d525913811..5bda42694b 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,10 +48,12 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -105,12 +108,16 @@ int main(int argc, char *argv[])
 
         #include "alphaEqnSubCycle.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 5bf3ca9418..84c1fe45b6 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -18,7 +18,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -29,13 +29,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H b/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
index 3640d73adc..d4e332ff38 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
@@ -1,14 +1,6 @@
-#   include "readTimeControls.H"
-#   include "readPISOControls.H"
+    #include "readTimeControls.H"
 
-    bool correctPhi = true;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
-
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+    bool correctPhi =
+        pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+    bool checkMeshCourantNo =
+        pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 98cd759b4f..a297b5fcab 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Description
     The momentum and other fluid properties are of the "mixture" and a single
     momentum equation is solved.
 
-    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
     For a two-fluid approach see twoPhaseEulerFoam.
 
@@ -44,6 +44,7 @@ Description
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "interpolationTable.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
@@ -60,6 +63,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -80,12 +84,16 @@ int main(int argc, char *argv[])
 
         #include "alphaEqnSubCycle.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index db1eccc5a5..d75a76766a 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -29,13 +29,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H b/applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H
index 72ae38ba06..155c0527b7 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H
@@ -22,7 +22,7 @@
 
     pZones.addResistance(UEqn);
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index f0949e290e..5059e55973 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Description
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "porousZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,11 +54,13 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createPorousZones.H"
     #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -68,7 +71,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
@@ -82,12 +84,16 @@ int main(int argc, char *argv[])
 
         #include "alphaEqnSubCycle.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
index 81e09f9b90..d816f8aced 100644
--- a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
@@ -1,12 +1,6 @@
-label nAlphaCorr
-(
-    readLabel(piso.lookup("nAlphaCorr"))
-);
-
-label nAlphaSubCycles
-(
-    readLabel(piso.lookup("nAlphaSubCycles"))
-);
+const dictionary& pimpleDict = pimple.dict();
+label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
+label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {
@@ -19,7 +13,7 @@ if (nAlphaSubCycles > 1)
         !(++alphaSubCycle).end();
     )
     {
-#       include "alphaEqns.H"
+        #include "alphaEqns.H"
         rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
     }
 
@@ -27,7 +21,7 @@ if (nAlphaSubCycles > 1)
 }
 else
 {
-#       include "alphaEqns.H"
+    #include "alphaEqns.H"
 }
 
 interface.correct();
diff --git a/applications/solvers/multiphase/interMixingFoam/createFields.H b/applications/solvers/multiphase/interMixingFoam/createFields.H
index 8ef77a5546..60543dc826 100644
--- a/applications/solvers/multiphase/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interMixingFoam/createFields.H
@@ -150,7 +150,7 @@
     (
         p,
         p_rgh,
-        mesh.solutionDict().subDict("PISO"),
+        mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
index 5da002283a..496e2772ea 100644
--- a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "threePhaseMixture.H"
 #include "threePhaseInterfaceProperties.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,12 +46,14 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,7 +62,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
@@ -74,12 +76,17 @@ int main(int argc, char *argv[])
         #include "alphaEqnsSubCycle.H"
 
         #define twoPhaseProperties threePhaseProperties
-        #include "UEqn.H"
 
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         #include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index 042f6f18a1..b8fe82ff14 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -17,7 +17,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
index eea578e6e5..b68160e62d 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
@@ -11,15 +11,11 @@ surfaceScalarField rhoPhi
 );
 
 {
-    label nAlphaCorr
-    (
-        readLabel(piso.lookup("nAlphaCorr"))
-    );
+    const dictionary& pimpleDict = pimple.dict();
 
-    label nAlphaSubCycles
-    (
-        readLabel(piso.lookup("nAlphaSubCycles"))
-    );
+    label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
+
+    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
@@ -44,7 +40,7 @@ surfaceScalarField rhoPhi
         #include "alphaEqn.H"
     }
 
-    if (nOuterCorr == 1)
+    if (pimple.nOuterCorr() == 1)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index 381e80b881..60f15ab553 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -98,7 +98,7 @@
     (
         p,
         p_rgh,
-        mesh.solutionDict().subDict("PISO"),
+        mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index eb03daea93..2e94017da2 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -56,27 +56,25 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "../interFoam/correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    pimpleControl pimple(mesh);
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -92,7 +90,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index c15d76312c..6cfa59cc7c 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -23,7 +23,7 @@
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -35,18 +35,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index 9c1e97eae3..85198665c6 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "multiphaseMixture.H"
 #include "turbulenceModel.H"
 #include "MRFZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,7 +45,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createMRFZones.H"
@@ -59,7 +62,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
@@ -73,12 +75,16 @@ int main(int argc, char *argv[])
         rho = mixture.rho();
         #include "zonePhaseVolumes.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H
index b47e711192..3b860b975a 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H
@@ -17,7 +17,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index 8be2e12bf6..9975ebac83 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -20,7 +20,7 @@
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -31,13 +31,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
index a2145bc0a2..43ee1e3eb9 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
index e3ebd10743..78fa0ec309 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@@ -76,7 +76,7 @@
     (
         p,
         p_rgh,
-        mesh.solutionDict().subDict("PISO"),
+        mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index feb3bf00a4..f9d92768ec 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "fvCFD.H"
 #include "multiphaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,10 +44,12 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -57,7 +60,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
@@ -70,12 +72,16 @@ int main(int argc, char *argv[])
         mixture.solve();
         rho = mixture.rho();
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
+         // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "pEqn.H"
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index eecec78206..3e78c97ab3 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -245,28 +245,15 @@ void Foam::multiphaseMixture::solve()
 
     const Time& runTime = mesh_.time();
 
-    label nAlphaSubCycles
-    (
-        readLabel
-        (
-            mesh_.solutionDict().subDict("PISO").lookup("nAlphaSubCycles")
-        )
-    );
+    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
 
-    label nAlphaCorr
-    (
-        readLabel(mesh_.solutionDict().subDict("PISO").lookup("nAlphaCorr"))
-    );
+    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
-    bool cycleAlpha
-    (
-        Switch(mesh_.solutionDict().subDict("PISO").lookup("cycleAlpha"))
-    );
+    label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
 
-    scalar cAlpha
-    (
-        readScalar(mesh_.solutionDict().subDict("PISO").lookup("cAlpha"))
-    );
+    bool cycleAlpha(Switch(pimpleDict.lookup("cycleAlpha")));
+
+    scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
 
 
     volScalarField& alpha = phases_.first();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index 712fcfc9a1..78bfc89bb4 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -30,13 +30,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H
index cd3985b364..fc750765b7 100644
--- a/applications/solvers/multiphase/settlingFoam/UEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/UEqn.H
@@ -14,7 +14,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index f359e95254..7b62411d2d 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<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     fvScalarMatrix p_rghEqn
     (
@@ -27,7 +27,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
     p_rghEqn.solve();
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi -= p_rghEqn.flux();
     }
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index c91229abd7..c588ee44c3 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "Switch.H"
 #include "plasticViscosity.H"
 #include "yieldStress.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,20 +51,21 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (pimple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
 
         #include "rhoEqn.H"
 
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
         {
             #include "calcVdj.H"
 
@@ -74,7 +76,7 @@ int main(int argc, char *argv[])
             #include "correctViscosity.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
index 5d9a317e08..3dbb374774 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index 42e85e1db9..01b6442e9d 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<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -25,18 +25,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
-            (
-                p_rgh.select
-                (
-                    pimple.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= p_rghEqn.flux();
         }
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index fb06b7fed3..89ab5aebb3 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -44,7 +44,6 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
-    #include "readPIMPLEControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
@@ -59,11 +58,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
-        pimple.read();
 
         runTime++;
 
@@ -79,7 +76,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index 00581cd302..95395d97d4 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -234,4 +234,4 @@
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index b956336389..a345c4e53b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,6 +1,5 @@
-#   include "readTimeControls.H"
-#   include "readPISOControls.H"
+    #include "readTimeControls.H"
 
-    int nAlphaCorr(readInt(piso.lookup("nAlphaCorr")));
+    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
 
-    Switch correctAlpha(piso.lookup("correctAlpha"));
+    Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index 6366738579..0e5b8c3776 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "MRFZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,16 +47,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
@@ -82,7 +83,7 @@ int main(int argc, char *argv[])
             adjustPhi(phi, U, p);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
             {
                 fvScalarMatrix pEqn
                 (
@@ -92,7 +93,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (nonOrth == simple.nNonOrthCorr())
                 {
                     phi -= pEqn.flux();
                 }