diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index 2f11a8db1a..3c4da27187 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -53,6 +53,18 @@ Description
     PDR (porosity/distributed resistance) modelling is included to handle
     regions containing blockages which cannot be resolved by the mesh.
 
+    The fields used by this solver are:
+
+    betav:  Volume porosity
+    Lobs:   Average diameter of obstacle in cell (m)
+    Aw:     Obstacle surface area per unit volume (1/m)
+    CR:     Drag tensor (1/m)
+    CT:     Turbulence generation parameter (1/m)
+    Nv:     Number of obstacles in cell per unit volume (m^-2)
+    nsv:    Tensor whose diagonal indicates the number to substract from
+            Nv to get the number of obstacles crossing the flow in each
+            direction.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
@@ -64,6 +76,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "bound.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +94,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -90,7 +105,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -98,26 +112,34 @@ int main(int argc, char *argv[])
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "bEqn.H"
-            #include "ftEqn.H"
-            #include "huEqn.H"
-            #include "hEqn.H"
+            #include "UEqn.H"
 
-            if (!ign.ignited())
+            // --- PISO loop
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
-                hu == h;
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "huEqn.H"
+                #include "hEqn.H"
+
+                if (!ign.ignited())
+                {
+                    hu == h;
+                }
+
+                #include "pEqn.H"
             }
 
-            #include "pEqn.H"
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "\nExecutionTime = "
diff --git a/applications/solvers/combustion/PDRFoam/UEqn.H b/applications/solvers/combustion/PDRFoam/UEqn.H
index deb00a3121..70912c7cfa 100644
--- a/applications/solvers/combustion/PDRFoam/UEqn.H
+++ b/applications/solvers/combustion/PDRFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     volSymmTensorField invA(inv(I*UEqn.A() + drag->Dcu()));
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 9488e60e72..3c6a6341f9 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = invA & 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*invA, 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*invA, 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/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index 26c1cbd060..fc3658269a 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -88,10 +88,11 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        #include "rhoEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "rhoEqn.H"
             #include "UEqn.H"
 
             #include "ftEqn.H"
@@ -110,7 +111,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         rho = thermo.rho();
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index 89cc360d9f..bb1f3f2018 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.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 "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,13 +52,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"
@@ -71,16 +73,23 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "hEqn.H"
-            #include "pEqn.H"
-        }
+            #include "UEqn.H"
 
-        turbulence->correct();
+            // --- PISO loop
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            {
+                #include "hEqn.H"
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 467f7ef9d7..94ab928999 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
@@ -79,10 +79,6 @@ int main(int argc, char *argv[])
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
         mesh.move();
-        const_cast<volPointInterpolation&>
-        (
-            volPointInterpolation::New(mesh)
-        ).updateMesh();
 
         dieselSpray.evolve();
 
@@ -122,9 +118,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "logSummary.H"
         #include "spraySummary.H"
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index 35d3a304b2..7fd7410247 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C
@@ -101,6 +101,7 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
+        // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "UEqn.H"
@@ -112,9 +113,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "spraySummary.H"
 
diff --git a/applications/solvers/combustion/engineFoam/UEqn.H b/applications/solvers/combustion/engineFoam/UEqn.H
index 64dd52b1b0..eff9e1a5d2 100644
--- a/applications/solvers/combustion/engineFoam/UEqn.H
+++ b/applications/solvers/combustion/engineFoam/UEqn.H
@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
     );
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 4fedd9b92c..b92da02651 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.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
@@ -58,6 +58,7 @@ Description
 #include "Switch.H"
 #include "OFstream.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,13 +76,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"
@@ -94,11 +96,11 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
+            #include "UEqn.H"
+
             #include "ftEqn.H"
             #include "bEqn.H"
             #include "huEqn.H"
@@ -109,10 +111,17 @@ int main(int argc, char *argv[])
                 hu == h;
             }
 
-            #include "pEqn.H"
-        }
+            // --- PISO loop
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "logSummary.H"
 
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 78a90c2357..1b75a93398 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/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
     (
@@ -12,18 +12,21 @@ 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
         (
             fvm::ddt(psi, p)
-          + fvm::div(phid, p, "div(phid,p)")
+          + fvm::div(phid, p)
           - 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();
         }
@@ -34,7 +37,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
         (
@@ -43,9 +46,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/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index 4bef0f80ea..a517eccde1 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -81,9 +81,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 14539b6fe6..ca156c3f7f 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -81,9 +81,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         if (runTime.write())
         {
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index 539fe9d938..9c22e63c51 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -83,9 +83,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 746fb1224f..a5e92c8479 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -83,7 +83,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 5195655a33..ad9e04cc7e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -88,7 +88,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
index c8ed4e2cbd..96c07a1756 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.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
@@ -34,10 +34,11 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nReading field U\n" << endl;
     volVectorField U
@@ -57,9 +58,9 @@ int main(int argc, char *argv[])
 
     moleculeCloud molecules(mesh, pot);
 
-#   include "temperatureAndPressureVariables.H"
+    #include "temperatureAndPressureVariables.H"
 
-#   include "readmdEquilibrationDict.H"
+    #include "readmdEquilibrationDict.H"
 
     label nAveragingSteps = 0;
 
@@ -67,18 +68,17 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-
         nAveragingSteps++;
 
         Info<< "Time = " << runTime.timeName() << endl;
 
         molecules.evolve();
 
-#       include "meanMomentumEnergyAndNMols.H"
+        #include "meanMomentumEnergyAndNMols.H"
 
-#       include "temperatureAndPressure.H"
+        #include "temperatureAndPressure.H"
 
-#       include "temperatureEquilibration.H"
+        #include "temperatureEquilibration.H"
 
         runTime.write();
 
@@ -96,3 +96,6 @@ int main(int argc, char *argv[])
 
     return 0;
 }
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
index 6bd90a8128..779172803d 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.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
@@ -34,10 +34,11 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nReading field U\n" << endl;
     volVectorField U
@@ -57,7 +58,7 @@ int main(int argc, char *argv[])
 
     moleculeCloud molecules(mesh, pot);
 
-#   include "temperatureAndPressureVariables.H"
+    #include "temperatureAndPressureVariables.H"
 
     label nAveragingSteps = 0;
 
@@ -65,16 +66,15 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-
         nAveragingSteps++;
 
         Info<< "Time = " << runTime.timeName() << endl;
 
         molecules.evolve();
 
-#       include "meanMomentumEnergyAndNMols.H"
+        #include "meanMomentumEnergyAndNMols.H"
 
-#       include "temperatureAndPressure.H"
+        #include "temperatureAndPressure.H"
 
         runTime.write();
 
@@ -92,3 +92,6 @@ int main(int argc, char *argv[])
 
     return 0;
 }
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index adb08c1ee5..ce177dfeb9 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -95,7 +95,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index ca7724f101..3c3eb76ac0 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -88,11 +88,14 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
-
-            rho = thermo.rho();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
index 6b5901eac8..3bedfc756d 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
@@ -101,7 +101,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 215947c1e8..8d7bcc372d 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -79,7 +79,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index f3131fbf63..a86a3977b3 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
 
         #include "CourantNo.H"
 
+        // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index ac0cfb7cf8..3d0eabcef8 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -103,9 +103,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index 35394576f8..e52431f1e8 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/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/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 27921d6396..63d746fe79 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -32,7 +32,7 @@
       + massSource.SuTot()
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -42,17 +42,10 @@
 
         pEqn.solve
         (
-            mesh.solver
-            (
-                p.select
-                (
-                    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/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index 35dfc7e616..638d83da0c 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.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
@@ -49,6 +49,7 @@ Description
 #include "porousZones.H"
 #include "timeActivatedExplicitSource.H"
 #include "SLGThermo.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,6 +71,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -77,7 +80,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "readAdditionalSolutionControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -90,17 +92,25 @@ int main(int argc, char *argv[])
 
         #include "chemistry.H"
         #include "rhoEqn.H"
-        #include "UEqn.H"
-        #include "YEqn.H"
-        #include "hsEqn.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"
+            #include "YEqn.H"
+            #include "hsEqn.H"
 
-        turbulence->correct();
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 6d0a023ef0..6ae930c7a7 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -98,9 +98,12 @@ int main(int argc, char *argv[])
                 {
                     #include "pEqn.H"
                 }
-            }
 
-            turbulence->correct();
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
 
             rho = thermo.rho();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 0525585c17..ae2a06a363 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -92,11 +92,14 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
-
-            rho = thermo.rho();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
+        rho = thermo.rho();
+
         if (runTime.write())
         {
             chemistry.dQ()().write();
diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
index 8c509ab72c..d028a1376f 100644
--- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
+++ b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.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 "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -58,26 +61,32 @@ int main(int argc, char *argv[])
         #include "readBubbleFoamControls.H"
         #include "CourantNo.H"
 
-        #include "alphaEqn.H"
-        #include "liftDragCoeffs.H"
-        #include "UEqns.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 "alphaEqn.H"
+            #include "liftDragCoeffs.H"
+            #include "UEqns.H"
 
-            if (correctAlpha)
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
-                #include "alphaEqn.H"
+                #include "pEqn.H"
+
+                if (correctAlpha)
+                {
+                    #include "alphaEqn.H"
+                }
+            }
+
+            #include "DDtU.H"
+
+            if (pimple.turbCorr())
+            {
+                #include "kEpsilon.H"
             }
         }
 
-        #include "DDtU.H"
-
-        #include "kEpsilon.H"
-
         #include "write.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/multiphase/bubbleFoam/createFields.H b/applications/solvers/multiphase/bubbleFoam/createFields.H
index efb41613a6..9b1d9aa2ae 100644
--- a/applications/solvers/multiphase/bubbleFoam/createFields.H
+++ b/applications/solvers/multiphase/bubbleFoam/createFields.H
@@ -192,4 +192,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/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H
index faae92d397..453575b19b 100644
--- a/applications/solvers/multiphase/bubbleFoam/pEqn.H
+++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H
@@ -42,7 +42,7 @@
         alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -50,9 +50,13 @@
         );
 
         pEqn.setReference(pRefCell, pRefValue);
-        pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
+
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
diff --git a/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H b/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
index 4619e78f31..248f608f5c 100644
--- a/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
+++ b/applications/solvers/multiphase/bubbleFoam/readBubbleFoamControls.H
@@ -1,5 +1,4 @@
-    #include "readPISOControls.H"
 
-    int nAlphaCorr(readInt(pisoDict.lookup("nAlphaCorr")));
+    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
 
-    Switch correctAlpha(pisoDict.lookup("correctAlpha"));
+    Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index d0b15b3c11..64eb38f71e 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -78,9 +78,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 146293e72c..d8204d7f2d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -88,12 +88,15 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         rho = alpha1*rho1 + alpha2*rho2;
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index 8eff1637e2..14fa1cfe7a 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -93,9 +93,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 5bda42694b..84d6552f02 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -118,9 +118,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 8104cdde82..d74283cea5 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -93,9 +93,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 5059e55973..2d1b874c15 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -94,9 +94,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
index 496e2772ea..498d180aa0 100644
--- a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
@@ -87,12 +87,15 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         #include "continuityErrs.H"
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index 85198665c6..e731c2579b 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -85,9 +85,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index f9d92768ec..533466ef1f 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -82,9 +82,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index ac45f62edf..483cb276cd 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -81,7 +81,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            #include "kEpsilon.H"
+            if (pimple.turbCorr())
+            {
+                #include "kEpsilon.H"
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 89ab5aebb3..39587a007b 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -81,7 +81,10 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();
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/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 82fc1280c1..53d11f9e06 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -58,7 +58,7 @@
         alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -66,9 +66,13 @@
         );
 
         pEqn.setReference(pRefCell, pRefValue);
-        pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
+
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readPPProperties.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readPPProperties.H
index c1b319cbc9..a655736c54 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readPPProperties.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readPPProperties.H
@@ -10,27 +10,12 @@
         )
     );
 
-    scalar preAlphaExp
-    (
-        readScalar(ppProperties.lookup("preAlphaExp"))
-    );
+    scalar preAlphaExp(readScalar(ppProperties.lookup("preAlphaExp")));
 
-    scalar alphaMax
-    (
-        readScalar(ppProperties.lookup("alphaMax"))
-    );
+    scalar alphaMax(readScalar(ppProperties.lookup("alphaMax")));
 
-    scalar expMax
-    (
-        readScalar(ppProperties.lookup("expMax"))
-    );
+    scalar expMax(readScalar(ppProperties.lookup("expMax")));
 
-    dimensionedScalar g0
-    (
-        ppProperties.lookup("g0")
-    );
+    dimensionedScalar g0(ppProperties.lookup("g0"));
 
-    Switch packingLimiter
-    (
-        ppProperties.lookup("packingLimiter")
-    );
+    Switch packingLimiter(ppProperties.lookup("packingLimiter"));
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index 7116c7450d..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"
 
-    int nAlphaCorr(readInt(pisoDict.lookup("nAlphaCorr")));
+    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
 
-    Switch correctAlpha(pisoDict.lookup("correctAlpha"));
+    Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 24f55327ee..17cf2e7ec0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.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
@@ -42,6 +42,8 @@ Description
 #include "phaseModel.H"
 #include "kineticTheoryModel.H"
 
+#include "pimpleControl.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -58,6 +60,8 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -71,27 +75,34 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "alphaEqn.H"
-
-        #include "liftDragCoeffs.H"
-
-        #include "UEqns.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 "alphaEqn.H"
 
-            if (correctAlpha && corr<nCorr-1)
+            #include "liftDragCoeffs.H"
+
+            #include "UEqns.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
-                #include "alphaEqn.H"
+                #include "pEqn.H"
+
+                if (correctAlpha && !pimple.finalIter())
+                {
+                    #include "alphaEqn.H"
+                }
+            }
+
+            #include "DDtU.H"
+
+            if (pimple.turbCorr())
+            {
+                #include "kEpsilon.H"
             }
         }
 
-        #include "DDtU.H"
-
-        #include "kEpsilon.H"
-
         #include "write.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index 367c09f716..cd5d46c824 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "pimpleControl.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -43,6 +44,8 @@ void Foam::pimpleControl::read()
     const dictionary& pimpleDict = dict();
     nOuterCorr_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
     nCorr_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    turbOnFinalIterOnly_ =
+        pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
 }
 
 
@@ -118,7 +121,8 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh)
     solutionControl(mesh, "PIMPLE"),
     nOuterCorr_(0),
     nCorr_(0),
-    corr_(0)
+    corr_(0),
+    turbOnFinalIterOnly_(true)
 {
     read();
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index a7007bd0ff..fe8be71ec4 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -63,6 +63,9 @@ protected:
             //- Current PIMPLE corrector
             label corr_;
 
+            //- Flag to indicate whether to only solve turbulence on final iter
+            bool turbOnFinalIterOnly_;
+
 
     // Protected Member Functions
 
@@ -130,6 +133,9 @@ public:
                 const label nonOrth
             ) const;
 
+            //- Helper function to identify whether to solve for turbulence
+            inline bool turbCorr() const;
+
 
     // Member Operators
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index 4c761dfe23..62c30f8452 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -110,6 +110,12 @@ inline bool Foam::pimpleControl::finalInnerIter
 }
 
 
+inline bool Foam::pimpleControl::turbCorr() const
+{
+    return !turbOnFinalIterOnly_ || finalIter();
+}
+
+
 inline void Foam::pimpleControl::operator++(int)
 {
     if (finalIter())
diff --git a/src/thermophysicalModels/properties/liquidProperties/C3H6O/C3H6O.C b/src/thermophysicalModels/properties/liquidProperties/C3H6O/C3H6O.C
index c4b8254bb5..e499fb845f 100644
--- a/src/thermophysicalModels/properties/liquidProperties/C3H6O/C3H6O.C
+++ b/src/thermophysicalModels/properties/liquidProperties/C3H6O/C3H6O.C
@@ -55,7 +55,7 @@ Foam::C3H6O::C3H6O()
         1.9774e+4
     ),
     rho_(71.426784, 0.2576, 508.2, 0.29903),
-    pv_(70.72, -5.685, -7.351, 6.3e-06, 2.0),
+    pv_(70.72, -5685, -7.351, 6.3e-06, 2.0),
     hl_(508.20, 846590.909090909, 1.036, -1.294, 0.672, 0.0),
     Cp_
     (
diff --git a/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSolution b/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSolution
index 8a927992fe..b981b69820 100755
--- a/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSolution
+++ b/tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSolution
@@ -22,6 +22,13 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-05;
+        relTol          0.1;
+    };
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-05;
         relTol          0;
     };
 
@@ -30,10 +37,25 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-6;
+        relTol          0.1;
+    };
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-6;
         relTol          0;
     };
 
-    "(U|ft|fu|b|Xi|Su|h|hu|R|k|epsilon)"
+    "(b|Xi|ft|h|hu|k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(b|Xi|ft|h|hu|k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -42,7 +64,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
index a99a6a1eb9..e4df9498a8 100644
--- a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
+++ b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
@@ -32,15 +32,7 @@ solvers
         relTol          0;
     }
 
-    "(Yi|k|epsilon)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
-    }
-
-   "(U|hs)"
+    "(k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -48,12 +40,32 @@ solvers
         relTol          0.1;
     }
 
-   "(U|hs)Final"
+    "(k|epsilon)Final"
+    {
+        $k;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "(U|hs)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|hs)Final"
     {
         $U;
         tolerance       1e-06;
         relTol          0;
     }
+
+    Yi
+    {
+        $hsFinal;
+    }
 }
 
 PIMPLE
diff --git a/tutorials/combustion/engineFoam/kivaTest/system/fvSolution b/tutorials/combustion/engineFoam/kivaTest/system/fvSolution
index 7f7303e2fc..bf662ae7b7 100644
--- a/tutorials/combustion/engineFoam/kivaTest/system/fvSolution
+++ b/tutorials/combustion/engineFoam/kivaTest/system/fvSolution
@@ -22,26 +22,47 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-06;
         relTol          0;
     }
 
     rho
     {
-        $p;
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
         tolerance       1e-05;
         relTol          0;
     }
 
-    "(U|ft|fu|b|Xi|Su|h|hu|R|k|epsilon)"
+    "(U|ft|h|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-05;
+        relTol          0.1;
+    }
+    "(U|ft|h|k|epsilon)Final"
+    {
+        $U;
+        tolerance       1e-05;
         relTol          0;
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 1;
diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
index 54fb2634d0..f2c453ff5c 100755
--- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
+++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
@@ -47,7 +47,15 @@ solvers
         relTol           0.0;
     }
 
-    "(U|hs)"
+    "(U|hs||k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|hs||k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -55,20 +63,9 @@ solvers
         relTol          0;
     }
 
-    "(U|hs)Final"
+    Yi
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
-    }
-
-    "(Yi|k|epsilon)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
+        $hsFinal;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
index 32e542337a..8eacc7bc9c 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
@@ -17,11 +17,18 @@ FoamFile
 
 solvers
 {
-    "(rho|G)"
+    rho
     {
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-05;
         relTol          0;
     }
 
@@ -29,7 +36,14 @@ solvers
     {
         solver          PBiCG;
         preconditioner  DILU;
-        tolerance       1e-05;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|hs|k|epsilon)Final"
+    {
+        $U;
+        tolerance       1e-06;
         relTol          0;
     }
 
@@ -56,14 +70,11 @@ solvers
 
     "(Yi|O2|N2|H2O)"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
+        $hsFinal;
     }
 }
 
-PISO
+PIMPLE
 {
     transonic       no;
     nOuterCorrectors 1;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
index af936fdba1..b913d158bf 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
@@ -17,55 +17,64 @@ FoamFile
 
 solvers
 {
-    p
-    {
-        solver           GAMG;
-
-        tolerance        0;
-        relTol           0.01;
-
-        smoother         DICGaussSeidel;
-        nPreSweeps       0;
-        nPostSweeps      2;
-
-        cacheAgglomeration true;
-
-        nCellsInCoarsestLevel 10;
-        agglomerator     faceAreaPair;
-        mergeLevels      1;
-
-        maxIter          20;
-    };
-    pFinal
-    {
-        $p
-        tolerance        1e-6;
-        relTol           0;
-    };
-    "(rho|G)"
+    rho
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-06;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-05;
         relTol          0;
-    };
-    "(Yi|hs)"
+    }
+
+    "(U|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
-    };
-    "(U|k|omega)"
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(U|hs|k|epsilon)Final"
     {
-        solver          smoothSolver;
+        $U;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    p
+    {
+        solver          GAMG;
+        tolerance       0;
+        relTol          0.1;
         smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    pFinal
+    {
+        $p;
         tolerance       1e-06;
         relTol          0;
-    };
+    }
+
+    "(Yi|O2|N2|H2O)"
+    {
+        $hsFinal;
+    }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
index a235f4ecbe..c008af3c1e 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
@@ -17,57 +17,64 @@ FoamFile
 
 solvers
 {
-    p
-    {
-        solver           GAMG;
-        tolerance        0;
-        relTol           0.01;
-
-        smoother         DICGaussSeidel;
-        nPreSweeps       0;
-        nPostSweeps      2;
-
-        cacheAgglomeration true;
-
-        nCellsInCoarsestLevel 10;
-        agglomerator     faceAreaPair;
-        mergeLevels      1;
-
-        maxIter          50;
-    };
-    pFinal
-    {
-        $p
-        tolerance        1e-6;
-        relTol           0;
-    };
-    "(rho|G)"
+    rho
     {
         solver          PCG;
         preconditioner  DIC;
-        maxIter         200;
-        tolerance       1e-06;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-05;
         relTol          0;
-    };
-    "(Yi|hs)"
+    }
+
+    "(U|hs|k|omega)"
     {
         solver          PBiCG;
         preconditioner  DILU;
-        maxIter         200;
         tolerance       1e-06;
-        relTol          0;
-    };
-    "(U|k|omega)"
+        relTol          0.1;
+    }
+
+    "(U|hs|k|omega)Final"
     {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        maxIter         200;
+        $U;
         tolerance       1e-06;
         relTol          0;
-    };
+    }
+
+    p
+    {
+        solver          GAMG;
+        tolerance       0;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "(Yi|O2|N2|H2O)"
+    {
+        $hsFinal;
+    }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     3;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSolution b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSolution
index 72274eb02e..9f903209b0 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/fvSolution
@@ -32,7 +32,7 @@ solvers
         relTol          0;
     }
 
-    "(U|hs)"
+    "(U|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -40,18 +40,13 @@ solvers
         relTol          0.1;
     }
 
-    "(U|hs)Final"
+    "(U|hs|k|epsilon)Final"
     {
         $U;
         relTol          0;
         tolerance       1e-05;
     }
 
-    "(k|epsilon)"
-    {
-        $UFinal;
-    }
-
     p_rgh
     {
         solver          GAMG;
diff --git a/tutorials/multiphase/bubbleFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/bubbleFoam/bubbleColumn/system/fvSolution
index 22f5e23dec..c22cf4cff5 100644
--- a/tutorials/multiphase/bubbleFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/bubbleFoam/bubbleColumn/system/fvSolution
@@ -22,27 +22,33 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-10;
+        relTol          0.1;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-10;
         relTol          0;
     }
 
-    "(Ua|Ub|k|epsilon)"
+    alpha
     {
         solver          PBiCG;
         preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
+        tolerance       1e-10;
+        relTol          0.1;
     }
 
-    "(alpha|beta)"
+    alphaFinal
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        $alpha;
         tolerance       1e-10;
         relTol          0;
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
index 78eb84c984..e2298856f0 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSolution
@@ -55,15 +55,7 @@ solvers
         relTol          0;
     }
 
-    "(k|omega)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0;
-    }
-
-    U
+    "(U|k|omega)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -71,7 +63,7 @@ solvers
         relTol          0.1;
     }
 
-    UFinal
+    "(U|k|omega)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
index ec33c03abb..f326f43350 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSolution
@@ -53,7 +53,7 @@ solvers
         relTol          0;
     }
 
-    U
+    "(U|k|omega)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -61,20 +61,12 @@ solvers
         relTol          0.1;
     }
 
-    UFinal
+    "(U|k|omega)Final"
     {
         $U;
         tolerance       1e-08;
         relTol          0;
     }
-
-    "(k|omega)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0;
-    }
 }
 
 PIMPLE
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
index 78eb84c984..e2298856f0 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSolution
@@ -55,15 +55,7 @@ solvers
         relTol          0;
     }
 
-    "(k|omega)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-08;
-        relTol          0;
-    }
-
-    U
+    "(U|k|omega)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -71,7 +63,7 @@ solvers
         relTol          0.1;
     }
 
-    UFinal
+    "(U|k|omega)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
index 36882c9e4f..359560da4e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
@@ -91,7 +91,16 @@ solvers
         maxIter         20;
     }
 
-    U
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-6;
+        relTol          0.1;
+        nSweeps         1;
+    }
+
+    "(U|k|epsilon)Final"
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
@@ -99,14 +108,6 @@ solvers
         relTol          0;
         nSweeps         1;
     }
-
-    "(k|epsilon|omega|R|nuTilda)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-8;
-        relTol          0;
-    }
 }
 
 PIMPLE
diff --git a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSolution b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSolution
index f088132683..0271dcc369 100644
--- a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSolution
+++ b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSolution
@@ -77,16 +77,16 @@ solvers
         maxIter         20;
     }
 
-    U
+    "(U|k)"
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
         tolerance       1e-06;
-        relTol          0;
+        relTol          0.1;
         nSweeps         1;
     }
 
-    "(k|B|nuTilda)"
+    "(U|k)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSolution b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSolution
index daffd13542..3b4404dafd 100644
--- a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSolution
@@ -41,7 +41,7 @@ solvers
         relTol          0;
     }
 
-    U
+    "(U|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,7 +49,7 @@ solvers
         relTol          0;
     }
 
-    "(k|epsilon|R|nuTilda)"
+    "(U|k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
index 50d3809fbe..3fc3db5f29 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution
@@ -21,7 +21,7 @@ solvers
     {
         solver          GAMG;
         tolerance       1e-08;
-        relTol          0;
+        relTol          0.1;
         smoother        DIC;
         nPreSweeps      0;
         nPostSweeps     2;
@@ -32,7 +32,22 @@ solvers
         mergeLevels     1;
     }
 
-    "(Ua|Ub|Theta|k|epsilon)"
+    pFinal
+    {
+        $p;
+        tolerance       1e-08;
+        relTol          0;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -40,7 +55,15 @@ solvers
         relTol          0;
     }
 
-    "(alpha|beta)"
+    alpha
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-10;
+        relTol          0.1;
+    }
+
+    alphaFinal
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,7 +72,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
index ac04ffc4e3..3fc3db5f29 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/system/fvSolution
@@ -21,7 +21,7 @@ solvers
     {
         solver          GAMG;
         tolerance       1e-08;
-        relTol          0;
+        relTol          0.1;
         smoother        DIC;
         nPreSweeps      0;
         nPostSweeps     2;
@@ -32,7 +32,22 @@ solvers
         mergeLevels     1;
     }
 
-    "(Ua|Ub|Theta|k|epsilon)"
+    pFinal
+    {
+        $p;
+        tolerance       1e-08;
+        relTol          0;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -40,7 +55,15 @@ solvers
         relTol          0;
     }
 
-    "(alpha|beta)"
+    alpha
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-10;
+        relTol          0.1;
+    }
+
+    alphaFinal
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,7 +72,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
@@ -59,15 +82,5 @@ PISO
     pRefValue       0;
 }
 
-relaxationFactors
-{
-    Ua              1;
-    Ub              1;
-    alpha           1;
-    beta            1;
-    Theta           1;
-    "(k|epsilon)"   1;
-}
-
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
index 87e79cf796..3fc3db5f29 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/system/fvSolution
@@ -21,7 +21,7 @@ solvers
     {
         solver          GAMG;
         tolerance       1e-08;
-        relTol          0;
+        relTol          0.1;
         smoother        DIC;
         nPreSweeps      0;
         nPostSweeps     2;
@@ -32,7 +32,22 @@ solvers
         mergeLevels     1;
     }
 
-    "(Ua|Ub|Theta|k|epsilon)"
+    pFinal
+    {
+        $p;
+        tolerance       1e-08;
+        relTol          0;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(k|epsilon)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -40,7 +55,15 @@ solvers
         relTol          0;
     }
 
-    "(alpha|beta)"
+    alpha
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-10;
+        relTol          0.1;
+    }
+
+    alphaFinal
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,12 +72,12 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     nAlphaCorr      2;
-    correctAlpha    no;
+    correctAlpha    yes;
     pRefCell        0;
     pRefValue       0;
 }