WIP: TUT: tutorial adjustment for backward compatibility

WIP: TUT: tutorial adjustments for backward compatibility
WIP: Updated lagrangian solvers followinf library changes
2023-01-31 16:14:43 +00:00 · 2023-01-20 16:56:11 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00
247 changed files with 3854 additions and 1848 deletions
--- a/applications/solvers/combustion/fireFoam/createClouds.H
+++ b/applications/solvers/combustion/fireFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
+reactingCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingCloud.H"
+#include "reactingCloud.H"
 #include "surfaceFilmModel.H"
 #include "pyrolysisModelCollection.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
@ -42,11 +42,11 @@ Description
 #include "CorrectPhi.H"

 #ifdef MPPIC
-    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCloud
+    #include "kinematicCloud.H"
+    #define kinematicTypeCloud kinematicCloud
 #else
-    #include "basicKinematicCollidingCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCollidingCloud
+    #include "kinematicCollidingCloud.H"
+    #define kinematicTypeCloud kinematicCollidingCloud
 #endif

 int main(int argc, char *argv[])
@ -88,7 +88,7 @@ int main(int argc, char *argv[])
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Store the particle positions
-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

@ -111,15 +111,15 @@ int main(int argc, char *argv[])
        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

-        kinematicCloud.evolve();
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -12,6 +12,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
@ -11,8 +11,11 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
@ -43,11 +43,11 @@ Description
 #include "pimpleControl.H"

 #ifdef MPPIC
-    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCloud
+    #include "kinematicCloud.H"
+    #define kinematicTypeCloud kinematicCloud
 #else
-    #include "basicKinematicCollidingCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCollidingCloud
+    #include "kinematicCollidingCloud.H"
+    #define kinematicTypeCloud kinematicCollidingCloud
 #endif

 int main(int argc, char *argv[])
@ -91,16 +91,16 @@ int main(int argc, char *argv[])
        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
--- a/applications/solvers/lagrangian/DPMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/Make/options
@ -10,8 +10,11 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/createFields.H
+++ b/applications/solvers/lagrangian/DPMFoam/createFields.H
@ -125,13 +125,13 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicTypeCloud kinematicCloud
+kinematicTypeCloud kCloud
 (
    kinematicCloudName,
+    g,
    rhoc,
    Uc,
-    muc,
-    g
+    muc
 );

 // Particle fraction upper limit
@ -139,13 +139,13 @@ scalar alphacMin
 (
    1.0
  - (
-        kinematicCloud.particleProperties().subDict("constantProperties")
+        kCloud.particleProperties().subDict("constantProperties")
       .get<scalar>("alphaMax")
    )
 );

 // Update alphac from the particle locations
-alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+alphac = max(1.0 - kCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -6,7 +6,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@ -35,7 +34,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lcoalCombustion\
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@ -37,8 +37,8 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicThermoCloud.H"
-#include "coalCloud.H"
+#include "thermoCloud.H"
+#include "reactingMultiphaseCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "fvOptions.H"
--- a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
@ -1,19 +1,19 @@
 Info<< "\nConstructing coal cloud" << endl;
-coalCloud coalParcels
+reactingMultiphaseCloud coalParcels
 (
    "coalCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );

 Info<< "\nConstructing limestone cloud" << endl;
-basicThermoCloud limestoneParcels
+thermoCloud limestoneParcels
 (
    "limestoneCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/createFields.H
@ -63,13 +63,13 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicCollidingCloud kinematicCloud
+kinematicCollidingCloud kCloud
 (
    kinematicCloudName,
+    g,
    rhoInf,
    U,
-    mu,
-    g
+    mu
 );

 IOobject Hheader
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
@ -41,7 +41,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -76,19 +76,19 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

        U.correctBoundaryConditions();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C
@ -40,7 +40,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -75,13 +75,13 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
@ -8,6 +8,12 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
--- a/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
@ -6,7 +6,7 @@
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
-        parcels.SU(U, true)
+        scalar(1)/rhoInfValue*parcels.SU(U)
      + fvOptions(U)
    );

--- a/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
@ -4,13 +4,12 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-
-basicKinematicCloud parcels
+kinematicCloud parcels
 (
    kinematicCloudName,
+    g,
    rhoInf,
    U,
-    muc,
-    g
+    muc
 );

--- a/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
@ -40,7 +40,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "surfaceFilmModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
@ -67,7 +67,6 @@ int main(int argc, char *argv[])
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
-    #include "createRegionControls.H"
    #include "createUfIfPresent.H"

    turbulence->validate();
@ -95,7 +94,7 @@ int main(int argc, char *argv[])
        // Do any mesh changes
        mesh.update();

-        if (solvePrimaryRegion && mesh.changing())
+        if (pimple.solveFlow() && mesh.changing())
        {
            MRF.update();

@ -120,7 +119,7 @@ int main(int argc, char *argv[])
        parcels.evolve();
        surfaceFilm.evolve();

-        if (solvePrimaryRegion)
+        if (pimple.solveFlow())
        {
            // --- PIMPLE loop
            while (pimple.loop())
--- a/applications/solvers/lagrangian/parcelFoam/Make/files
+++ b/applications/solvers/lagrangian/parcelFoam/Make/files
@ -0,0 +1,3 @@
+parcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/parcelFoam
--- a/applications/solvers/lagrangian/parcelFoam/Make/options
+++ b/applications/solvers/lagrangian/parcelFoam/Make/options
@ -1,49 +1,48 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude

-LIB_LIBS = \
+
+EXE_LIBS = \
    -lfiniteVolume \
+    -lfvOptions \
    -lmeshTools \
+    -lsampling \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -latmosphericModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsurfaceFilmDerivedFvPatchFields \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -ldistributionModels \
-    -lspecie \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
-    -lreactionThermophysicalModels \
-    -lSLGThermo \
-    -lradiationModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
-    -lregionModels \
-    -lsurfaceFilmModels \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
    -lregionFaModels \
-    -lfiniteArea
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/parcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/parcelFoam/UEqn.H
@ -0,0 +1,22 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U) + fvm::div(phi, U)
+      + MRF.DDt(U)
+      + turbulence->divDevReff(U)
+     ==
+        scalar(1)/rhoInfValue*parcels.SU(U)
+      + fvOptions(U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+
+        fvOptions.correct(U);
+    }
--- a/applications/solvers/lagrangian/parcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/parcelFoam/createClouds.H
@ -0,0 +1,2 @@
+auto parcels = parcelCloudModelList(g, rhoInf, U, muc);
+
--- a/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
@ -0,0 +1 @@
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/parcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFields.H
@ -0,0 +1,85 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+singlePhaseTransportModel laminarTransport(U, phi);
+
+dimensionedScalar rhoInfValue
+(
+    "rhoInf",
+    dimDensity,
+    laminarTransport
+);
+
+volScalarField rhoInf
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    rhoInfValue
+);
+
+volScalarField muc
+(
+    IOobject
+    (
+        "muc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    rhoInf*laminarTransport.nu()
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+mesh.setFluxRequired(p.name());
+
+#include "createMRF.H"
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+#include "createFvOptions.H"
+
+
--- a/applications/solvers/lagrangian/parcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/parcelFoam/pEqn.H
@ -0,0 +1,58 @@
+
+volScalarField rAU(1.0/UEqn.A());
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
+
+if (pimple.ddtCorr())
+{
+    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
+}
+
+MRF.makeRelative(phiHbyA);
+
+if (p.needReference())
+{
+    fvc::makeRelative(phiHbyA, U);
+    adjustPhi(phiHbyA, U, p);
+    fvc::makeAbsolute(phiHbyA, U);
+}
+
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, U, phiHbyA, rAU, MRF);
+
+// Non-orthogonal pressure corrector loop
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAU, p)
+     ==
+        fvc::div(phiHbyA)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA - pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+// Correct rhoUf if the mesh is moving
+fvc::correctUf(Uf, U, phi);
+
+// Make the fluxes relative to the mesh motion
+fvc::makeRelative(phi, U);
+
--- a/applications/solvers/lagrangian/parcelFoam/parcelFoam.C
+++ b/applications/solvers/lagrangian/parcelFoam/parcelFoam.C
@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    parcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for incompressible, turbulent flow with kinematic,
+    particle cloud, and surface film modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "surfaceFilmModel.H"
+#include "parcelCloudModelList.H"
+#include "fvOptions.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Transient solver for incompressible, turbulent flow"
+        " with kinematic particle clouds"
+        " and surface film modelling."
+    );
+
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createUfIfPresent.H"
+
+    turbulence->validate();
+
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+        #include "CourantNo.H"
+        #include "setMultiRegionDeltaT.H"
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Store the particle positions
+        parcels.storeGlobalPositions();
+
+        // Do any mesh changes
+        mesh.update();
+
+        if (pimple.solveFlow() && mesh.changing())
+        {
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & Uf();
+
+                #include "../../incompressible/pimpleFoam/correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, U);
+            }
+
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
+        }
+
+        parcels.evolve();
+        surfaceFilm.evolve();
+
+        if (pimple.solveFlow())
+        {
+            // --- PIMPLE loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    laminarTransport.correct();
+                    turbulence->correct();
+                }
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingTypeCloud parcels
+reactingTypeCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
@ -37,8 +37,8 @@ Description

 \*---------------------------------------------------------------------------*/

-#define CLOUD_BASE_TYPE HeterogeneousReacting
-#define CLOUD_BASE_TYPE_NAME "HeterogeneousReacting"
+#define CLOUD_BASE_TYPE heterogeneousReacting
+#define CLOUD_BASE_TYPE_NAME "heterogeneousReacting"

 #include "reactingParcelFoam.C"

--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@ -53,12 +53,12 @@ Description
 #include "cloudMacros.H"

 #ifndef CLOUD_BASE_TYPE
-    #define CLOUD_BASE_TYPE ReactingMultiphase
+    #define CLOUD_BASE_TYPE reactingMultiphase
    #define CLOUD_BASE_TYPE_NAME "reacting"
 #endif

 #include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
-#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -81,7 +81,6 @@ int main(int argc, char *argv[])
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
-    #include "createRegionControls.H"
    #include "initContinuityErrs.H"
    #include "createRhoUfIfPresent.H"

@ -105,7 +104,7 @@ int main(int argc, char *argv[])
        // so that it can be mapped and used in correctPhi
        // to ensure the corrected phi has the same divergence
        autoPtr<volScalarField> divrhoU;
-        if (solvePrimaryRegion && correctPhi)
+        if (pimple.solveFlow() && correctPhi)
        {
            divrhoU.reset
            (
@ -133,7 +132,7 @@ int main(int argc, char *argv[])

        // Store momentum to set rhoUf for introduced faces.
        autoPtr<volVectorField> rhoU;
-        if (solvePrimaryRegion && rhoUf.valid())
+        if (pimple.solveFlow() && rhoUf.valid())
        {
            rhoU.reset(new volVectorField("rhoU", rho*U));
        }
@ -144,7 +143,7 @@ int main(int argc, char *argv[])
        // Do any mesh changes
        mesh.update();

-        if (solvePrimaryRegion && mesh.changing())
+        if (pimple.solveFlow() && mesh.changing())
        {
            gh = (g & mesh.C()) - ghRef;
            ghf = (g & mesh.Cf()) - ghRef;
@ -172,7 +171,7 @@ int main(int argc, char *argv[])
        parcels.evolve();
        surfaceFilm.evolve();

-        if (solvePrimaryRegion)
+        if (pimple.solveFlow())
        {
            if (pimple.nCorrPIMPLE() <= 1)
            {
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -7,7 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing " << CLOUD_BASE_TYPE_NAME << " cloud" << endl;
-basicReactingTypeCloud parcels
+reactingTypeCloud parcels
 (
    word(CLOUD_BASE_TYPE_NAME) & "Cloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
@ -45,13 +45,13 @@ Description
 #include "cloudMacros.H"

 #ifndef CLOUD_BASE_TYPE
-    #define CLOUD_BASE_TYPE ReactingMultiphase
+    #define CLOUD_BASE_TYPE reactingMultiphase
    //#define CLOUD_BASE_TYPE_NAME "reactingMultiphase" Backwards compat
    #define CLOUD_BASE_TYPE_NAME "reacting"
 #endif

 #include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
-#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
@ -0,0 +1,41 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        rho*(U&g)
+      + parcels.Sh(he)
+      + surfaceFilm.Sh()
+      + radiation->Sh(thermo, he)
+      + Qdot
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/Make/files
+++ b/applications/solvers/lagrangian/rhoParcelFoam/Make/files
@ -0,0 +1,3 @@
+rhoParcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/rhoParcelFoam
--- a/applications/solvers/lagrangian/rhoParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/rhoParcelFoam/Make/options
@ -1,52 +1,58 @@
 EXE_INC = \
+    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-    -I$(LIB_SRC)/lagrangian/turbulence/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam

-LIB_LIBS = \
+EXE_LIBS = \
    -lfiniteVolume \
-    -lsurfMesh \
+    -lfvOptions \
    -lmeshTools \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -llagrangianTurbulence \
-    -ldistributionModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -lsampling \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lthermoTools \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
    -lreactionThermophysicalModels \
+    -lthermophysicalProperties \
    -lSLGThermo \
-    -lradiationModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
+    -lchemistryModel \
    -lregionModels \
+    -lradiationModels \
    -lsurfaceFilmModels \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
+    -lsurfaceFilmDerivedFvPatchFields \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lODE \
+    -lcombustionModels \
    -lregionFaModels \
-    -lfiniteArea
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
@ -0,0 +1,34 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        parcels.SU(U)
+      + fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }
--- a/applications/solvers/lagrangian/rhoParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/YEqn.H
@ -0,0 +1,51 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+{
+    combustion->correct();
+    Qdot = combustion->Qdot();
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+              + fvOptions(rho, Yi)
+              + combustion->R(Yi)
+              + surfaceFilm.Srho(i)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createClouds.H
@ -0,0 +1,3 @@
+Info<< "\nConstructing cloud" << endl;
+
+auto parcels = parcelCloudModelList(g, rho, U, slgThermo);
--- a/applications/solvers/lagrangian/rhoParcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createFieldRefs.H
@ -0,0 +1,5 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+const label inertIndex(composition.species().find(inertSpecie));
+
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/rhoParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createFields.H
@ -0,0 +1,132 @@
+#include "createRDeltaT.H"
+
+Info<< "Reading thermophysical properties\n" << endl;
+autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
+rhoReactionThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+const word inertSpecie(thermo.get<word>("inertSpecie"));
+if
+(
+    !composition.species().found(inertSpecie)
+  && composition.species().size() > 0
+)
+{
+    FatalIOErrorIn(args.executable().c_str(), thermo)
+        << "Inert specie " << inertSpecie << " not found in available species "
+        << composition.species()
+        << exit(FatalIOError);
+}
+
+Info<< "Creating field rho\n" << endl;
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+volScalarField& p = thermo.p();
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating combustion model\n" << endl;
+autoPtr<CombustionModel<rhoReactionThermo>> combustion
+(
+    CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
+);
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
+mesh.setFluxRequired(p_rgh.name());
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField Qdot
+(
+    IOobject
+    (
+        "Qdot",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
+);
+
+#include "createDpdt.H"
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+#include "createFvOptions.H"
--- a/applications/solvers/lagrangian/rhoParcelFoam/createMeshesPostProcess.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createMeshesPostProcess.H
@ -0,0 +1,35 @@
+#include "createMesh.H"
+
+dictionary filmDict;
+
+IOobject io
+(
+    "surfaceFilmProperties",
+    mesh.time().constant(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::NO_WRITE,
+    false
+);
+
+if (io.typeHeaderOk<IOdictionary>())
+{
+    IOdictionary propDict(io);
+
+    filmDict = std::move(propDict);
+
+    const word filmRegionName = filmDict.get<word>("region");
+
+    fvMesh filmMesh
+    (
+        IOobject
+        (
+            filmRegionName,
+            runTime.timeName(),
+            runTime,
+            IOobject::MUST_READ
+        )
+    );
+}
+
+
--- a/applications/solvers/lagrangian/rhoParcelFoam/createSurfaceFilmModel.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createSurfaceFilmModel.H
@ -0,0 +1,6 @@
+Info<< "\nConstructing surface film model" << endl;
+
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
+(
+    regionModels::surfaceFilmModel::New(mesh, g)
+);
--- a/applications/solvers/lagrangian/rhoParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/pEqn.H
@ -0,0 +1,96 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
+    )
+  + phig
+);
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+ ==
+    parcels.Srho()
+  + surfaceFilm.Srho()
+  + fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix p_rghEqn
+    (
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
+    );
+
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }
+}
+
+p = p_rgh + rho*gh;
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+    p_rgh = p - rho*gh;
+}
+else if (pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Correct rhoUf if the mesh is moving
+fvc::correctRhoUf(rhoUf, rho, U, phi);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+
+    if (mesh.moving())
+    {
+        dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
+    }
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/rhoEqn.H
@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho(rho)
+      + surfaceFilm.Srho()
+      + fvOptions(rho)
+    );
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+}
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/rhoParcelFoam/rhoParcelFoam.C
+++ b/applications/solvers/lagrangian/rhoParcelFoam/rhoParcelFoam.C
@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2020 OpenFOAM Foundation
+    Copyright (C) 2018-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoParcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for compressible, turbulent flow with a reacting,
+    multiphase particle cloud, and surface film modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "turbulentFluidThermoModel.H"
+#include "surfaceFilmModel.H"
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
+#include "radiationModel.H"
+#include "SLGThermo.H"
+#include "fvOptions.H"
+#include "pimpleControl.H"
+#include "pressureControl.H"
+#include "CorrectPhi.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "parcelCloudModelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Transient solver for compressible, turbulent flow"
+        " with lagrangian parcel clouds and surface film modelling."
+    );
+
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "initContinuityErrs.H"
+    #include "createRhoUfIfPresent.H"
+
+    turbulence->validate();
+
+    if (!LTS)
+    {
+        #include "compressibleCourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+
+        // Store divrhoU from the previous mesh
+        // so that it can be mapped and used in correctPhi
+        // to ensure the corrected phi has the same divergence
+        autoPtr<volScalarField> divrhoU;
+        if (pimple.solveFlow() && correctPhi)
+        {
+            divrhoU.reset
+            (
+                new volScalarField
+                (
+                    "divrhoU",
+                    fvc::div(fvc::absolute(phi, rho, U))
+                )
+            );
+        }
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "compressibleCourantNo.H"
+            #include "setMultiRegionDeltaT.H"
+        }
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Store momentum to set rhoUf for introduced faces.
+        autoPtr<volVectorField> rhoU;
+        if (pimple.solveFlow() && rhoUf.valid())
+        {
+            rhoU.reset(new volVectorField("rhoU", rho*U));
+        }
+
+        // Store the particle positions
+        parcels.storeGlobalPositions();
+
+        // Do any mesh changes
+        mesh.update();
+
+        if (pimple.solveFlow() && mesh.changing())
+        {
+            gh = (g & mesh.C()) - ghRef;
+            ghf = (g & mesh.Cf()) - ghRef;
+
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & rhoUf();
+
+                #include "../../compressible/rhoPimpleFoam/correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, rho, U);
+            }
+
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
+        }
+
+        parcels.evolve();
+        surfaceFilm.evolve();
+
+        if (pimple.solveFlow())
+        {
+            if (pimple.nCorrPIMPLE() <= 1)
+            {
+                #include "rhoEqn.H"
+            }
+
+            // --- PIMPLE loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+                #include "YEqn.H"
+                #include "EEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+
+            rho = thermo.rho();
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -23,38 +23,33 @@ License
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

-Class
-    Foam::basicThermoCloud
+Global
+    setMultiRegionDeltaT

 Description
-    Cloud class to introduce thermodynamic parcels
+    Reset the timestep to maintain a constant maximum Courant numbers.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.

 \*---------------------------------------------------------------------------*/

-#ifndef basicThermoCloud_H
-#define basicThermoCloud_H
-
-#include "ThermoCloud.H"
-#include "basicThermoParcel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
+if (adjustTimeStep)
 {
-    typedef ThermoCloud
-    <
-        KinematicCloud
-        <
-            Cloud
-            <
-                basicThermoParcel
-            >
-        >
-    > basicThermoCloud;
+    const scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
+    const scalar deltaTFact =
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }

-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
 // ************************************************************************* //
--- a/applications/solvers/lagrangian/rhoParcelFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/setRDeltaT.H
@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    // Maximum flow Courant number
+    scalar maxCo(pimpleDict.get<scalar>("maxCo"));
+
+    // Maximum time scale
+    scalar maxDeltaT(pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT));
+
+    // Smoothing parameter (0-1) when smoothing iterations > 0
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    // Damping coefficient (1-0)
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTDampingCoeff", 0.2)
+    );
+
+    // Maximum change in cell temperature per iteration
+    // (relative to previous value)
+    scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
+
+
+    Info<< "Time scales min/max:" << endl;
+
+    // Cache old reciprocal time scale field
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Flow time scale
+    {
+        rDeltaT.ref() =
+        (
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V()*rho())
+        );
+
+        // Limit the largest time scale
+        rDeltaT.max(1/maxDeltaT);
+
+        Info<< "    Flow        = "
+            << gMin(1/rDeltaT.primitiveField()) << ", "
+            << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+
+    // Reaction source time scale
+    {
+        volScalarField::Internal rDeltaTT
+        (
+            mag
+            (
+                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+              + Qdot()
+            )
+           /(
+               alphaTemp
+              *rho()
+              *thermo.Cp()()()
+              *T()
+           )
+        );
+
+        Info<< "    Temperature = "
+            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            rDeltaTT
+        );
+    }
+
+    // Update the boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    // Spatially smooth the time scale field
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+    }
+
+    Info<< "    Overall     = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
@ -7,7 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@ -47,7 +46,6 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lsampling \
-    -lcoalCombustion \
    -lregionFaModels \
    -lfiniteArea \
    -lfaOptions
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing coal cloud" << endl;
-coalCloud parcels
+reactingMultiphaseCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/simpleCoalParcelFoam.C
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/simpleCoalParcelFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "coalCloud.H"
+#include "reactingMultiphaseCloud.H"
 #include "rhoReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -35,7 +35,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/createClouds.H
+++ b/applications/solvers/lagrangian/sprayFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicSprayCloud parcels
+sprayCloud parcels
 (
    "sprayCloud",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -39,7 +39,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/engineFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/engineFoam.C
@ -37,7 +37,7 @@ Description
 #include "engineTime.H"
 #include "engineMesh.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
@ -34,7 +34,6 @@ EXE_LIBS = \
    -lthermoTools \
    -llagrangian \
    -llagrangianIntermediate \
-    -llagrangianSpray \
    -llagrangianTurbulence \
    -lspecie \
    -lcompressibleTransportModels \
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/simpleSprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/simpleSprayFoam.C
@ -35,7 +35,7 @@ Description

 \*---------------------------------------------------------------------------*/

-#define CLOUD_BASE_TYPE Spray
+#define CLOUD_BASE_TYPE spray
 #define CLOUD_BASE_TYPE_NAME "spray"

 #include "simpleReactingParcelFoam.C"
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -40,7 +40,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
@ -38,7 +38,7 @@ Description
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
 #include "turbulenceModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
@ -55,11 +55,11 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicCloud kinematicCloud
+kinematicCloud kCloud
 (
    kinematicCloudName,
+    g,
    rho,
    U,
-    thermo.mu(),
-    g
+    thermo.mu()
 );
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
@ -37,7 +37,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -73,14 +73,14 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

        U.correctBoundaryConditions();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
@ -39,7 +39,7 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -75,8 +75,8 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/multiphase/MPPICInterFoam/MPPICInterFoam.C
+++ b/applications/solvers/multiphase/MPPICInterFoam/MPPICInterFoam.C
@ -52,7 +52,7 @@ Description
 #include "CorrectPhi.H"
 #include "fvcSmooth.H"

-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -111,12 +111,12 @@ int main(int argc, char *argv[])

        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();

         Info<< "Continuous phase-1 volume fraction = "
@ -130,7 +130,7 @@ int main(int argc, char *argv[])
        alphaPhic = alphacf*phi;
        alphacRho = alphac*rho;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(U));
+        fvVectorMatrix cloudSU(kCloud.SU(U));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -12,11 +12,14 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
--- a/applications/solvers/multiphase/MPPICInterFoam/createFields.H
+++ b/applications/solvers/multiphase/MPPICInterFoam/createFields.H
@ -147,13 +147,13 @@ volScalarField alphacRho(alphac*rho);
 alphacRho.oldTime();

 Info<< "Constructing kinematicCloud " << endl;
-basicKinematicCloud kinematicCloud
+kinematicCloud kCloud
 (
    "kinematicCloud",
+    g,
    rho,
    U,
-    mu,
-    g
+    mu
 );

 // Particle fraction upper limit
@ -161,13 +161,13 @@ scalar alphacMin
 (
    1.0
  - (
-        kinematicCloud.particleProperties().subDict("constantProperties")
+        kCloud.particleProperties().subDict("constantProperties")
       .get<scalar>("alphaMax")
    )
 );

 // Update alphac from the particle locations
-alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+alphac = max(1.0 - kCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/src/functionObjects/lagrangian/Make/options
+++ b/src/functionObjects/lagrangian/Make/options
@ -12,7 +12,14 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude

 LIB_LIBS = \
    -lfiniteVolume \
--- a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
+++ b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
@ -28,7 +28,7 @@ License

 #include "cloudInfo.H"
 #include "cloud.H"
-#include "kinematicCloud.H"
+#include "parcelCloudModel.H"
 #include "dictionary.H"
 #include "mathematicalConstants.H"
 #include "PstreamReduceOps.H"
@ -147,7 +147,7 @@ bool Foam::functionObjects::cloudInfo::performAction(unsigned request)

        const word& cloudName = names()[cloudi];

-        const auto* kinCloudPtr = obr_.cfindObject<kinematicCloud>(cloudName);
+        const auto* kinCloudPtr = obr_.cfindObject<parcelCloudModel>(cloudName);

        if (!kinCloudPtr)
        {
--- a/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C
+++ b/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C
@ -92,10 +92,10 @@ Foam::functionObjects::icoUncoupledKinematicCloud::icoUncoupledKinematicCloud
    kinematicCloud_
    (
        kinematicCloudName_,
+        g_,
        rho_,
        U_,
-        mu_,
-        g_
+        mu_
    )
 {}

--- a/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.H
+++ b/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.H
@ -76,7 +76,7 @@ SourceFiles
 #define functionObjects_icoUncoupledKinematicCloud_H

 #include "fvMeshFunctionObject.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"
 #include "volFields.H"
 #include "uniformDimensionedFields.H"

@ -122,7 +122,7 @@ class icoUncoupledKinematicCloud
        word kinematicCloudName_;

        //- The kinematic cloud
-        basicKinematicCollidingCloud kinematicCloud_;
+        kinematicCollidingCloud kinematicCloud_;


    // Private member functions
--- a/src/lagrangian/Allwmake
+++ b/src/lagrangian/Allwmake
@ -9,9 +9,7 @@ wmake $targetType basic
 wmake $targetType solidParticle
 wmake $targetType intermediate
 wmake $targetType turbulence
-wmake $targetType spray
 wmake $targetType DSMC
-wmake $targetType coalCombustion

 molecularDynamics/Allwmake $targetType $*

--- a/src/lagrangian/coalCombustion/Make/files
+++ b/src/lagrangian/coalCombustion/Make/files
@ -1,5 +0,0 @@
-/* Coal parcel and sub-models */
-coalParcel/makeCoalParcelSubmodels.C
-coalCloudList/coalCloudList.C
-
-LIB = $(FOAM_LIBBIN)/libcoalCombustion
--- a/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
@ -1,264 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2012-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvMatrices.H"
-#include "volFields.H"
-#include "DimensionedField.H"
-
-Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
-Foam::coalCloudList::UTrans() const
-{
-    tmp<volVectorField::Internal> tfld
-    (
-        new volVectorField::Internal
-        (
-            IOobject
-            (
-                "UTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedVector(dimMass*dimVelocity, Zero)
-        )
-    );
-
-    volVectorField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).UTrans();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
-(
-    volVectorField& U
-) const
-{
-    tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
-    fvVectorMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).SU(U);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::hsTrans() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "hsTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimEnergy, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).hsTrans();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
-(
-    volScalarField& hs
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).Sh(hs);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
-(
-    const label ii,
-    volScalarField& Yi
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).SYi(ii, Yi);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::rhoTrans() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimMass, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        forAll(operator[](i).rhoTrans(), j)
-        {
-            fld += operator[](i).rhoTrans()[j];
-        }
-    }
-
-    return tfld;
-}
-
-
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::Srho() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimDensity/dimTime, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).Srho();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::Srho
-(
-    const label i
-) const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimDensity/dimTime, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, j)
-    {
-        fld += operator[](j).Srho(i);
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
-(
-    volScalarField& rho
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).Srho(rho);
-    }
-
-    return tfvm;
-}
-
-
-// ************************************************************************* //
--- a/src/lagrangian/coalCombustion/coalParcel/coalParcel.H
+++ b/src/lagrangian/coalCombustion/coalParcel/coalParcel.H
@ -1,74 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::coalParcel
-
-Description
-    Definition of coal parcel
-
-SourceFiles
-    CoalParcel.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef coalParcel_H
-#define coalParcel_H
-
-#include "contiguous.H"
-#include "particle.H"
-#include "KinematicParcel.H"
-#include "ThermoParcel.H"
-#include "ReactingParcel.H"
-#include "ReactingMultiphaseParcel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef ReactingMultiphaseParcel
-    <
-        ReactingParcel
-        <
-            ThermoParcel
-            <
-                KinematicParcel
-                <
-                    particle
-                >
-            >
-        >
-    > coalParcel;
-
-    //- Non-contiguous data for coalParcel
-    template<> struct is_contiguous<coalParcel> : std::false_type {};
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
--- a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
+++ b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
@ -1,59 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef makeCoalParcelSurfaceReactionModels_H
-#define makeCoalParcelSurfaceReactionModels_H
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "NoSurfaceReaction.H"
-#include "COxidationDiffusionLimitedRate.H"
-#include "COxidationIntrinsicRate.H"
-#include "COxidationKineticDiffusionLimitedRate.H"
-#include "COxidationHurtMitchell.H"
-#include "COxidationMurphyShaddix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define makeCoalParcelSurfaceReactionModels(CloudType)                         \
-                                                                               \
-    makeSurfaceReactionModelType(COxidationDiffusionLimitedRate, CloudType);   \
-    makeSurfaceReactionModelType                                               \
-    (                                                                          \
-        COxidationKineticDiffusionLimitedRate,                                 \
-        CloudType                                                              \
-    );                                                                         \
-    makeSurfaceReactionModelType(COxidationIntrinsicRate, CloudType);          \
-    makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);           \
-    makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@ -2,59 +2,75 @@ PARCELS=parcels
 BASEPARCELS=$(PARCELS)/baseClasses
 DERIVEDPARCELS=$(PARCELS)/derived

+
 CLOUDS=clouds
-BASECLOUDS=$(CLOUDS)/baseClasses
+TEMPLATECLOUDS=$(CLOUDS)/Templates
+$(TEMPLATECLOUDS)/KinematicCloud/KinematicCloudName.C
+$(TEMPLATECLOUDS)/ThermoCloud/ThermoCloudName.C
+$(TEMPLATECLOUDS)/ReactingCloud/ReactingCloudName.C
+$(TEMPLATECLOUDS)/ReactingMultiphaseCloud/ReactingMultiphaseCloudName.C
+$(TEMPLATECLOUDS)/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudName.C
+$(TEMPLATECLOUDS)/CollidingCloud/CollidingCloudName.C
+$(TEMPLATECLOUDS)/MPPICCloud/MPPICCloudName.C
+$(TEMPLATECLOUDS)/SprayCloud/SprayCloudName.C
+
+
 DERIVEDCLOUDS=$(CLOUDS)/derived
-
-
-/* Cloud base classes */
-$(BASECLOUDS)/kinematicCloud/kinematicCloud.C
-$(BASECLOUDS)/thermoCloud/thermoCloud.C
-$(BASECLOUDS)/reactingCloud/reactingCloud.C
-$(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
-$(BASECLOUDS)/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
+$(DERIVEDCLOUDS)/kinematicCloud/kinematicCloud.C
+$(DERIVEDCLOUDS)/thermoCloud/thermoCloud.C
+$(DERIVEDCLOUDS)/reactingCloud/reactingCloud.C
+$(DERIVEDCLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+$(DERIVEDCLOUDS)/heterogeneousReactingCloud/heterogeneousReactingCloud.C
+$(DERIVEDCLOUDS)/kinematicCollidingCloud/kinematicCollidingCloud.C
+$(DERIVEDCLOUDS)/kinematicMPPICCloud/kinematicMPPICCloud.C
+$(DERIVEDCLOUDS)/sprayCloud/sprayCloud.C


 /* kinematic parcel sub-models */
-KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
-$(KINEMATICPARCEL)/defineBasicKinematicParcel.C
-$(KINEMATICPARCEL)/makeBasicKinematicParcelSubmodels.C
+KINEMATICPARCEL=$(DERIVEDPARCELS)/kinematicParcel
+$(KINEMATICPARCEL)/defineKinematicParcel.C
+$(KINEMATICPARCEL)/makeKinematicParcelSubmodels.C


 /* kinematic colliding parcel sub-models */
-KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/basicKinematicCollidingParcel
-$(KINEMATICCOLLIDINGPARCEL)/defineBasicKinematicCollidingParcel.C
-$(KINEMATICCOLLIDINGPARCEL)/makeBasicKinematicCollidingParcelSubmodels.C
+KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/kinematicCollidingParcel
+$(KINEMATICCOLLIDINGPARCEL)/defineKinematicCollidingParcel.C
+$(KINEMATICCOLLIDINGPARCEL)/makeKinematicCollidingParcelSubmodels.C


 /* thermo parcel sub-models */
-THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
-$(THERMOPARCEL)/defineBasicThermoParcel.C
-$(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
+THERMOPARCEL=$(DERIVEDPARCELS)/thermoParcel
+$(THERMOPARCEL)/defineThermoParcel.C
+$(THERMOPARCEL)/makeThermoParcelSubmodels.C


 /* reacting parcel sub-models */
-REACTINGPARCEL=$(DERIVEDPARCELS)/basicReactingParcel
-$(REACTINGPARCEL)/defineBasicReactingParcel.C
-$(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
+REACTINGPARCEL=$(DERIVEDPARCELS)/reactingParcel
+$(REACTINGPARCEL)/defineReactingParcel.C
+$(REACTINGPARCEL)/makeReactingParcelSubmodels.C


 /* reacting multiphase parcel sub-models */
-REACTINGMPPARCEL=$(DERIVEDPARCELS)/basicReactingMultiphaseParcel
-$(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
-$(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
+REACTINGMPPARCEL=$(DERIVEDPARCELS)/reactingMultiphaseParcel
+$(REACTINGMPPARCEL)/defineReactingMultiphaseParcel.C
+$(REACTINGMPPARCEL)/makeReactingMultiphaseParcelSubmodels.C


-/* heterogeous reacting parcel sub-models */
-REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/basicHeterogeneousReactingParcel
-$(REACTINGHETERMPPARCEL)/defineBasicHeterogeneousReactingParcel.C
-$(REACTINGHETERMPPARCEL)/makeBasicHeterogeneousReactingParcelSubmodels.C
+/* heterogeneous reacting parcel sub-models */
+REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/heterogeneousReactingParcel
+$(REACTINGHETERMPPARCEL)/defineHeterogeneousReactingParcel.C
+$(REACTINGHETERMPPARCEL)/makeHeterogeneousReactingParcelSubmodels.C


 /* kinematic MPPIC parcel sub-models */
-KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel
-$(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
-$(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
+KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/kinematicMPPICParcel
+$(KINEMATICMPPICPARCEL)/defineKinematicMPPICParcel.C
+$(KINEMATICMPPICPARCEL)/makeKinematicMPPICParcelSubmodels.C
+
+/* spray parcel sub-models */
+SPRAYPARCEL=$(DERIVEDPARCELS)/sprayParcel
+$(SPRAYPARCEL)/defineSprayParcel.C
+$(SPRAYPARCEL)/makeSprayParcelSubmodels.C


 /* bolt-on models */
@ -130,4 +146,10 @@ conversion/ensight/ensightOutputCloud.C
 conversion/vtk/foamVtkLagrangianWriter.C


+clouds/baseCloudInterface.C
+clouds/parcelCloudModel.C
+clouds/parcelCloudModelNew.C
+clouds/parcelCloudModelList/parcelCloudModelList.C
+
+
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate
--- a/src/lagrangian/intermediate/clouds/BaseCloud.H
+++ b/src/lagrangian/intermediate/clouds/BaseCloud.H
@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BaseCloud
+
+Description
+    Base template for parcel clouds. Inserts the baseCloudInterface
+    virtualisation layer into the class. Also defines default zero-return
+    source methods to enable the less functional clouds to be used in more
+    complex situations.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_BaseCloud_H
+#define Foam_BaseCloud_H
+
+#include "Cloud.H"
+#include "baseCloudInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class SLGThermo;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+class BaseCloud
+:
+    public Cloud<ParticleType>,
+    virtual public baseCloudInterface
+{
+    // Private data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+
+public:
+
+    // Constructors
+
+        BaseCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            const bool readFields = true
+        )
+        :
+            Cloud<ParticleType>(rho.mesh(), cloudName, false),
+            mesh_(rho.mesh())
+        {}
+
+
+        BaseCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
+            const bool readFields = true
+        )
+        :
+            Cloud<ParticleType>(rho.mesh(), cloudName, false),
+            mesh_(rho.mesh())
+        {}
+
+        //- Copy constructor with new name
+        BaseCloud
+        (
+            BaseCloud<ParticleType>& c,
+            const word& name
+        )
+        :
+            Cloud<ParticleType>(c.mesh(), name, c),
+            mesh_(c.mesh_)
+        {}
+
+        //- Copy constructor with new name - creates bare cloud
+        BaseCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const BaseCloud<ParticleType>& c
+        )
+        :
+            Cloud<ParticleType>(mesh, name, IDLList<ParticleType>()),
+            mesh_(mesh)
+        {}
+
+    //- Destructor
+    virtual ~BaseCloud() = default;
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return reference to the mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+
+        // Sources
+
+            // Momentum
+
+                //- Momentum transfer [kg m/s]
+                virtual tmp<volVectorField::Internal> UTrans() const
+                {
+                    return volVectorField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), ":UTrans"),
+                        mesh_,
+                        dimensionedVector(dimMass*dimVelocity, Zero)
+                    );
+                }
+
+                //- Momentum transfer coefficient [kg]
+                virtual tmp<volScalarField::Internal> UCoeff() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "UCoeff"),
+                        mesh_,
+                        dimensionedScalar(dimMass, Zero)
+                    );
+                }
+
+                //- Return momentum source term [kg m/s^2]
+                virtual tmp<fvVectorMatrix> SU(volVectorField& U) const
+                {
+                    return tmp<fvVectorMatrix>::New(U, dimMass*dimAcceleration);
+                }
+
+
+            // Energy
+
+                //- Return sensible enthalpy transfer [J]
+                virtual tmp<volScalarField::Internal> hsTrans() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "hsTrans"),
+                        mesh_,
+                        dimensionedScalar(dimEnergy, Zero)
+                    );
+                }
+
+                //- Return sensible enthalpy transfer coefficient [J/K]
+                virtual tmp<volScalarField::Internal> hsCoeff() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "hsCoeff"),
+                        mesh_,
+                        dimensionedScalar(dimEnergy/dimTemperature, Zero)
+                    );
+                }
+
+                //- Return sensible enthalpy source term [J/s]
+                virtual tmp<fvScalarMatrix> Sh(volScalarField& hs) const
+                {
+                    return tmp<fvScalarMatrix>::New(hs, dimEnergy/dimTime);
+                }
+
+
+            // Radiation - overrides thermoCloud virtual abstract members
+
+                //- Return equivalent particulate emission [kg/m/s^3]
+                virtual tmp<volScalarField> Ep() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_Ep"),
+                        mesh_,
+                        dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+                    );
+                }
+
+                //- Return equivalent particulate absorption [1/m]
+                virtual tmp<volScalarField> ap() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_ap"),
+                        mesh_,
+                        dimensionedScalar(dimless/dimLength, Zero)
+                    );
+                }
+
+                //- Return equivalent particulate scattering factor [1/m]
+                virtual tmp<volScalarField> sigmap() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_sigmap"),
+                        mesh_,
+                        dimensionedScalar(dimless/dimLength, Zero)
+                    );
+                }
+
+
+            // Mass
+
+                //- Return mass source term for specie [kg/s]
+                virtual tmp<fvScalarMatrix> SYi
+                (
+                    const label speciei,
+                    volScalarField& Yi
+                ) const
+                {
+                    return tmp<fvScalarMatrix>::New(Yi, dimMass/dimTime);
+                }
+
+                //- Return total mass source term [kg/m^3/s]
+                virtual tmp<volScalarField::Internal> Srho() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "Srho"),
+                        mesh_,
+                        dimensionedScalar(dimDensity/dimTime, Zero)
+                    );
+                }
+
+                //- Return total mass source term [kg/s]
+                virtual tmp<fvScalarMatrix> Srho(volScalarField& rho) const
+                {
+                    return tmp<fvScalarMatrix>::New(rho, dimMass/dimTime);
+                }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/baseClasses/thermoCloud/thermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/baseClasses/thermoCloud/thermoCloud.H
@ -5,7 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020-2021 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -24,65 +25,67 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::thermoCloud
+    Foam::ParcelCloudModel

 Description
-    Virtual abstract base class for templated ThermoCloud
+    Outermost template for parcel clouds. Adds the parcelCloudModel
+    virtualisation layer and forwards the methods required by that layer.

 SourceFiles
-    thermoCloud.C
+    ParcelCloudModel.C

 \*---------------------------------------------------------------------------*/

-#ifndef thermoCloud_H
-#define thermoCloud_H
+#ifndef Foam_ParcelCloudModel_H
+#define Foam_ParcelCloudModel_H

-#include "volFields.H"
+#include "parcelCloudModel.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 namespace Foam
 {

-/*---------------------------------------------------------------------------*\
-                       Class thermoCloud Declaration
-\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-class thermoCloud
+template<class CloudType>
+class ParcelCloudModel
+:
+    public parcelCloudModel,
+    public CloudType
 {
 public:

-    //- Runtime type information
-    TypeName("thermoCloud");
-
-
    // Constructors

-        //- Null constructor
-        thermoCloud() = default;
-
-
-    //- Destructor
-    virtual ~thermoCloud() = default;
+        //- Inherit base constructors
+        using CloudType::CloudType;


    // Member Functions

-        // Access
+        //- Evolve the cloud
+        virtual void evolveME() final
+        {
+            // CloudType::evolve();
+            if (this->solution().canEvolve())
+            {
+                typedef typename CloudType::particleType parcelType;

-            // Radiation
+                typename parcelType::trackingData td(*this);

-                //- Return tmp equivalent particulate emission
-                virtual tmp<volScalarField> Ep() const = 0;
+                this->solve(*this, td);
+            }
+        }

-                //- Return tmp equivalent particulate absorption
-                virtual tmp<volScalarField> ap() const = 0;
-
-                //- Return tmp equivalent particulate scattering factor
-                virtual tmp<volScalarField> sigmap() const = 0;
+        //- Call this before a topology change. Stores the particles global
+        //- positions in the database for use during mapping.
+        virtual void storeGlobalPositionsME() const final
+        {
+            CloudType::storeGlobalPositions();
+        }
 };

-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 } // End namespace Foam
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
@ -28,6 +28,7 @@ License
 #include "CollidingCloud.H"
 #include "CollisionModel.H"
 #include "NoCollision.H"
+#include "subCycleTime.H"

 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //

@ -88,14 +89,14 @@ template<class CloudType>
 Foam::CollidingCloud<CloudType>::CollidingCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, mu, g, false),
+    CloudType(cloudName, g, rho, U, mu, false),
    constProps_(this->particleProperties()),
    collisionModel_(nullptr)
 {
@ -148,13 +149,6 @@ Foam::CollidingCloud<CloudType>::CollidingCloud
 {}


-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::CollidingCloud<CloudType>::~CollidingCloud()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
@ -58,6 +58,8 @@ namespace Foam
 template<class CloudType>
 class CollisionModel;

+TemplateName(CollidingCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class CollidingCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -65,7 +67,8 @@ class CollisionModel;
 template<class CloudType>
 class CollidingCloud
 :
-    public CloudType
+    public CloudType,
+    public CollidingCloudName
 {
 public:

@ -109,8 +112,7 @@ protected:
        // References to the cloud sub-models

            //- Collision model
-            autoPtr<CollisionModel<CollidingCloud<CloudType>>>
-                collisionModel_;
+            autoPtr<CollisionModel<CollidingCloud<CloudType>>> collisionModel_;


        // Initialisation
@ -142,10 +144,10 @@ public:
        CollidingCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
            const volScalarField& mu,
-            const dimensionedVector& g,
            bool readFields = true
        );

@ -184,7 +186,7 @@ public:


    //- Destructor
-    virtual ~CollidingCloud();
+    virtual ~CollidingCloud() = default;


    // Member Functions
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloudName.C
@ -5,7 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -25,15 +26,13 @@ License

 \*---------------------------------------------------------------------------*/

-#include "coalCloud.H"
-#include "makeCoalParcelSurfaceReactionModels.H"
+#include "CollidingCloud.H"

-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

 namespace Foam
 {
-    makeCoalParcelSurfaceReactionModels(coalCloud);
+defineTypeNameAndDebug(CollidingCloudName, 0);
 }

-
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@ -30,6 +30,7 @@ License
 #include "integrationScheme.H"
 #include "interpolation.H"
 #include "subCycleTime.H"
+#include "SLGThermo.H"

 #include "InjectionModelList.H"
 #include "DispersionModel.H"
@ -342,15 +343,14 @@ template<class CloudType>
 Foam::KinematicCloud<CloudType>::KinematicCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(rho.mesh(), cloudName, false),
-    kinematicCloud(),
+    CloudType(cloudName, g, rho, U, mu, false),
    cloudCopyPtr_(nullptr),
    mesh_(rho.mesh()),
    particleProperties_
@ -480,6 +480,29 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
 }


+template<class CloudType>
+Foam::KinematicCloud<CloudType>::KinematicCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    KinematicCloud<CloudType>
+    (
+        cloudName,
+        g,
+        rho,
+        U,
+        thermo.thermo().mu(),
+        readFields
+    )
+{}
+
+
 template<class CloudType>
 Foam::KinematicCloud<CloudType>::KinematicCloud
 (
@ -488,7 +511,6 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
 )
 :
    CloudType(c.mesh_, name, c),
-    kinematicCloud(),
    cloudCopyPtr_(nullptr),
    mesh_(c.mesh_),
    particleProperties_(c.particleProperties_),
@ -561,8 +583,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
    const KinematicCloud<CloudType>& c
 )
 :
-    CloudType(mesh, name, IDLList<parcelType>()),
-    kinematicCloud(),
+    CloudType(mesh, name, c),
    cloudCopyPtr_(nullptr),
    mesh_(mesh),
    particleProperties_
@ -687,8 +708,8 @@ void Foam::KinematicCloud<CloudType>::restoreState()
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::resetSourceTerms()
 {
-    UTrans().field() = Zero;
-    UCoeff().field() = 0.0;
+    UTransRef().field() = Zero;
+    UCoeffRef().field() = 0.0;
 }


@ -725,8 +746,8 @@ void Foam::KinematicCloud<CloudType>::relaxSources
    const KinematicCloud<CloudType>& cloudOldTime
 )
 {
-    this->relax(UTrans_(), cloudOldTime.UTrans(), "U");
-    this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U");
+    this->relax(UTransRef(), cloudOldTime.UTrans_(), "U");
+    this->relax(UCoeffRef(), cloudOldTime.UCoeff_(), "U");
 }


--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@ -59,7 +59,6 @@ SourceFiles

 #include "particle.H"
 #include "Cloud.H"
-#include "kinematicCloud.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "Random.H"
@ -80,6 +79,8 @@ namespace Foam

 class integrationScheme;

+class SLGThermo;
+
 template<class CloudType>
 class InjectionModelList;

@ -105,6 +106,8 @@ template<class CloudType>
 class IsotropyModel;


+TemplateName(KinematicCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class KinematicCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -113,7 +116,7 @@ template<class CloudType>
 class KinematicCloud
 :
    public CloudType,
-    public kinematicCloud
+    public KinematicCloudName
 {
 public:

@ -309,15 +312,24 @@ public:


    // Constructors
+        KinematicCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            bool readFields = true
+        );

        //- Construct given carrier gas fields
        KinematicCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
-            const volScalarField& mu,
-            const dimensionedVector& g,
+            const SLGThermo& thermo,
            bool readFields = true
        );

@ -527,33 +539,27 @@ public:
                // Momentum

                    //- Return reference to momentum source
-                    inline volVectorField::Internal& UTrans();
+                    inline volVectorField::Internal& UTransRef();

                    //- Return const reference to momentum source
-                    inline const volVectorField::Internal&
-                        UTrans() const;
+                    inline tmp<volVectorField::Internal> UTrans() const;

                     //- Return coefficient for carrier phase U equation
-                    inline volScalarField::Internal& UCoeff();
+                    inline volScalarField::Internal& UCoeffRef();

                    //- Return const coefficient for carrier phase U equation
-                    inline const volScalarField::Internal&
-                        UCoeff() const;
+                    inline tmp<volScalarField::Internal> UCoeff() const;

                    //- Return tmp momentum source term (compressible)
-                    inline tmp<fvVectorMatrix> SU
-                    (
-                        volVectorField& U,
-                        bool incompressible = false
-                    ) const;
+                    inline tmp<fvVectorMatrix> SU(volVectorField& U) const;


        // Check

            //- Total number of parcels
-            virtual label nParcels() const
+            label nParcels() const
            {
-                return CloudType::nParcels();
+                return this->size();
            }

            //- Total mass in system
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@ -436,14 +436,14 @@ Foam::KinematicCloud<CloudType>::cellLengthScale() const

 template<class CloudType>
 inline Foam::DimensionedField<Foam::vector, Foam::volMesh>&
-Foam::KinematicCloud<CloudType>::UTrans()
+Foam::KinematicCloud<CloudType>::UTransRef()
 {
    return *UTrans_;
 }


 template<class CloudType>
-inline const Foam::DimensionedField<Foam::vector, Foam::volMesh>&
+inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
 Foam::KinematicCloud<CloudType>::UTrans() const
 {
    return *UTrans_;
@ -452,14 +452,14 @@ Foam::KinematicCloud<CloudType>::UTrans() const

 template<class CloudType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::KinematicCloud<CloudType>::UCoeff()
+Foam::KinematicCloud<CloudType>::UCoeffRef()
 {
    return *UCoeff_;
 }


 template<class CloudType>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
 Foam::KinematicCloud<CloudType>::UCoeff() const
 {
    return *UCoeff_;
@ -468,7 +468,7 @@ Foam::KinematicCloud<CloudType>::UCoeff() const

 template<class CloudType>
 inline Foam::tmp<Foam::fvVectorMatrix>
-Foam::KinematicCloud<CloudType>::SU(volVectorField& U, bool incompressible)
+Foam::KinematicCloud<CloudType>::SU(volVectorField& U)
 const
 {
    if (debug)
@ -479,38 +479,28 @@ const
            << max(UCoeff()).value() << endl;
    }

-    dimensionSet dim(dimForce);
-    if (incompressible)
-    {
-        dim.reset(dimForce/dimDensity);
-    }
-
    if (solution_.coupled())
    {
+        const auto& deltaT = this->db().time().deltaT();
+
        if (solution_.semiImplicit("U"))
        {
-            volScalarField::Internal
-                Vdt(mesh_.V()*this->db().time().deltaT());
-
-            if (incompressible)
-            {
-                Vdt.dimensions() *= dimDensity;
-            }
+            volScalarField::Internal Vdt(mesh_.V()*deltaT);

            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
-            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dim));
+            auto tfvm = tmp<fvVectorMatrix>::New(U, dimForce);
            fvVectorMatrix& fvm = tfvm.ref();

-            fvm.source() = -UTrans()/(this->db().time().deltaT());
+            fvm.source() = -UTrans()/deltaT;

            return tfvm;
        }
    }

-    return tmp<fvVectorMatrix>::New(U, dim);
+    return tmp<fvVectorMatrix>::New(U, dimForce);
 }


--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "KinematicCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(KinematicCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@ -31,6 +31,7 @@ License
 #include "DampingModel.H"
 #include "IsotropyModel.H"
 #include "TimeScaleModel.H"
+#include "SLGThermo.H"

 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //

@ -70,14 +71,14 @@ template<class CloudType>
 Foam::MPPICCloud<CloudType>::MPPICCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, mu, g, false),
+    CloudType(cloudName, g, rho, U, mu, false),
    packingModel_(nullptr),
    dampingModel_(nullptr),
    isotropyModel_(nullptr)
@ -102,6 +103,29 @@ Foam::MPPICCloud<CloudType>::MPPICCloud
 }


+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    MPPICCloud<CloudType>
+    (
+        cloudName,
+        g,
+        rho,
+        U,
+        thermo.thermo().mu(),
+        readFields
+    )
+{}
+
+
 template<class CloudType>
 Foam::MPPICCloud<CloudType>::MPPICCloud
 (
@ -131,13 +155,6 @@ Foam::MPPICCloud<CloudType>::MPPICCloud
 {}


-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::MPPICCloud<CloudType>::~MPPICCloud()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
@ -41,7 +41,6 @@ SourceFiles
 #ifndef MPPICCloud_H
 #define MPPICCloud_H

-#include "particle.H"
 #include "Cloud.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
@ -55,6 +54,8 @@ namespace Foam

 // Forward declaration of classes

+class SLGThermo;
+
 template<class CloudType>
 class PackingModel;

@ -64,6 +65,8 @@ class DampingModel;
 template<class CloudType>
 class IsotropyModel;

+TemplateName(MPPICCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class MPPICCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -71,7 +74,8 @@ class IsotropyModel;
 template<class CloudType>
 class MPPICCloud
 :
-    public CloudType
+    public CloudType,
+    public MPPICCloudName
 {
 public:

@ -114,12 +118,10 @@ protected:
            autoPtr<PackingModel<MPPICCloud<CloudType>>> packingModel_;

            //- Damping model
-            autoPtr<DampingModel<MPPICCloud<CloudType>>>
-                dampingModel_;
+            autoPtr<DampingModel<MPPICCloud<CloudType>>> dampingModel_;

            //- Exchange model
-            autoPtr<IsotropyModel<MPPICCloud<CloudType>>>
-                isotropyModel_;
+            autoPtr<IsotropyModel<MPPICCloud<CloudType>>> isotropyModel_;


        // Initialisation
@ -136,10 +138,21 @@ public:
        MPPICCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
            const volScalarField& mu,
+            bool readFields = true
+        );
+
+        //- Construct given carrier gas fields
+        MPPICCloud
+        (
+            const word& cloudName,
            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
            bool readFields = true
        );

@ -178,7 +191,7 @@ public:


    //- Destructor
-    virtual ~MPPICCloud();
+    virtual ~MPPICCloud() = default;


    // Member Functions
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MPPICCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(MPPICCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@ -91,15 +91,14 @@ template<class CloudType>
 Foam::ReactingCloud<CloudType>::ReactingCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
-    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, g, thermo, false),
-    reactingCloud(),
+    CloudType(cloudName, g, rho, U, thermo, false),
    cloudCopyPtr_(nullptr),
    constProps_(this->particleProperties()),
    compositionModel_(nullptr),
@ -155,7 +154,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
 )
 :
    CloudType(c, name),
-    reactingCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(c.constProps_),
    compositionModel_(c.compositionModel_->clone()),
@ -195,7 +193,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
 )
 :
    CloudType(mesh, name, c),
-    reactingCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(),
    compositionModel_(c.compositionModel_->clone()),
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@ -46,7 +46,11 @@ SourceFiles
 #ifndef ReactingCloud_H
 #define ReactingCloud_H

-#include "reactingCloud.H"
+#include "className.H"
+#include "autoPtr.H"
+#include "Cloud.H"
+#include "volFields.H"
+#include "fvMatricesFwd.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -55,12 +59,16 @@ namespace Foam

 // Forward declaration of classes

+class SLGThermo;
+
 template<class CloudType>
 class CompositionModel;

 template<class CloudType>
 class PhaseChangeModel;

+TemplateName(ReactingCloud);
+
 /*---------------------------------------------------------------------------*\
                      Class ReactingCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -69,7 +77,7 @@ template<class CloudType>
 class ReactingCloud
 :
    public CloudType,
-    public reactingCloud
+    public ReactingCloudName
 {
 public:

@ -160,9 +168,9 @@ public:
        ReactingCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
-            const dimensionedVector& g,
            const SLGThermo& thermo,
            bool readFields = true
        );
@ -236,16 +244,14 @@ public:
                //- Mass

                    //- Return reference to mass source for field i
-                    inline volScalarField::Internal&
-                        rhoTrans(const label i);
+                    inline volScalarField::Internal& rhoTrans(const label i);

                    //- Return const access to mass source fields
                    inline const PtrList<volScalarField::Internal>&
                        rhoTrans() const;

                    //- Return reference to mass source fields
-                    inline PtrList<volScalarField::Internal>&
-                        rhoTrans();
+                    inline PtrList<volScalarField::Internal>& rhoTrans();

                    //- Return mass source term for specie i - specie eqn
                    inline tmp<fvScalarMatrix> SYi
@ -255,8 +261,10 @@ public:
                    ) const;

                    //- Return tmp mass source for field i - fully explicit
-                    inline tmp<volScalarField::Internal>
-                        Srho(const label i) const;
+                    inline tmp<volScalarField::Internal> Srho
+                    (
+                        const label i
+                    ) const;

                    //- Return tmp total mass source for carrier phase
                    //  - fully explicit
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@ -25,6 +25,8 @@ License

 \*---------------------------------------------------------------------------*/

+#include "fvmSup.H"
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ReactingCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(ReactingCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
@ -61,15 +61,14 @@ template<class CloudType>
 Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
-    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, g, thermo, false),
-    reactingHeterogeneousCloud(),
+    CloudType(cloudName, g, rho, U, thermo, false),
    cloudCopyPtr_(nullptr),
    heterogeneousReactionModel_(nullptr)
 {
@ -94,7 +93,6 @@ Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
 )
 :
    CloudType(c, name),
-    reactingHeterogeneousCloud(),
    cloudCopyPtr_(nullptr),
    heterogeneousReactionModel_(c.heterogeneousReactionModel_->clone())
 {}
@ -109,7 +107,6 @@ Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
 )
 :
    CloudType(mesh, name, c),
-    reactingHeterogeneousCloud(),
    cloudCopyPtr_(nullptr),
    heterogeneousReactionModel_(c.heterogeneousReactionModel_->clone())
 {}
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.H
@ -44,7 +44,11 @@ SourceFiles
 #ifndef ReactingHeterogeneousCloud_H
 #define ReactingHeterogeneousCloud_H

-#include "reactingHeterogeneousCloud.H"
+#include "className.H"
+#include "autoPtr.H"
+#include "Cloud.H"
+#include "volFields.H"
+#include "fvMatricesFwd.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -52,19 +56,23 @@ namespace Foam
 {

 // Forward declaration of classes
+
+class SLGThermo;
+
 template<class CloudType>
 class HeterogeneousReactingModel;

+TemplateName(ReactingHeterogeneousCloud);

 /*---------------------------------------------------------------------------*\
-                      Class ReactingHeterogeneousCloud Declaration
+                 Class ReactingHeterogeneousCloud Declaration
 \*---------------------------------------------------------------------------*/

 template<class CloudType>
 class ReactingHeterogeneousCloud
 :
    public CloudType,
-    public reactingHeterogeneousCloud
+    public ReactingHeterogeneousCloudName
 {
 public:

@ -136,9 +144,9 @@ public:
        ReactingHeterogeneousCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
-            const dimensionedVector& g,
            const SLGThermo& thermo,
            bool readFields = true
        );
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ReactingHeterogeneousCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(ReactingHeterogeneousCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@ -78,15 +78,14 @@ template<class CloudType>
 Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
-    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, g, thermo, false),
-    reactingMultiphaseCloud(),
+    CloudType(cloudName, g, rho, U, thermo, false),
    cloudCopyPtr_(nullptr),
    constProps_(this->particleProperties()),
    devolatilisationModel_(nullptr),
@ -120,7 +119,6 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
 )
 :
    CloudType(c, name),
-    reactingMultiphaseCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(c.constProps_),
    devolatilisationModel_(c.devolatilisationModel_->clone()),
@ -139,7 +137,6 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
 )
 :
    CloudType(mesh, name, c),
-    reactingMultiphaseCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(),
    devolatilisationModel_(nullptr),
@ -149,13 +146,6 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
 {}


-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::ReactingMultiphaseCloud<CloudType>::~ReactingMultiphaseCloud()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Kutalmis Bercin	55d5f99047	WIP: TUT: tutorial adjustment for backward compatibility	2023-01-31 16:14:43 +00:00
Kutalmis Bercin	d439440453	WIP: TUT: tutorial adjustments for backward compatibility	2023-01-20 16:56:11 +00:00
Andrew Heather	6c66cc58a6	WIP: Updated lagrangian solvers followinf library changes	2023-01-10 13:35:42 +00:00
Andrew Heather	265b3dd2f8	WIP: deprecated lagrangian coal combustion - models added to reactingMultiphase	2023-01-10 13:35:42 +00:00
Andrew Heather	b17a9128a1	ENH: Added new parcelFoam solver based on kinematicParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	31d4148173	STYLE: lagrangian turbulence library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	13b8a3c02a	STYLE: lagrangian intermediate library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	715831610d	WIP: Added new rhoParcelFoam solver based on reactingParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	ab42d1a1e6	WIP: migrated coalCloudList to parcelCloudModelList	2023-01-10 13:35:42 +00:00
Andrew Heather	13cf489f63	WIP: clouds - moved spray cloud and parcels to main intermediate library	2023-01-10 13:35:42 +00:00
Andrew Heather	3277fb931b	WIP: clouds - updated dependent code	2023-01-10 13:35:40 +00:00
Andrew Heather	9e83d31138	WIP: clouds - refactoring	2023-01-10 13:34:47 +00:00
Andrew Heather	878f1bf773	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	84348dc3f5	WIP: clouds - added debug - based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	1f7c038ad7	DUMMY	2023-01-10 13:34:47 +00:00
Andrew Heather	f20361d44a	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:33 +00:00
Andrew Heather	580a003d4d	WIP: Removed base classes	2023-01-10 13:34:22 +00:00
				`@ -0,0 +1,2 @@`
				`auto parcels = parcelCloudModelList(g, rhoInf, U, muc);`
				`@ -0,0 +1 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`