diff --git a/applications/solvers/combustion/XiFoam/UEqn.H b/applications/solvers/combustion/XiFoam/UEqn.H
index 9697c6e1ed..1e626d75b8 100644
--- a/applications/solvers/combustion/XiFoam/UEqn.H
+++ b/applications/solvers/combustion/XiFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index b551d2dd8b..26c1cbd060 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.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
@@ -55,6 +55,7 @@ Description
 #include "laminarFlameSpeed.H"
 #include "ignition.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,6 +73,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -79,7 +82,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -87,7 +89,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "rhoEqn.H"
             #include "UEqn.H"
@@ -103,7 +105,7 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 4168eb0e34..bef433d388 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/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
         (
@@ -24,9 +24,12 @@ if (transonic)
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +44,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
         (
@@ -50,9 +53,12 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/dieselEngineFoam/UEqn.H b/applications/solvers/combustion/dieselEngineFoam/UEqn.H
index 47912e8dd5..c6caf5989c 100644
--- a/applications/solvers/combustion/dieselEngineFoam/UEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/UEqn.H
@@ -8,7 +8,7 @@
       + dieselSpray.momentumSource()
     );
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 92bd3ba1f6..467f7ef9d7 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.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 "volPointInterpolation.H"
 #include "thermoPhysicsTypes.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,13 +62,14 @@ int main(int argc, char *argv[])
     #include "setInitialDeltaT.H"
     #include "startSummary.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readEngineTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -106,16 +108,17 @@ int main(int argc, char *argv[])
 
         chemistrySh = kappa*chemistry.Sh()();
 
-        #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        #include "rhoEqn.H"
+
+        for (pimple.start(); pimple.loop(); pimple++)
         {
+            #include "UEqn.H"
             #include "YEqn.H"
             #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/combustion/dieselEngineFoam/pEqn.H b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
index 049db9d250..42228cc7b0 100644
--- a/applications/solvers/combustion/dieselEngineFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField A(UEqn.A());
 U = UEqn.H()/A;
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -12,7 +12,7 @@ if (transonic)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -23,9 +23,12 @@ if (transonic)
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -36,7 +39,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -47,9 +50,12 @@ else
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index c5ea609ed7..35d3a304b2 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.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
@@ -41,6 +41,7 @@ Description
 #include "OFstream.H"
 #include "Switch.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,13 +59,14 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -98,15 +100,15 @@ int main(int argc, char *argv[])
         chemistrySh = kappa*chemistry.Sh()();
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
+            #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/dieselFoam/pEqn.H b/applications/solvers/combustion/dieselFoam/pEqn.H
index cb9f95c03f..5c4bd5ddf6 100644
--- a/applications/solvers/combustion/dieselFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselFoam/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
         (
@@ -26,9 +26,12 @@ if (transonic)
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -43,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
         (
@@ -54,9 +57,12 @@ else
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
index 9697c6e1ed..1e626d75b8 100644
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 4168eb0e34..bef433d388 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/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
         (
@@ -24,9 +24,12 @@ if (transonic)
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +44,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
         (
@@ -50,9 +53,12 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 42fe69390a..14539b6fe6 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.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 "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +52,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -58,7 +61,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -68,14 +70,14 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index 197cb83ff9..f3131fbf63 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.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 "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,10 +55,9 @@ int main(int argc, char *argv[])
     {
         Info<< "\n Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
-        for (int ucorr=0; ucorr<nOuterCorr; ucorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
 
@@ -67,7 +69,7 @@ int main(int argc, char *argv[])
 
             hUEqn.relax();
 
-            if (momentumPredictor)
+            if (pimple.momentumPredictor())
             {
                 if (rotating)
                 {
@@ -87,7 +89,7 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 volScalarField rAU(1.0/hUEqn.A());
                 surfaceScalarField ghrAUf(magg*fvc::interpolate(h*rAU));
@@ -107,7 +109,7 @@ int main(int argc, char *argv[])
                     + fvc::ddtPhiCorr(rAU, h, hU, phi)
                     - phih0;
 
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
                 {
                     fvScalarMatrix hEqn
                     (
@@ -116,16 +118,15 @@ int main(int argc, char *argv[])
                       - fvm::laplacian(ghrAUf, h)
                     );
 
-                    if (ucorr < nOuterCorr-1 || corr < nCorr-1)
-                    {
-                        hEqn.solve();
-                    }
-                    else
-                    {
-                        hEqn.solve(mesh.solver(h.name() + "Final"));
-                    }
+                    hEqn.solve
+                    (
+                        mesh.solver
+                        (
+                            h.select(pimple.finalInnerIter(corr, nonOrth))
+                        )
+                    );
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (nonOrth == pimple.nNonOrthCorr())
                     {
                         phi += hEqn.flux();
                     }
diff --git a/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H b/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
index bf70e2babd..4619e78f31 100644
--- a/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
+++ b/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
@@ -1,5 +1,5 @@
-#   include "readPISOControls.H"
+    #include "readPISOControls.H"
 
-    int nAlphaCorr(readInt(piso.lookup("nAlphaCorr")));
+    int nAlphaCorr(readInt(pisoDict.lookup("nAlphaCorr")));
 
-    Switch correctAlpha(piso.lookup("correctAlpha"));
+    Switch correctAlpha(pisoDict.lookup("correctAlpha"));
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index 0c6129ebe7..d0b15b3c11 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.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
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        for (int outerCorr=0; outerCorr<pimple.nOuterCorr(); outerCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "rhoEqn.H"
             #include "gammaPsi.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
index 7645b2d4aa..9161c80563 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
@@ -1,13 +1,7 @@
 {
-    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")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
@@ -36,7 +30,7 @@
         #include "alphaEqns.H"
     }
 
-    if (oCorr == 0)
+    if (pimple.corr() == 0)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
index a973c23dfe..532d2bcc05 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
@@ -1,13 +1,7 @@
 {
-    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")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
@@ -38,7 +32,7 @@
         #include "alphaEqns.H"
     }
 
-    if (oCorr == 0)
+    if (pimple.corr() == 0)
     {
         interface.correct();
     }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 161d30f8a6..2adf446fda 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.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
@@ -54,6 +54,9 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -61,8 +64,6 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    pimpleControl pimple(mesh);
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
 
@@ -116,7 +117,7 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         // --- Outer-corrector loop
-        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "alphaEqnsSubCycle.H"
 
@@ -125,7 +126,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<oimple.nCorr(); corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
index f681ba6594..c3b3726d12 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
@@ -42,7 +42,7 @@
 
     adjustPhi(phi, U, pcorr);
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pcorrEqn
         (
@@ -51,7 +51,7 @@
 
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
index c93b1c0bfa..38bacbf4ea 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
@@ -4,7 +4,7 @@
 
     label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
-    if (nAlphaSubCycles > 1 && nOuterCorr != 1)
+    if (nAlphaSubCycles > 1 && pimple.nOuterCorr() != 1)
     {
         FatalErrorIn(args.executable())
             << "Sub-cycling alpha is only allowed for PISO, "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index fda51d8546..146293e72c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.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
@@ -51,14 +51,15 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    pimpleControl pimple(mesh);
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -74,7 +75,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Outer-corrector loop
-        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "alphaEqnsSubCycle.H"
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 3378fe4b5d..28c42abe79 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -4,7 +4,7 @@
 
     tmp<fvScalarMatrix> p_rghEqnComp;
 
-    if (transonic)
+    if (pimple.transonic())
     {
         p_rghEqnComp =
         (
diff --git a/applications/solvers/multiphase/compressibleInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
index f66da8e5f9..954fcb27f2 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
@@ -10,7 +10,7 @@
         readLabel(pimple.dict().lookup("nAlphaSubCycles"))
     );
 
-    if (nAlphaSubCycles > 1 && nOuterCorr != 1)
+    if (nAlphaSubCycles > 1 && pimple.nOuterCorr() != 1)
     {
         FatalErrorIn(args.executable())
             << "Sub-cycling alpha is only allowed for PISO, "
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index a297b5fcab..8104cdde82 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -70,7 +70,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
index 74ebde2b32..9f4e800092 100644
--- a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.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
@@ -174,7 +174,8 @@ Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
     (
         readScalar
         (
-            mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
+            mixture.U().mesh().solutionDict().subDict("PIMPLE").
+                lookup("cAlpha")
         )
     ),
     sigma12_(mixture.lookup("sigma12")),
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 3e78c97ab3..096141d75a 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.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
diff --git a/applications/solvers/multiphase/settlingFoam/createFields.H b/applications/solvers/multiphase/settlingFoam/createFields.H
index eaa21a8065..057487f7ac 100644
--- a/applications/solvers/multiphase/settlingFoam/createFields.H
+++ b/applications/solvers/multiphase/settlingFoam/createFields.H
@@ -366,7 +366,7 @@
     (
         p,
         p_rgh,
-        mesh.solutionDict().subDict("PISO"),
+        mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index c588ee44c3..ac45f62edf 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        for (int oCorr=0; oCorr<pimple.nOuterCorr(); oCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "calcVdj.H"
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index 95395d97d4..00581cd302 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("PIMPLE"), pRefCell, pRefValue);
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index a345c4e53b..7116c7450d 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,5 +1,6 @@
     #include "readTimeControls.H"
+    #include "readPISOControls.H"
 
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
+    int nAlphaCorr(readInt(pisoDict.lookup("nAlphaCorr")));
 
-    Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
+    Switch correctAlpha(pisoDict.lookup("correctAlpha"));