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/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;
 }