diff --git a/applications/solvers/combustion/XiFoam/UEqn.H b/applications/solvers/combustion/XiFoam/UEqn.H
index 9697c6e1e..1e626d75b 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 c652e2624..f4b93a51f 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -55,6 +55,7 @@ Description
 #include "laminarFlameSpeed.H"
 #include "ignition.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,6 +65,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
@@ -79,7 +83,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -87,7 +90,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++)
+        while (pimple.loop())
         {
 #           include "rhoEqn.H"
 #           include "UEqn.H"
@@ -103,7 +106,7 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            while (pimple.correct())
             {
 #              include "pEqn.H"
             }
diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H
index a40c9e9ef..a92694355 100644
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/applications/solvers/combustion/XiFoam/createFields.H
@@ -132,3 +132,5 @@
     fields.add(b);
     fields.add(h);
     fields.add(hu);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 280842ecc..970a11b8f 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,7 +26,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -52,7 +52,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index 747b035f1..1ee53dde6 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -35,6 +35,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createEngineMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 #   include "createTimeControls.H"
@@ -58,7 +62,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
@@ -75,7 +78,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct())
         {
 #           include "hEqn.H"
 #           include "pEqn.H"
diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H
index 7b1c086df..28f1f7d64 100644
--- a/applications/solvers/combustion/coldEngineFoam/createFields.H
+++ b/applications/solvers/combustion/coldEngineFoam/createFields.H
@@ -57,3 +57,5 @@
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/engineFoam/UEqn.H b/applications/solvers/combustion/engineFoam/UEqn.H
index 64dd52b1b..3016e1b9a 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 (piso.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 56b487bf8..0b55614b1 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@@ -57,6 +57,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "OFstream.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,6 +67,9 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createEngineMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -81,7 +85,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
@@ -97,7 +100,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct()) 
         {
 #           include "ftEqn.H"
 #           include "bEqn.H"
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 39b496731..ed86b9c97 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (piso.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++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -23,7 +23,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -34,7 +34,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -45,7 +45,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
index 617e2c137..85149d415 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSchemes
@@ -44,7 +44,7 @@ divSchemes
     div(U)          Gauss linear;
     div((Su*grad(b))) Gauss linear;
     div((U+((Su*Xi)*grad(b)))) Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -70,10 +70,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSolution b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSolution
index 14c6270cb..bb58ac024 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSolution
+++ b/tutorials/combustion/XiFoam/les/pitzDaily/system/fvSolution
@@ -122,7 +122,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nOuterCorrectors 2;
     nCorrectors     1;
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
index 9e60a7502..923b1c8b5 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSchemes
@@ -44,7 +44,7 @@ divSchemes
     div(U)          Gauss linear;
     div((Su*grad(b))) Gauss linear;
     div((U+((Su*Xi)*grad(b)))) Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -70,10 +70,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSolution b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSolution
index aa029ac96..f6268e03c 100644
--- a/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSolution
+++ b/tutorials/combustion/XiFoam/les/pitzDaily3D/system/fvSolution
@@ -122,7 +122,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nOuterCorrectors 2;
     nCorrectors     1;
diff --git a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
index fa004306d..96e263d29 100644
--- a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
+++ b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes
@@ -50,7 +50,7 @@ divSchemes
     div(U)          Gauss linear;
     div((Su*grad(b))) Gauss linear;
     div((U+((Su*Xi)*grad(b)))) Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -77,10 +77,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSolution b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSolution
index 7c6380a8e..ad27ea05d 100644
--- a/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSolution
+++ b/tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSolution
@@ -122,7 +122,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/combustion/engineFoam/kivaTest/constant/polyMesh/boundary b/tutorials/combustion/engineFoam/kivaTest/constant/polyMesh/boundary
index fd71ea09e..a50d04c25 100644
--- a/tutorials/combustion/engineFoam/kivaTest/constant/polyMesh/boundary
+++ b/tutorials/combustion/engineFoam/kivaTest/constant/polyMesh/boundary
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
diff --git a/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes b/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
index 6e689868b..15ad9d5f3 100644
--- a/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
+++ b/tutorials/combustion/engineFoam/kivaTest/system/fvSchemes
@@ -43,7 +43,7 @@ divSchemes
     div(U)          Gauss linear;
     div((Su*grad(b))) Gauss linear;
     div((U+((Su*Xi)*grad(b)))) Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -59,6 +59,7 @@ laplacianSchemes
     laplacian(alphaEff,h) Gauss linear limited 0.5;
     laplacian(alphaEff,hu) Gauss linear limited 0.5;
     laplacian(alphaEff,ft) Gauss linear limited 0.5;
+    laplacian(alphaEff,b) Gauss linear limited 0.5;
 }
 
 interpolationSchemes
@@ -71,10 +72,4 @@ snGradSchemes
     default         limited 0.5;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //