diff --git a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
index 3735f340bb..18631da635 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
@@ -76,7 +76,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
         )
     ),
     Ubar_(dict_.lookup("Ubar")),
-    gradPini_(readScalar(dict_.lookup("gradPini"))),
+    gradPini_(dict_.lookup("gradPini")),
     gradP_(gradPini_),
     flowDir_(Ubar_/mag(Ubar_)),
     cellSource_(dict_.lookup("cellSource")),
@@ -121,7 +121,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
         propsDict.lookup("gradient") >> gradP_;
     }
 
-    Info<< "    Initial pressure gradient = " << gradP_ << endl;
+    Info<< "    Initial pressure gradient = " << gradP_ << nl << endl;
 }
 
 
@@ -143,7 +143,7 @@ Foam::pressureGradientExplicitSource::Su() const
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedVector("zero", dimVelocity/dimTime, vector::zero)
+            dimensionedVector("zero", gradP_.dimensions(), vector::zero)
         )
     );
 
@@ -153,7 +153,7 @@ Foam::pressureGradientExplicitSource::Su() const
     {
         label cellI = iter.key();
 
-        sourceField[cellI] = flowDir_*gradP_;
+        sourceField[cellI] = flowDir_*gradP_.value();
     }
 
     return tSource;
@@ -201,10 +201,10 @@ void Foam::pressureGradientExplicitSource::update()
     }
 
     // Update pressure gradient
-    gradP_ += gradPplus;
+    gradP_.value() += gradPplus;
 
     Info<< "Uncorrected Ubar = " << magUbarAve << tab
-        << "Pressure gradient = " << gradP_ << endl;
+        << "Pressure gradient = " << gradP_.value() << endl;
 
     writeGradP();
 }
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H
index ce322332ee..5046df3f1c 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H
@@ -73,10 +73,10 @@ class pressureGradientExplicitSource
         vector Ubar_;
 
         //- Initial pressure gradient
-        scalar gradPini_;
+        dimensionedScalar gradPini_;
 
         //- Pressure gradient
-        scalar gradP_;
+        dimensionedScalar gradP_;
 
         //- Flow direction
         vector flowDir_;
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index cfd9bef857..8a8723a539 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -36,43 +36,32 @@ License
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>
-Foam::scalar Foam::KinematicCloud<ParcelType>::setNumberOfParticles
+void Foam::KinematicCloud<ParcelType>::addNewParcel
 (
-    const label nParcels,
-    const scalar pDiameter,
-    const scalar pVolumeFraction,
-    const scalar pRho,
-    const scalar pVolume
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
 )
 {
-    scalar nP = 0.0;
-    switch (parcelBasis_)
-    {
-        case pbMass:
-        {
-            nP = pVolumeFraction*massTotal_/nParcels
-               /(pRho*mathematicalConstant::pi/6.0*pow(pDiameter, 3));
-            break;
-        }
-        case pbNumber:
-        {
-            nP = pVolumeFraction*massTotal_/(pRho*pVolume);
-            break;
-        }
-        default:
-        {
-            nP = 0.0;
-            FatalErrorIn
-            (
-                "Foam::KinematicCloud<ParcelType>::setNumberOfParticles"
-                "(const label, const scalar, const scalar, const scalar, "
-                "const scalar)"
-            )<< "Unknown parcelBasis type" << nl
-             << exit(FatalError);
-        }
-    }
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        parcelTypeId_,
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        constProps_
+    );
 
-    return nP;
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
 }
 
 
@@ -107,14 +96,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
     rndGen_(label(0)),
-    time0_(this->db().time().value()),
-    parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
-    parcelBasis_(pbNumber),
-    massTotal_
-    (
-        dimensionedScalar(particleProperties_.lookup("massTotal")).value()
-    ),
-    massInjected_(0.0),
     rho_(rho),
     U_(U),
     mu_(mu),
@@ -160,9 +141,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             particleProperties_.subDict("integrationSchemes")
         )
     ),
-    nInjections_(0),
-    nParcelsAdded_(0),
-    nParcelsAddedTotal_(0),
     UTrans_
     (
         IOobject
@@ -191,27 +169,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
         mesh_,
         dimensionedScalar("zero",  dimensionSet(1, 0, -1, 0, 0), 0.0)
     )
-{
-    if (parcelBasisType_ == "mass")
-    {
-        parcelBasis_ = pbMass;
-    }
-    else if (parcelBasisType_ == "number")
-    {
-        parcelBasis_ = pbNumber;
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "Foam::KinematicCloud<ParcelType>::KinematicCloud"
-            "(const word&, const volScalarField&"
-            ", const volVectorField&, const volScalarField&, const "
-            "dimensionedVector&)"
-        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
-         << exit(FatalError);
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -265,7 +223,12 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         g_.value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (coupled_)
     {
@@ -276,160 +239,16 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = constProps_.rho0();
-
-    this->injection().prepareForNextTimeStep(time0_, time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - time0_
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - time0_
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                parcelTypeId_,
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                constProps_
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::injectParcel(ParcelType* p)
-{
-    addParticle(p);
-    nParcelsAdded_++;
-    nParcelsAddedTotal_++;
-    massInjected_ += p->mass()*p->nParticle();
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::postInjectCheck()
-{
-    if (nParcelsAdded_)
-    {
-        Pout<< "\n--> Cloud: " << this->name() << nl
-            << "    Added " << nParcelsAdded_
-            <<  " new parcels" << nl << endl;
-    }
-
-    // Reset parcel counters
-    nParcelsAdded_ = 0;
-
-    // Set time for start of next injection
-    time0_ = this->db().time().value();
-
-    // Increment number of injections
-    nInjections_++;
-}
-
-
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::info() const
 {
     Info<< "Cloud name: " << this->name() << nl
         << "    Parcels added during this run   = "
-        << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
+        << returnReduce(this->injection().nParcelsAddedTotal(), sumOp<label>())
+            << nl
         << "    Mass introduced during this run = "
-        << returnReduce(massInjected_, sumOp<scalar>()) << nl
+        << returnReduce(this->injection().massInjected(), sumOp<scalar>())
+            << nl
         << "    Current number of parcels       = "
         << returnReduce(this->size(), sumOp<label>()) << nl
         << "    Current mass in system          = "
@@ -445,7 +264,7 @@ void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
     (
         this->db().time().path()/"parcelPositions_"
       + this->name() + "_"
-      + name(this->nInjections_) + ".obj"
+      + name(this->injection().nInjections()) + ".obj"
     );
 
     forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 252d16a0f1..d606b38c3d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -85,19 +85,6 @@ class KinematicCloud
     public kinematicCloud
 {
 
-public:
-
-    // Enumerations
-
-        //- Parcel basis representation options
-        //  i.e constant number of particles OR constant mass per parcel
-        enum parcelBasis
-        {
-            pbNumber,
-            pbMass
-        };
-
-
 private:
 
     // Private data
@@ -126,22 +113,6 @@ private:
         //- Random number generator - used by some injection routines
         Random rndGen_;
 
-        //- Time at beginning of timestep
-        scalar time0_;
-
-
-        // Injection properties
-
-            //- Parcel basis
-            const word parcelBasisType_;
-            parcelBasis parcelBasis_;
-
-            //- Total mass to inject [kg]
-            scalar massTotal_;
-
-            //- Total mass injected to date [kg]
-            scalar massInjected_;
-
 
         // References to the carrier gas fields
 
@@ -161,9 +132,8 @@ private:
             const dimensionedVector& g_;
 
 
-        // Interpolation
-
-            dictionary interpolationSchemes_;
+        //- Interpolation schemes dictionary
+        dictionary interpolationSchemes_;
 
 
         // References to the cloud sub-models
@@ -190,17 +160,6 @@ private:
             autoPtr<vectorIntegrationScheme> UIntegrator_;
 
 
-        // Counters
-
-            //- Number of injections counter
-            label nInjections_;
-
-            //- Running counters describing parcels added during each
-            //  injection
-            label nParcelsAdded_;
-            label nParcelsAddedTotal_;
-
-
         // Sources
 
             //- Momentum
@@ -219,30 +178,6 @@ private:
         void operator=(const KinematicCloud&);
 
 
-protected:
-
-    // Protected member functions
-
-        //- Set the number of particles per parcel
-        scalar setNumberOfParticles
-        (
-            const label nParcels,
-            const scalar pDiameter,
-            const scalar pVolumeFraction,
-            const scalar pRho,
-            const scalar pVolume
-        );
-
-        //- Inject more parcels
-        void inject();
-
-        //- Inject parcel if it is valid - delete otherwise
-        void injectParcel(ParcelType* p);
-
-        //- Post-injection checks
-        void postInjectCheck();
-
-
 public:
 
     // Constructors
@@ -286,12 +221,6 @@ public:
             //- Return refernce to the random object
             inline Random& rndGen();
 
-            //- Return the start of injection interval time
-            inline scalar time0() const;
-
-            //- Return a reference to the mass of particles to introduce
-            inline scalar massTotal() const;
-
 
             // References to the carrier gas fields
 
@@ -380,15 +309,6 @@ public:
             void dumpParticlePositions() const;
 
 
-            // Counters
-
-                //- Return the number of injections
-                inline label nInjections() const;
-
-                //- Return the total number parcels added
-                inline label nParcelsAddedTotal() const;
-
-
             // Fields
 
                 //- Return the particle volume fraction field
@@ -402,6 +322,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index b354566c62..2b8cca12d1 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -155,41 +155,6 @@ Foam::KinematicCloud<ParcelType>::UIntegrator() const
 }
 
 
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nInjections() const
-{
-    return nInjections_;
-}
-
-
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nParcelsAddedTotal() const
-{
-    return nParcelsAddedTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::time0() const
-{
-    return time0_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massTotal() const
-{
-    return massTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInjected() const
-{
-    return massInjected_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInSystem() const
 {
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 93f3bb4f19..b574435cb8 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -29,6 +29,42 @@ License
 #include "MassTransferModel.H"
 #include "SurfaceReactionModel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ReactingCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        composition().YGas0(),
+        composition().YLiquid0(),
+        composition().YSolid0(),
+        composition().YMixture0(),
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -174,7 +210,12 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -185,123 +226,4 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::ReactingCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = this->constProps().rho0();
-
-    this->injection().prepareForNextTimeStep(this->time0(), time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - this->time0()
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - this->time0()
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = this->setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                this->parcelTypeId(),
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                composition().YGas0(),
-                composition().YLiquid0(),
-                composition().YSolid0(),
-                composition().YMixture0(),
-                this->constProps()
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->db().time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index 3e07ed56f8..0e3f97320a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -114,12 +114,6 @@ class ReactingCloud
         void operator=(const ReactingCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    void inject();
-
-
 public:
 
     //- Runtime type information
@@ -199,6 +193,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index eec7a00723..8a284024ff 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -30,6 +30,38 @@ License
 #include "interpolationCellPoint.H"
 #include "ThermoParcel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ThermoCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -167,7 +199,12 @@ void Foam::ThermoCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject(td);
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -178,16 +215,4 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoCloud<ParcelType>::inject
-(
-    TrackingData& td
-)
-{
-    // Injection is same as for KinematicCloud<ParcelType>
-    KinematicCloud<ParcelType>::inject(td);
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index 19c617a18e..84603dd4ff 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -116,13 +116,6 @@ class ThermoCloud
         void operator=(const ThermoCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    template<class TrackingData>
-    void inject(TrackingData& td);
-
-
 public:
 
     //- Runtime type information
@@ -173,7 +166,7 @@ public:
 
             // Modelling options
 
-                //- Radiation flag
+                 //- Radiation flag
                 inline bool radiation() const;
 
 
@@ -208,6 +201,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index 9b521d5ab7..9d54beba25 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -193,7 +193,7 @@ protected:
             scalar cp_;
 
 
-        // Call-based quantities
+        // Cell-based quantities
 
             //- Temperature [K]
             scalar Tc_;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 465949a120..3ef7e2a873 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -25,109 +25,22 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "InjectionModel.H"
+#include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::InjectionModel
-(
-    const dictionary& dict,
-    CloudType& owner,
-    const word& type
-)
-:   dict_(dict),
-    owner_(owner),
-    coeffDict_(dict.subDict(type + "Coeffs")),
-    SOI_(readScalar(coeffDict_.lookup("SOI"))),
-    volumeTotal_(0.0),
-    timeStep0_(0.0),
-    nParcels_(0),
-    volume_(0.0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::~InjectionModel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class CloudType>
-const CloudType& Foam::InjectionModel<CloudType>::owner() const
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-CloudType& Foam::InjectionModel<CloudType>::owner()
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::dict() const
-{
-    return dict_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::coeffDict() const
-{
-    return coeffDict_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::timeStart() const
-{
-    return SOI_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeTotal() const
-{
-    return volumeTotal_;
-}
-
-
-template<class CloudType>
-Foam::label Foam::InjectionModel<CloudType>::nParcels() const
-{
-    return nParcels_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volume() const
-{
-    return volume_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeFraction() const
-{
-    return volume_/volumeTotal_;
-}
-
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class CloudType>
 void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 (
     const scalar time0,
-    const scalar time1
+    const scalar time1,
+    label& nParcels,
+    scalar& volume
 )
 {
     // Initialise values
-    nParcels_ = 0;
-    volume_ = 0.0;
+    nParcels = 0;
+    volume = 0.0;
 
     // Return if not started injection event
     if (time1 < SOI_)
@@ -141,13 +54,13 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
     scalar t1 = time1 - SOI_;
 
     // Number of parcels to inject
-    nParcels_ = nParcelsToInject(t0, t1);
+    nParcels = nParcelsToInject(t0, t1);
 
     // Volume of parcels to inject
-    volume_ = volumeToInject(t0, t1);
+    volume = volumeToInject(t0, t1);
 
     // Hold previous time if no parcels, but non-zero volume fraction
-    if ((nParcels_ == 0) && (volume_ > 0.0))
+    if ((nParcels == 0) && (volume > 0.0))
     {
         // hold value of timeStep0_
     }
@@ -210,6 +123,203 @@ void Foam::InjectionModel<CloudType>::findInjectorCellAndPosition
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
+(
+    const label nParcels,
+    const scalar diameter,
+    const scalar volumeFraction,
+    const scalar rho,
+    const scalar volume
+)
+{
+    scalar nP = 0.0;
+    switch (parcelBasis_)
+    {
+        case pbMass:
+        {
+            nP = volumeFraction*massTotal_/nParcels
+               /(rho*mathematicalConstant::pi/6.0*pow3(diameter));
+            break;
+        }
+        case pbNumber:
+        {
+            nP = volumeFraction*massTotal_/(rho*volume);
+            break;
+        }
+        default:
+        {
+            nP = 0.0;
+            FatalErrorIn
+            (
+                "void Foam::InjectionModel<CloudType>::setNumberOfParticles"
+                "(const label, const scalar, const scalar, const scalar, "
+                "const scalar)"
+            )<< "Unknown parcelBasis type" << nl
+             << exit(FatalError);
+        }
+    }
+
+    return nP;
+}
+
+
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::postInjectCheck()
+{
+    if (nParcelsAdded_ > 0)
+    {
+        Pout<< "\n--> Cloud: " << owner_.name() << nl
+            << "    Added " << nParcelsAdded_
+            <<  " new parcels" << nl << endl;
+    }
+
+    // Increment total number of parcels added
+    nParcelsAddedTotal_ += nParcelsAdded_;
+
+    // Reset parcel counters
+    nParcelsAdded_ = 0;
+
+    // Update time for start of next injection
+    time0_ = owner_.db().time().value();
+
+    // Increment number of injections
+    nInjections_++;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::InjectionModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:   dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs")),
+    SOI_(readScalar(coeffDict_.lookup("SOI"))),
+    volumeTotal_(0.0),
+    massTotal_(dimensionedScalar(coeffDict_.lookup("massTotal")).value()),
+    massInjected_(0.0),
+    nInjections_(0),
+    nParcelsAdded_(0),
+    nParcelsAddedTotal_(0),
+    parcelBasisType_(coeffDict_.lookup("parcelBasisType")),
+    parcelBasis_(pbNumber),
+    time0_(owner.db().time().value()),
+    timeStep0_(0.0)
+{
+    if (parcelBasisType_ == "mass")
+    {
+        parcelBasis_ = pbMass;
+    }
+    else if (parcelBasisType_ == "number")
+    {
+        parcelBasis_ = pbNumber;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::InjectionModel<CloudType>::InjectionModel"
+            "(const dictionary&, CloudType&, const word&)"
+        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
+         << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::~InjectionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+template<class TrackData>
+void Foam::InjectionModel<CloudType>::inject(TrackData& td)
+{
+    const scalar time = owner_.db().time().value();
+    const scalar continuousDt = owner_.db().time().deltaT().value();
+
+    // Prepare for next time step
+    nParcelsAdded_ = 0;
+    label nParcels = 0;
+    scalar volume = 0.0;
+    prepareForNextTimeStep(time0_, time, nParcels, volume);
+
+    // Return if no parcels are required
+    if (nParcels == 0)
+    {
+        postInjectCheck();
+        return;
+    }
+
+    // Particle density given by constant properties
+    const scalar rho = td.constProps().rho0();
+
+    // Volume fraction to introduce during this timestep
+    const scalar volFraction = volumeFraction(volume);
+
+    // Duration of injection period during this timestep
+    const scalar deltaT = min
+    (
+        continuousDt,
+        min(time - SOI_, timeEnd() - time0_)
+    );
+
+    // Pad injection time if injection starts during this timestep
+    const scalar padTime = max(0.0, SOI_ - time0_);
+
+    // Introduce new parcels linearly with time
+    for (label iParcel=0; iParcel<nParcels; iParcel++)
+    {
+        // Calculate the pseudo time of injection for parcel 'iParcel'
+        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
+
+        // Determine injected parcel properties
+        vector pos = position(iParcel, timeInj, owner_.meshInfo());
+
+        // Diameter of parcels
+        scalar d = d0(iParcel, timeInj);
+
+        // Number of particles per parcel
+        scalar nP = setNumberOfParticles
+        (
+            nParcels,
+            d,
+            volFraction,
+            rho,
+            volume
+        );
+
+        // Velocity of parcels
+        vector U = velocity(iParcel, timeInj, owner_.meshInfo());
+
+        // Determine the injection cell
+        label cellI = -1;
+        findInjectorCellAndPosition(cellI, pos);
+
+        if (cellI >= 0)
+        {
+            scalar dt = time - timeInj;
+            td.cloud().addNewParcel(pos, cellI, d, U, nP, dt);
+
+            massInjected_ += nP*rho*mathematicalConstant::pi*pow3(d)/6.0;
+            nParcelsAdded_++;
+        }
+    }
+
+    postInjectCheck();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewInjectionModel.C"
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index 780aa2d776..791631fb99 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -55,6 +55,21 @@ template<class CloudType>
 class InjectionModel
 {
 
+public:
+
+    // Enumerations
+
+        //- Parcel basis representation options
+        //  i.e constant number of particles OR constant mass per parcel
+        enum parcelBasis
+        {
+            pbNumber,
+            pbMass
+        };
+
+
+private:
+
     // Private data
 
         //- The cloud dictionary
@@ -74,24 +89,47 @@ protected:
         // Global injection properties
 
             //- Start of injection [s]
-            scalar SOI_;
+            const scalar SOI_;
 
             //- Total volume of parcels to introduce [m^3]
             //  Initialised in the individual injection models
             scalar volumeTotal_;
 
+            //- Total mass to inject [kg]
+            const scalar massTotal_;
+
+            //- Total mass injected to date [kg]
+            scalar massInjected_;
+
+
+        // Counters
+
+            //- Number of injections counter
+            label nInjections_;
+
+            //- Running counter of parcels added during each injection
+            label nParcelsAdded_;
+
+            //- Running counter of total number of parcels added
+            label nParcelsAddedTotal_;
+
 
         // Injection properties per Lagrangian time step
 
+            // Parcel basis
+
+                //- Parcel basis name
+                const word parcelBasisType_;
+
+                //- Parcel basis enumeration
+                parcelBasis parcelBasis_;
+
+            //- Continuous phase time at start of injection time step [s]
+            scalar time0_;
+
             //- Time at start of injection time step [s]
             scalar timeStep0_;
 
-            //- Number of parcels to introduce []
-            label nParcels_;
-
-            //- Volume of parcels to introduce [m^3]
-            scalar volume_;
-
 
     // Protected member functions
 
@@ -110,6 +148,37 @@ protected:
         ) const = 0;
 
 
+        //- Determine properties for next time step/injection interval
+        void prepareForNextTimeStep
+        (
+            const scalar time0,
+            const scalar time1,
+            label& nParcels,
+            scalar& volume
+        );
+
+        //- Find the cell that contains the injector position
+        //  Will modify position slightly towards the owner cell centroid
+        virtual void findInjectorCellAndPosition
+        (
+            label& cellI,
+            vector& position
+        );
+
+        //- Set number of particles to inject given parcel properties
+        scalar setNumberOfParticles
+        (
+            const label nParcels,
+            const scalar diameter,
+            const scalar volumeFraction,
+            const scalar rho,
+            const scalar volume
+        );
+
+        //- Post injection checks
+        void postInjectCheck();
+
+
 public:
 
     //- Runtime type information
@@ -140,33 +209,31 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~InjectionModel();
+    //- Destructor
+    virtual ~InjectionModel();
 
 
-    // Selector
-
-        static autoPtr<InjectionModel<CloudType> > New
-        (
-            const dictionary& dict,
-            CloudType& owner
-        );
+    //- Selector
+    static autoPtr<InjectionModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
 
 
     // Access
 
+        //- Return the owner cloud dictionary
+        inline const dictionary& dict() const;
+
         //- Return const access the owner cloud object
-        const CloudType& owner() const;
+        inline const CloudType& owner() const;
 
         //- Return non-const access the owner cloud object for manipulation
-        CloudType& owner();
-
-        //- Return the dictionary
-        const dictionary& dict() const;
+        inline CloudType& owner();
 
         //- Return the coefficients dictionary
-        const dictionary& coeffDict() const;
+        inline const dictionary& coeffDict() const;
 
 
     // Member Functions
@@ -178,44 +245,41 @@ public:
         // Global information
 
             //- Return the start-of-injection time
-            scalar timeStart() const;
+            inline scalar timeStart() const;
 
             //- Return the total volume to be injected across the event
-            scalar volumeTotal() const;
+            inline scalar volumeTotal() const;
+
+            //- Return mass of particles to introduce
+            inline scalar massTotal() const;
+
+            //- Return mass of particles injected (cummulative)
+            inline scalar massInjected() const;
 
             //- Return the end-of-injection time
             virtual scalar timeEnd() const = 0;
 
+            // Counters
 
-        // Per Lagrangian time step properties
+                //- Return the number of injections
+                inline label nInjections() const;
 
-            //- Determine properties for next time step/injection interval
-            void prepareForNextTimeStep
-            (
-                const scalar time0,
-                const scalar time1
-            );
+                //- Return the total number parcels added
+                inline label nParcelsAddedTotal() const;
 
-            //- Return the number of parcels to introduce
-            label nParcels() const;
 
-            //- Return the volume of parcels to introduce
-            scalar volume() const;
+        // Per-injection event functions
+
+            //- Main injection loop
+            template<class TrackData>
+            void inject(TrackData& td);
 
             //- Return the volume fraction to introduce
-            scalar volumeFraction() const;
+            inline scalar volumeFraction(const scalar volume) const;
 
 
         // Injection geometry
 
-            //- Find the cell that contains the injector position
-            //  Will modify position slightly towards the owner cell centroid
-            virtual void findInjectorCellAndPosition
-            (
-                label& cellI,
-                vector& position
-            );
-
             //- Return the injection position
             virtual vector position
             (
@@ -241,6 +305,10 @@ public:
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "InjectionModelI.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModelI.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModelI.H
new file mode 100644
index 0000000000..fd4ef4d8c8
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModelI.H
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "InjectionModel.H"
+
+template<class CloudType>
+const Foam::dictionary& Foam::InjectionModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const CloudType& Foam::InjectionModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+CloudType& Foam::InjectionModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::InjectionModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::timeStart() const
+{
+    return SOI_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeTotal() const
+{
+    return volumeTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::massTotal() const
+{
+    return massTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::massInjected() const
+{
+    return massInjected_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::nInjections() const
+{
+    return nInjections_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::nParcelsAddedTotal() const
+{
+    return nParcelsAddedTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeFraction
+(
+    const scalar volume
+) const
+{
+    return volume/volumeTotal_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/foamCalcFunctions/Make/files b/src/postProcessing/foamCalcFunctions/Make/files
index 36dacdedb6..50dc36be0c 100644
--- a/src/postProcessing/foamCalcFunctions/Make/files
+++ b/src/postProcessing/foamCalcFunctions/Make/files
@@ -8,4 +8,6 @@ field/magGrad/magGrad.C
 field/div/div.C
 field/randomise/randomise.C
 
+basic/add/add.C
+
 LIB = $(FOAM_LIBBIN)/libfoamCalcFunctions
diff --git a/src/postProcessing/foamCalcFunctions/basic/add/add.C b/src/postProcessing/foamCalcFunctions/basic/add/add.C
new file mode 100644
index 0000000000..daf2aa6310
--- /dev/null
+++ b/src/postProcessing/foamCalcFunctions/basic/add/add.C
@@ -0,0 +1,287 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "add.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace calcTypes
+    {
+        defineTypeNameAndDebug(add, 0);
+        addToRunTimeSelectionTable(calcType, add, dictionary);
+    }
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::calcTypes::add::writeAddFields
+(
+    const Time& runTime,
+    const fvMesh& mesh,
+    const IOobject& baseFieldHeader
+)
+{
+    bool processed = false;
+
+    IOobject addFieldHeader
+    (
+        addFieldName_,
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    if (addFieldHeader.headerOk())
+    {
+        writeAddField<scalar>
+        (
+            baseFieldHeader,
+            addFieldHeader,
+            mesh,
+            processed
+        );
+        writeAddField<vector>
+        (
+            baseFieldHeader,
+            addFieldHeader,
+            mesh,
+            processed
+        );
+        writeAddField<sphericalTensor>
+        (
+            baseFieldHeader,
+            addFieldHeader,
+            mesh,
+            processed
+        );
+        writeAddField<symmTensor>
+        (
+            baseFieldHeader,
+            addFieldHeader,
+            mesh,
+            processed
+        );
+        writeAddField<tensor>
+        (
+            baseFieldHeader,
+            addFieldHeader,
+            mesh,
+            processed
+        );
+
+        if (!processed)
+        {
+            FatalError
+                << "Unable to process " << baseFieldName_
+                << " + " << addFieldName_ << nl
+                << "No call to add for fields of type "
+                << baseFieldHeader.headerClassName() << " + "
+                << addFieldHeader.headerClassName() << nl << nl
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        FatalErrorIn("calcTypes::add::writeAddFields()")
+            << "Unable to read add field: " << addFieldName_
+            << nl << exit(FatalError);
+    }
+}
+
+
+void Foam::calcTypes::add::writeAddValues
+(
+    const Time& runTime,
+    const fvMesh& mesh,
+    const IOobject& baseFieldHeader
+)
+{
+    bool processed = false;
+
+    writeAddValue<scalar>
+    (
+        baseFieldHeader,
+        addValueStr_,
+        mesh,
+        processed
+    );
+    writeAddValue<vector>
+    (
+        baseFieldHeader,
+        addValueStr_,
+        mesh,
+        processed
+    );
+    writeAddValue<sphericalTensor>
+    (
+        baseFieldHeader,
+        addValueStr_,
+        mesh,
+        processed
+    );
+    writeAddValue<symmTensor>
+    (
+        baseFieldHeader,
+        addValueStr_,
+        mesh,
+        processed
+    );
+    writeAddValue<tensor>
+    (
+        baseFieldHeader,
+        addValueStr_,
+        mesh,
+        processed
+    );
+
+    if (!processed)
+    {
+        FatalErrorIn("calcTypes::add::writeAddValues()")
+            << "Unable to process " << baseFieldName_
+            << " + " << addValueStr_ << nl
+            << "No call to add for fields of type "
+            << baseFieldHeader.headerClassName() << nl << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::calcTypes::add::add()
+:
+    calcType(),
+    baseFieldName_(""),
+    calcType_(FIELD),
+    addFieldName_(""),
+    addValueStr_(""),
+    resultName_("")
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::calcTypes::add::~add()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::calcTypes::add::init()
+{
+    argList::validArgs.append("add");
+    argList::validArgs.append("baseField");
+    argList::validOptions.insert("field", "fieldName");
+    argList::validOptions.insert("value", "valueString");
+    argList::validOptions.insert("resultName", "fieldName");
+}
+
+
+void Foam::calcTypes::add::preCalc
+(
+    const argList& args,
+    const Time& runTime,
+    const fvMesh& mesh
+)
+{
+    baseFieldName_ = args.additionalArgs()[1];
+
+    if (args.options().found("field"))
+    {
+        addFieldName_ = args.options()["field"];
+        calcType_ = FIELD;
+    }
+    else if (args.options().found("value"))
+    {
+        addValueStr_ = args.options()["value"];
+        calcType_ = VALUE;
+    }
+    else
+    {
+        FatalErrorIn("calcTypes::add::preCalc")
+            << "add requires either -field or -value option"
+            << nl << exit(FatalError);
+    }
+
+    if (args.options().found("resultName"))
+    {
+        resultName_ = args.options()["resultName"];
+    }
+}
+
+
+void Foam::calcTypes::add::calc
+(
+    const argList& args,
+    const Time& runTime,
+    const fvMesh& mesh
+)
+{
+    IOobject baseFieldHeader
+    (
+        baseFieldName_,
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    if (baseFieldHeader.headerOk())
+    {
+        switch (calcType_)
+        {
+            case FIELD:
+            {
+                writeAddFields(runTime, mesh, baseFieldHeader);
+                break;
+            }
+            case VALUE:
+            {
+                writeAddValues(runTime, mesh, baseFieldHeader);
+                break;
+            }
+            default:
+            {
+                FatalErrorIn("calcTypes::add::calc")
+                    << "unknown calcType " << calcType_ << nl
+                    << abort(FatalError);
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("calcTypes::add::calc")
+            << "Unable to read base field: " << baseFieldName_
+            << nl << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //
+
diff --git a/src/postProcessing/foamCalcFunctions/basic/add/add.H b/src/postProcessing/foamCalcFunctions/basic/add/add.H
new file mode 100644
index 0000000000..ae5f517618
--- /dev/null
+++ b/src/postProcessing/foamCalcFunctions/basic/add/add.H
@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::calcTypes::add
+
+Description
+    Adds and field or value to base field.
+
+    New field name specified by -resultName option, or automatically as:
+        <baseFieldName>_plus_<addFieldName>
+        <baseFieldName>_plus_value
+
+    Example usage:
+        add p -value 100000 -resultName pAbs
+        add U -field U0
+
+SourceFiles
+    add.C
+    writeAddField.C
+    writeAddValue.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef add_H
+#define add_H
+
+#include "calcType.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace calcTypes
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class add Declaration
+\*---------------------------------------------------------------------------*/
+
+class add
+:
+    public calcType
+{
+public:
+
+    enum calcTypes
+    {
+        FIELD,
+        VALUE
+    };
+
+
+private:
+
+    // Private data
+
+        //- Name of base field (to add to)
+        word baseFieldName_;
+
+        //- Calc type as given by enumerations above
+        calcTypes calcType_;
+
+        //- Name of field to add
+        word addFieldName_;
+
+        //- String representation of value to add
+        string addValueStr_;
+
+        //- Name of result field
+        word resultName_;
+
+
+    // Private Member Functions
+
+        // Output
+
+            //- Calc and output field additions
+            void writeAddFields
+            (
+                const Time& runTime,
+                const fvMesh& mesh,
+                const IOobject& baseFieldHeader
+            );
+
+            //- Calc and output field and value additions
+            void writeAddValues
+            (
+                const Time& runTime,
+                const fvMesh& mesh,
+                const IOobject& baseFieldHeader
+            );
+
+
+        //- Disallow default bitwise copy construct
+        add(const add&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const add&);
+
+
+protected:
+
+    // Member Functions
+
+        // Calculation routines
+
+            //- Initialise - typically setting static variables,
+            //  e.g. command line arguments
+            virtual void init();
+
+            //- Pre-time loop calculations
+            virtual void preCalc
+            (
+                const argList& args,
+                const Time& runTime,
+                const fvMesh& mesh
+            );
+
+            //- Time loop calculations
+            virtual void calc
+            (
+                const argList& args,
+                const Time& runTime,
+                const fvMesh& mesh
+            );
+
+
+        // I-O
+
+            //- Write add field
+            template<class Type>
+            void writeAddField
+            (
+                const IOobject& baseHeader,
+                const IOobject& addHeader,
+                const fvMesh& mesh,
+                bool& processed
+            );
+
+            //- Write add value
+            template<class Type>
+            void writeAddValue
+            (
+                const IOobject& baseHeader,
+                const string& valueStr,
+                const fvMesh& mesh,
+                bool& processed
+            );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("add");
+
+
+    // Constructors
+
+        //- Construct null
+        add();
+
+
+    // Destructor
+
+        virtual ~add();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace calcTypes
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "writeAddField.C"
+#   include "writeAddValue.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/foamCalcFunctions/basic/add/writeAddField.C b/src/postProcessing/foamCalcFunctions/basic/add/writeAddField.C
new file mode 100644
index 0000000000..78b51146a9
--- /dev/null
+++ b/src/postProcessing/foamCalcFunctions/basic/add/writeAddField.C
@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+void Foam::calcTypes::add::writeAddField
+(
+    const IOobject& baseHeader,
+    const IOobject& addHeader,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    if (resultName_ == "")
+    {
+        resultName_ = baseHeader.name() + "_plus_" + addHeader.name();
+    }
+
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if
+    (
+        baseHeader.headerClassName() == fieldType::typeName
+     && baseHeader.headerClassName() == addHeader.headerClassName()
+    )
+    {
+        Info<< "    Reading " << baseHeader.name() << endl;
+        fieldType baseField(baseHeader, mesh);
+
+        Info<< "    Reading " << addHeader.name() << endl;
+        fieldType addField(addHeader, mesh);
+
+        if (baseField.dimensions() == addField.dimensions())
+        {
+            Info<< "    Calculating " << resultName_ << endl;
+
+            fieldType newField
+            (
+                IOobject
+                (
+                    resultName_,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                baseField + addField
+            );
+            newField.write();
+        }
+        else
+        {
+            Info<< "    Cannot calculate " << resultName_ << nl
+                << "    - inconsistent dimensions: "
+                << baseField.dimensions() << " - " << addField.dimensions()
+                << endl;
+        }
+
+        processed = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/postProcessing/foamCalcFunctions/basic/add/writeAddValue.C b/src/postProcessing/foamCalcFunctions/basic/add/writeAddValue.C
new file mode 100644
index 0000000000..c62720549e
--- /dev/null
+++ b/src/postProcessing/foamCalcFunctions/basic/add/writeAddValue.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+void Foam::calcTypes::add::writeAddValue
+(
+    const IOobject& baseHeader,
+    const string& valueStr,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (baseHeader.headerClassName() == fieldType::typeName)
+    {
+        if (resultName_ == "")
+        {
+            resultName_ = baseHeader.name() + "_plus_value";
+        }
+
+        Type value;
+        IStringStream(valueStr)() >> value;
+
+        Info<< "    Reading " << baseHeader.name() << endl;
+            fieldType baseField(baseHeader, mesh);
+
+        fieldType newField
+        (
+            IOobject
+            (
+                resultName_,
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            baseField
+        );
+
+        Info<< "    Calculating " << resultName_ << endl;
+        newField == baseField
+            + dimensioned<Type>("value", baseField.dimensions(), value);
+        newField.write();
+
+        processed = true;
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.C b/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.C
index e30bb17f99..92b3f96366 100644
--- a/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.C
+++ b/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.C
@@ -170,6 +170,7 @@ Foam::fieldAverage::fieldAverage
     name_(name),
     obr_(obr),
     active_(true),
+    cleanRestart_(dict.lookupOrDefault<Switch>("cleanRestart", false)),
     faItems_(dict.lookup("fields")),
     meanScalarFields_(faItems_.size()),
     meanVectorFields_(faItems_.size()),
@@ -326,29 +327,44 @@ void Foam::fieldAverage::writeAveragingProperties() const
 
 void Foam::fieldAverage::readAveragingProperties()
 {
-    IFstream propsFile
-    (
-        obr_.time().path()/obr_.time().timeName()
-       /"uniform"/"fieldAveragingProperties"
-    );
-
-    if (!propsFile.good())
+    if (cleanRestart_)
     {
-        return;
+        Info<< "fieldAverage: starting averaging at time "
+            << obr_.time().timeName() << nl << endl;
     }
-
-    dictionary propsDict(dictionary::null, propsFile);
-
-    forAll(faItems_, i)
+    else
     {
-        const word& fieldName = faItems_[i].fieldName();
-        if (propsDict.found(fieldName))
-        {
-            dictionary fieldDict(propsDict.subDict(fieldName));
+        IFstream propsFile
+        (
+            obr_.time().path()/obr_.time().timeName()
+            /"uniform"/"fieldAveragingProperties"
+        );
 
-            totalIter_[i] = readLabel(fieldDict.lookup("totalIter"));
-            totalTime_[i] = readScalar(fieldDict.lookup("totalTime"));
+        if (!propsFile.good())
+        {
+            Info<< "fieldAverage: starting averaging at time "
+                << obr_.time().timeName() << nl << endl;
+            return;
         }
+
+        dictionary propsDict(dictionary::null, propsFile);
+
+        Info<< "fieldAverage: restarting averaging for fields:" << endl;
+        forAll(faItems_, i)
+        {
+            const word& fieldName = faItems_[i].fieldName();
+            if (propsDict.found(fieldName))
+            {
+                dictionary fieldDict(propsDict.subDict(fieldName));
+
+                totalIter_[i] = readLabel(fieldDict.lookup("totalIter"));
+                totalTime_[i] = readScalar(fieldDict.lookup("totalTime"));
+                Info<< "    " << fieldName
+                    << " iters = " << totalIter_[i]
+                    << " time = " << totalTime_[i] << endl;
+            }
+        }
+        Info<< endl;
     }
 }
 
diff --git a/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.H b/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.H
index e2437c2bde..7d9c4db00f 100644
--- a/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.H
+++ b/src/postProcessing/functionObjects/fieldAverage/fieldAverage/fieldAverage.H
@@ -37,6 +37,10 @@ Description
         // Where to load it from (if not already in solver)
         functionObjectLibs ("libfieldAverage.so");
 
+        // Whether to perform a clean restart, or start from previous
+        // averaging info if available
+        cleanRestart true;
+
         // Fields to be probed. runTime modifiable!
         fields
         (
@@ -80,6 +84,7 @@ SourceFiles
 
 #include "volFieldsFwd.H"
 #include "pointFieldFwd.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -114,6 +119,9 @@ protected:
         //- On/off switch
         bool active_;
 
+        //- Clean restart flag
+        Switch cleanRestart_;
+
         //- List of field average items, describing waht averages to be
         //  calculated and output
         List<fieldAverageItem> faItems_;
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options
index 4c9cff9d95..cffe6f344a 100644
--- a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options
@@ -6,7 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
 
 EXE_LIBS = \
     -llagrangian \
@@ -18,4 +18,5 @@ EXE_LIBS = \
     -lcombustionThermophysicalModels \
     -lspecie \
     -lradiation \
-    -lcompressibleRASModels
+    -lcompressibleRASModels \
+    -lcompressibleLESModels
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H
index a824308f13..c22eeb8cc3 100644
--- a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H
@@ -41,18 +41,11 @@
 
 
     Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
+    autoPtr<compressible::turbulenceModel> turbulence
     (
-        compressible::RASModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo()
-        )
+        compressible::turbulenceModel::New(rho, U, phi, thermo())
     );
 
-
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C
index 999140140b..e7f53341ab 100644
--- a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C
@@ -32,7 +32,7 @@ Description
 
 #include "fvCFD.H"
 #include "basicThermo.H"
-#include "compressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 
 #include "basicThermoCloud.H"
 #include "basicKinematicCloud.H"
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
index 6ae6ae3a91..d883b02f39 100644
--- a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
@@ -20,12 +20,6 @@ DragModel                                SphereDrag;
 DispersionModel                          StochasticDispersionRAS;
 WallInteractionModel                     StandardWallInteraction;
 
-// Parcel basis type
-parcelBasisType                          mass;
-
-// Total mass to inject
-massTotal  massTotal [ 1  0  0  0  0]    2.0e-4;
-
 // Minimum particle mass
 minParticleMass      minParticleMass     [ 1  0  0  0  0]     1.0e-15;
 
@@ -52,6 +46,12 @@ integrationSchemes
 
 ManualInjectionCoeffs
 {
+    // Parcel basis type
+    parcelBasisType                          mass;
+
+    // Total mass to inject
+    massTotal  massTotal [ 1  0  0  0  0]    2.0e-4;
+
     SOI                                  0.0;
     positionsFile                        kinematicCloud1Positions;
     U0                                   (0 0 0);
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
index 0ad5e425cd..795aeda2cd 100644
--- a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
@@ -23,12 +23,6 @@ HeatTransferModel                        RanzMarshall;
 
 radiation                                off;
 
-// Parcel basis type
-parcelBasisType                          mass;
-
-// Total mass to inject
-massTotal  massTotal [ 1  0  0  0  0]    1e-4;
-
 // Minimum particle mass
 minParticleMass      minParticleMass     [ 1  0  0  0  0]     1.0e-15;
 
@@ -62,6 +56,12 @@ integrationSchemes
 
 ManualInjectionCoeffs
 {
+    // Total mass to inject
+    massTotal  massTotal [ 1  0  0  0  0]    1e-4;
+
+    // Parcel basis type
+    parcelBasisType                          mass;
+
     SOI                                  0.0;
     positionsFile                        thermoCloud1Positions;
     U0                                   (0 0 0);