diff --git a/applications/solvers/cfdemSolverIB/Make/options b/applications/solvers/cfdemSolverIB/Make/options
index a51d1dae..b841d1a3 100755
--- a/applications/solvers/cfdemSolverIB/Make/options
+++ b/applications/solvers/cfdemSolverIB/Make/options
@@ -1,24 +1,27 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude
+
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
diff --git a/applications/solvers/cfdemSolverIB/cfdemSolverIB.C b/applications/solvers/cfdemSolverIB/cfdemSolverIB.C
index 9c930e63..9c6e41be 100755
--- a/applications/solvers/cfdemSolverIB/cfdemSolverIB.C
+++ b/applications/solvers/cfdemSolverIB/cfdemSolverIB.C
@@ -39,7 +39,8 @@ Contributions
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloudIB.H"
 #include "implicitCouple.H"
@@ -52,10 +53,9 @@ Contributions
 
 #include "cellSet.H"
 
-#if defined(version22)
-    #include "meshToMeshNew.H"
-    #include "fvIOoptionList.H"
-#endif
+#include "meshToMeshNew.H"
+#include "fvIOoptionList.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,14 +66,14 @@ int main(int argc, char *argv[])
     #include "createTime.H"
 
     #include "createDynamicFvMesh.H"
+  
+    #include "createControl.H"
 
     #include "createFields.H"
 
     #include "initContinuityErrs.H"
 
-    #if defined(version22)
     #include "createFvOptions.H"
-    #endif
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -91,7 +91,6 @@ int main(int argc, char *argv[])
         interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
         mesh.update(); //dyM
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -107,45 +106,33 @@ int main(int argc, char *argv[])
                 fvm::ddt(voidfraction,U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
-                #if defined(version22)
                 ==
                 fvOptions(U)
-                #endif
             );
 
             UEqn.relax();
 
-            #if defined(version22)
             fvOptions.constrain(UEqn);
-            #endif
 
-            if (momentumPredictor)
+            if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+           while (piso.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
                 surfaceScalarField rUAf(fvc::interpolate(rUA));
 
                 U = rUA*UEqn.H();
-                #ifdef version23
+
                 phi = (fvc::interpolate(U) & mesh.Sf())
                     + rUAf*fvc::ddtCorr(U, phi);
-                #else
-                phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
-                #endif
+
                 adjustPhi(phi, U, p);
 
-                #if defined(version22)
-                fvOptions.relativeFlux(phi);
-                #endif
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
                     // Pressure corrector
 
@@ -156,20 +143,9 @@ int main(int argc, char *argv[])
 
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (piso.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }
diff --git a/applications/solvers/cfdemSolverPiso/Make/options b/applications/solvers/cfdemSolverPiso/Make/options
index 68aecd17..91f9f553 100644
--- a/applications/solvers/cfdemSolverPiso/Make/options
+++ b/applications/solvers/cfdemSolverPiso/Make/options
@@ -1,19 +1,22 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
index a6a5e4f0..e779aac4 100644
--- a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
@@ -36,7 +36,8 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
@@ -51,6 +52,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -67,7 +69,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -85,7 +86,7 @@ int main(int argc, char *argv[])
 
        //Force Checks
        vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
        reduce(fImpTotal, sumOp<vector>());
        Info << "TotalForceExp: " << fTotal << endl;
        Info << "TotalForceImp: " << fImpTotal << endl;
@@ -111,17 +112,14 @@ int main(int argc, char *argv[])
                 );
 
                 UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
+                if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
                     solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
+                else if (piso.momentumPredictor())
                     solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
                     volScalarField rUA = 1.0/UEqn.A();
 
@@ -131,13 +129,10 @@ int main(int argc, char *argv[])
 
                     U = rUA*UEqn.H();
 
-                    #ifdef version23
+
                     phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
                         + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
+
                     surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
                     surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
 
@@ -145,7 +140,7 @@ int main(int argc, char *argv[])
                         rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
 
                     // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                    while (piso.correctNonOrthogonal())
                     {
                         // Pressure corrector
                         fvScalarMatrix pEqn
@@ -154,20 +149,9 @@ int main(int argc, char *argv[])
                         );
                         pEqn.setReference(pRefCell, pRefValue);
 
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
+                        pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                        if (nonOrth == nNonOrthCorr)
+                        if (piso.finalNonOrthogonalIter())
                         {
                             phiGes -= pEqn.flux();
                             phi = phiGes;
diff --git a/applications/solvers/cfdemSolverPisoMS/Make/options b/applications/solvers/cfdemSolverPisoMS/Make/options
index 84c16e2c..473ab040 100644
--- a/applications/solvers/cfdemSolverPisoMS/Make/options
+++ b/applications/solvers/cfdemSolverPisoMS/Make/options
@@ -1,16 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I ../cfdemSolverPiso \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C b/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
index 8e8e65b2..0d8e5137 100644
--- a/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
+++ b/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
@@ -36,7 +36,8 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloudMS.H"
 #include "implicitCouple.H"
@@ -50,6 +51,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -67,7 +69,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -103,15 +104,12 @@ int main(int argc, char *argv[])
 
             UEqn.relax();
 
-            if (momentumPredictor)
+            if (piso.momentumPredictor())
                 solve(UEqn);
 
             // --- PISO loop
 
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-            for (int corr=0; corr<nCorrSoph; corr++)
+            while (piso.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
 
@@ -129,7 +127,7 @@ int main(int argc, char *argv[])
                     rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
 
                 // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correctNonOrthogonal())
                 {
                     // Pressure corrector
                     fvScalarMatrix pEqn
@@ -138,20 +136,9 @@ int main(int argc, char *argv[])
                     );
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (piso.finalNonOrthogonalIter())
                     {
                         phiGes -= pEqn.flux();
                     }
diff --git a/applications/solvers/cfdemSolverPisoScalar/Make/options b/applications/solvers/cfdemSolverPisoScalar/Make/options
index 68aecd17..91f9f553 100644
--- a/applications/solvers/cfdemSolverPisoScalar/Make/options
+++ b/applications/solvers/cfdemSolverPisoScalar/Make/options
@@ -1,19 +1,22 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
index 4d964424..d2621d91 100644
--- a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
+++ b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
@@ -36,7 +36,8 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
@@ -50,6 +51,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -64,7 +66,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -114,17 +115,14 @@ int main(int argc, char *argv[])
                 );
 
                 UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
+                if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
                     solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
+                else if (piso.momentumPredictor())
                     solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
                     volScalarField rUA = 1.0/UEqn.A();
 
@@ -134,21 +132,16 @@ int main(int argc, char *argv[])
 
                     U = rUA*UEqn.H();
 
-                    #ifdef version23
                     phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
                         + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
+                    
                     surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
                     surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
 
                     if (modelType=="A")
                         rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
 
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                    while (piso.correctNonOrthogonal())
                     {
                         // Pressure corrector
                         fvScalarMatrix pEqn
@@ -157,20 +150,9 @@ int main(int argc, char *argv[])
                         );
                         pEqn.setReference(pRefCell, pRefValue);
 
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
+                        pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                        if (nonOrth == nNonOrthCorr)
+                        if (piso.finalNonOrthogonalIter())
                         {
                             phiGes -= pEqn.flux();
                             phi = phiGes;