diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/Make/options b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/Make/options
index cffe6f344a148d8b60b35ccb42990542b8c1cfcc..2f17d01ffb408994955aba9bdbc6836635eec207 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/Make/options
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/Make/options
@@ -4,9 +4,10 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -llagrangian \
@@ -15,7 +16,7 @@ EXE_LIBS = \
     -lmeshTools \
     -lthermophysicalFunctions \
     -lbasicThermophysicalModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lspecie \
     -lradiation \
     -lcompressibleRASModels \
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/UEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/UEqn.H
index 4d58a14da89574c12d4162e343f98eba5cf7085f..5f2e597a152ef5c08f811d1f374326db15a7ed2a 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/UEqn.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/UEqn.H
@@ -4,8 +4,8 @@
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
      ==
-        thermoCloud1.SU1()
-      + kinematicCloud1.SU1()
+        thermoCloud1.SU()
+      + kinematicCloud1.SU()
       + rho.dimensionedInternalField()*g
     );
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createClouds.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createClouds.H
index bdf17e19efcdb715b5f7da5c01e96052ae41645f..ff41634a4ef48e45c45841861a35d08dd35a0f9e 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createClouds.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createClouds.H
@@ -5,7 +5,7 @@
         rho,
         U,
         g,
-        thermo()
+        thermo
     );
 
     Info<< "Constructing kinematicCloud1" << endl;
@@ -14,7 +14,7 @@
         "kinematicCloud1",
         rho,
         U,
-        thermo().mu(),
+        thermo.mu(),
         g
     );
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
index b07398a0314073ca5d677d4e7fc25362463dc03b..cebd44cfd4041042584f6706484083b3a4dad9b0 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -19,7 +20,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -48,7 +49,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H
index 5359c9c2cecae01b52923ff14943cc4a1373d598..2b3b60ce980cd23324343c7647de47caba0658cf 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H
@@ -6,12 +6,12 @@
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
         DpDt
-      + thermoCloud1.Sh1()
+      + thermoCloud1.Sh()
     );
 
     hEqn.relax();
 
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
index b506245034010d76f0ad0fb87dd22c8b559f5597..9443f909a356699185bfbc8497cf46e26fd1ed3b 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -35,8 +35,8 @@ if (transonic)
 else
 {
     phi =
-        fvc::interpolate(rho)*
-        (
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
index 1082636303c6891d00c76806b9d085de7b097980..77715fd48c1905eb01642ddefb02aaac1338e0db 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 #include "basicThermoCloud.H"
@@ -41,62 +41,59 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createClouds.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "createClouds.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        Info<< "Evolving thermoCloud1" << endl;
         thermoCloud1.evolve();
         thermoCloud1.info();
 
-        Info<< "Evolving kinematicCloud1" << endl;
         kinematicCloud1.evolve();
         kinematicCloud1.info();
 
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
         // --- PIMPLE loop
         for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
         {
-#           include "UEqn.H"
+            #include "UEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "hEqn.H"
-#               include "pEqn.H"
+                #include "hEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/epsilon b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/epsilon
index ec63ed16726be65c6d5561d1b31ab88fd3ef8581..9fcf588c0c292537815ea7ef2c72563a85a14851 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/epsilon
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/epsilon
@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
+    location    "0";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -22,28 +23,28 @@ boundaryField
 {
     top
     {
-        type            zeroGradient;
+        type            compressible::epsilonWallFunction;
+        value           uniform 5390.5;
     }
-
     bottom
     {
-        type            zeroGradient;
+        type            compressible::epsilonWallFunction;
+        value           uniform 5390.5;
     }
-
     walls
     {
-        type            zeroGradient;
+        type            compressible::epsilonWallFunction;
+        value           uniform 5390.5;
     }
-
     symmetry
     {
         type            symmetryPlane;
     }
-
     frontAndBack
     {
         type            empty;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/k b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/k
index b2e4f58aab813ed1f1ba9667db76c30d0ad49f1c..9abb57fdb013a1d9dd1ac7859f9554dc0c642b71 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/k
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/0/k
@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
+    location    "0";
     object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -22,28 +23,28 @@ boundaryField
 {
     top
     {
-        type            zeroGradient;
+        type            compressible::kQRWallFunction;
+        value           uniform 37.5;
     }
-
     bottom
     {
-        type            zeroGradient;
+        type            compressible::kQRWallFunction;
+        value           uniform 37.5;
     }
-
     walls
     {
-        type            zeroGradient;
+        type            compressible::kQRWallFunction;
+        value           uniform 37.5;
     }
-
     symmetry
     {
         type            symmetryPlane;
     }
-
     frontAndBack
     {
         type            empty;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
index 4a7f8b7b2b5d0f5f7c8e979a07185eb634c9a8c5..8c69fb42b5812561e5fc746eceec85ae1481a22b 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
@@ -21,20 +21,27 @@ DragModel       SphereDrag;
 
 DispersionModel StochasticDispersionRAS;
 
-WallInteractionModel StandardWallInteraction;
+PatchInteractionModel StandardWallInteraction;
 
-minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
-
-rho0            rho0 [ 1 -3 0 0 0 ] 5000;
+PostProcessingModel none;
 
 coupled         true;
 
+cellValueSourceCorrection on;
+
 parcelTypeId    2;
 
+constantProperties
+{
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 1e-15;
+    minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
+    rho0            rho0 [ 1 -3 0 0 0 ] 5000;
+}
+
 interpolationSchemes
 {
     rho             cell;
-    U               cellPointFace;
+    U               cellPoint;
     mu              cell;
 }
 
@@ -43,12 +50,21 @@ integrationSchemes
     U               Euler;
 }
 
+particleForces
+{
+    gravity         on;
+    virtualMass     off;
+    Cvm             0.5;
+    pressureGradient off;
+    gradU           gradU;
+}
+
 ManualInjectionCoeffs
 {
-    parcelBasisType mass;
     massTotal       massTotal [ 1 0 0 0 0 ] 0.0002;
+    parcelBasisType mass;
     SOI             0;
-    positionsFile   kinematicCloud1Positions;
+    positionsFile   "kinematicCloud1Positions";
     U0              ( 0 0 0 );
     parcelPDF
     {
@@ -70,29 +86,10 @@ ConeInjectionCoeffs
     position        ( 0.25 0.25 0.05 );
     direction       ( 0 -1 0 );
     parcelsPerSecond 10000;
-    volumeFlowRate  Constant;
-    volumeFlowRateCoeffs
-    {
-        value           0.01;
-    }
-
-    Umag            Constant;
-    UmagCoeffs
-    {
-        value           50;
-    }
-
-    thetaInner      Constant;
-    thetaInnerCoeffs
-    {
-        value           0;
-    }
-
-    thetaOuter      Constant;
-    thetaOuterCoeffs
-    {
-        value           30;
-    }
+    volumeFlowRate  Constant 0.01;
+    Umag            Constant 50;
+    thetaInner      Constant 0;
+    thetaOuter      Constant 30;
 
     parcelPDF
     {
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/polyMesh/boundary b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/polyMesh/boundary
index b8848e2d9cd3e1dbe2fa2a9aabe159a7108753f7..b130d7bc0aa5e48e9a785745d71564a71ac9a626 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/polyMesh/boundary
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/polyMesh/boundary
@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       polyBoundaryMesh;
+    location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
index fd752a082c9d60ef27277908ec50e7bf425a3970..a7c19443e21a00f8964ab6b167621c7cd8090d33 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
@@ -21,33 +21,37 @@ DragModel       SphereDrag;
 
 DispersionModel StochasticDispersionRAS;
 
-WallInteractionModel StandardWallInteraction;
+PatchInteractionModel StandardWallInteraction;
 
 HeatTransferModel RanzMarshall;
 
-radiation       off;
-
-minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
-
-rho0            rho0 [ 1 -3 0 0 0 ] 2500;
-
-T0              T0 [ 0 0 0 1 0 ] 300;
+PostProcessingModel none;
 
-cp0             cp0 [ 0 2 -2 -1 0 ] 900;
-
-epsilon0        epsilon0 [ 0 0 0 0 0 ] 1;
-
-f0              f0 [ 0 0 0 0 0 ] 0.5;
+radiation       off;
 
 coupled         true;
 
+cellValueSourceCorrection on;
+
 parcelTypeId    1;
 
+constantProperties
+{
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 1e-15;
+    TMin            TMin [ 0 0 0 1 0 ] 200;
+    minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
+    rho0            rho0 [ 1 -3 0 0 0 ] 2500;
+    T0              T0 [ 0 0 0 1 0 ] 300;
+    cp0             cp0 [ 0 2 -2 -1 0 ] 900;
+    epsilon0        epsilon0 [ 0 0 0 0 0 ] 1;
+    f0              f0 [ 0 0 0 0 0 ] 0.5;
+}
+
 interpolationSchemes
 {
     rho             cell;
-    U               cellPointFace;
     mu              cell;
+    U               cellPoint;
     T               cell;
     Cp              cell;
 }
@@ -58,12 +62,21 @@ integrationSchemes
     T               Analytical;
 }
 
+particleForces
+{
+    gravity         on;
+    virtualMass     off;
+    Cvm             0.5;
+    pressureGradient off;
+    gradU           gradU;
+}
+
 ManualInjectionCoeffs
 {
     massTotal       massTotal [ 1 0 0 0 0 ] 0.0001;
     parcelBasisType mass;
     SOI             0;
-    positionsFile   thermoCloud1Positions;
+    positionsFile   "thermoCloud1Positions";
     U0              ( 0 0 0 );
     parcelPDF
     {
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..78f28eca604be4d8fc94d38d23ae0cbb679b18fb
--- /dev/null
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //