diff --git a/applications/solvers/lagrangian/DPMFoam/Allwclean b/applications/solvers/lagrangian/DPMFoam/Allwclean
index 390c1235d2..25d0b2955f 100755
--- a/applications/solvers/lagrangian/DPMFoam/Allwclean
+++ b/applications/solvers/lagrangian/DPMFoam/Allwclean
@@ -5,3 +5,4 @@ set -x
 
 wclean DPMTurbulenceModels
 wclean
+wclean MPPICFoam
diff --git a/applications/solvers/lagrangian/DPMFoam/Allwmake b/applications/solvers/lagrangian/DPMFoam/Allwmake
index 7cff009a34..6308a7052b 100755
--- a/applications/solvers/lagrangian/DPMFoam/Allwmake
+++ b/applications/solvers/lagrangian/DPMFoam/Allwmake
@@ -5,3 +5,4 @@ set -x
 
 wmake DPMTurbulenceModels
 wmake
+wmake MPPICFoam
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 6d4e26a441..9cadf55117 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -19,6 +19,8 @@ KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
 $(KINEMATICPARCEL)/defineBasicKinematicParcel.C
 $(KINEMATICPARCEL)/makeBasicKinematicParcelSubmodels.C
 
+
+/* kinematic colliding parcel sub-models */
 KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/basicKinematicCollidingParcel
 $(KINEMATICCOLLIDINGPARCEL)/defineBasicKinematicCollidingParcel.C
 $(KINEMATICCOLLIDINGPARCEL)/makeBasicKinematicCollidingParcelSubmodels.C
@@ -42,6 +44,12 @@ $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
 $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 
 
+/* kinematic MPPIC parcel sub-models */
+KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel
+$(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
+$(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
+
+
 /* bolt-on models */
 RADIATION=submodels/addOns/radiation
 $(RADIATION)/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
@@ -71,6 +79,24 @@ $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphase
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIO.C
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIOList.C
 
+MPPICPARTICLESTRESS=submodels/MPPIC/ParticleStressModels
+$(MPPICPARTICLESTRESS)/ParticleStressModel/ParticleStressModel.C
+$(MPPICPARTICLESTRESS)/HarrisCrighton/HarrisCrighton.C
+$(MPPICPARTICLESTRESS)/Lun/Lun.C
+$(MPPICPARTICLESTRESS)/exponential/exponential.C
+
+MPPICCORRECTIONLIMITING=submodels/MPPIC/CorrectionLimitingMethods
+$(MPPICCORRECTIONLIMITING)/CorrectionLimitingMethod/CorrectionLimitingMethod.C
+$(MPPICCORRECTIONLIMITING)/noCorrectionLimiting/noCorrectionLimiting.C
+$(MPPICCORRECTIONLIMITING)/absolute/absolute.C
+$(MPPICCORRECTIONLIMITING)/relative/relative.C
+
+MPPICTIMESCALE=submodels/MPPIC/TimeScaleModels
+$(MPPICTIMESCALE)/TimeScaleModel/TimeScaleModel.C
+$(MPPICTIMESCALE)/equilibrium/equilibrium.C
+$(MPPICTIMESCALE)/nonEquilibrium/nonEquilibrium.C
+$(MPPICTIMESCALE)/isotropic/isotropic.C
+
 
 /* integration schemes */
 IntegrationScheme/makeIntegrationSchemes.C
@@ -81,8 +107,13 @@ phaseProperties/phaseProperties/phaseProperties.C
 phaseProperties/phaseProperties/phasePropertiesIO.C
 phaseProperties/phasePropertiesList/phasePropertiesList.C
 
-/* Additional helper classes */
+
+/* additional helper classes */
 clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
 
 
+/* averaging methods */
+submodels/MPPIC/AveragingMethods/makeAveragingMethods.C
+
+
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
new file mode 100644
index 0000000000..f22cac7b0f
--- /dev/null
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@@ -0,0 +1,299 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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"
+#include "PackingModel.H"
+#include "ParticleStressModel.H"
+#include "DampingModel.H"
+#include "IsotropyModel.H"
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::setModels()
+{
+    packingModel_.reset
+    (
+        PackingModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+    dampingModel_.reset
+    (
+        DampingModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+    isotropyModel_.reset
+    (
+        IsotropyModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    const word& cloudName,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const volScalarField& mu,
+    const dimensionedVector& g,
+    bool readFields
+)
+:
+    CloudType(cloudName, rho, U, mu, g, false),
+    packingModel_(NULL),
+    dampingModel_(NULL),
+    isotropyModel_(NULL)
+{
+    if (this->solution().steadyState())
+    {
+        FatalErrorIn
+        (
+            "Foam::MPPICCloud<CloudType>::MPPICCloud"
+            "("
+                "const word&, "
+                "const volScalarField&, "
+                "const volVectorField&, "
+                "const volScalarField&, "
+                "const dimensionedVector&, "
+                "bool"
+            ")"
+        )   << "MPPIC modelling not available for steady state calculations"
+            << exit(FatalError);
+    }
+
+    if (this->solution().active())
+    {
+        setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this);
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    MPPICCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    packingModel_(c.packingModel_->clone()),
+    dampingModel_(c.dampingModel_->clone()),
+    isotropyModel_(c.isotropyModel_->clone())
+{}
+
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const MPPICCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    packingModel_(NULL),
+    dampingModel_(NULL),
+    isotropyModel_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::~MPPICCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<MPPICCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::restoreState()
+{
+    this->cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::template
+            TrackingData<MPPICCloud<CloudType> > td(*this);
+
+        this->solve(td);
+    }
+}
+
+
+template<class CloudType>
+template<class TrackData>
+void Foam::MPPICCloud<CloudType>::motion(TrackData& td)
+{
+    // Kinematic
+    // ~~~~~~~~~
+
+    // force calculation and tracking
+    td.part() = TrackData::tpLinearTrack;
+    CloudType::move(td, this->db().time().deltaTValue());
+
+
+    // Preliminary
+    // ~~~~~~~~~~~
+
+    // switch forces off so they are not applied in corrector steps
+    this->forces().setCalcNonCoupled(false);
+    this->forces().setCalcCoupled(false);
+
+
+    // Damping
+    // ~~~~~~~
+
+    if (dampingModel_->active())
+    {
+        // update averages
+        td.updateAverages(*this);
+
+        // memory allocation and eulerian calculations
+        dampingModel_->cacheFields(true);
+
+        // calculate the damping velocity corrections without moving the parcels
+        td.part() = TrackData::tpDampingNoTrack;
+        CloudType::move(td, this->db().time().deltaTValue());
+
+        // correct the parcel positions and velocities
+        td.part() = TrackData::tpCorrectTrack;
+        CloudType::move(td, this->db().time().deltaTValue());
+
+        // finalise and free memory
+        dampingModel_->cacheFields(false);
+    }
+
+
+    // Packing
+    // ~~~~~~~
+
+    if (packingModel_->active())
+    {
+        // same procedure as for damping
+        td.updateAverages(*this);
+        packingModel_->cacheFields(true);
+        td.part() = TrackData::tpPackingNoTrack;
+        CloudType::move(td, this->db().time().deltaTValue());
+        td.part() = TrackData::tpCorrectTrack;
+        CloudType::move(td, this->db().time().deltaTValue());
+        packingModel_->cacheFields(false);
+    }
+
+
+    // Isotropy
+    // ~~~~~~~~
+
+    if (isotropyModel_->active())
+    {
+        // update averages
+        td.updateAverages(*this);
+
+        // apply isotropy model
+        isotropyModel_->calculate();
+    }
+
+
+    // Final
+    // ~~~~~
+
+    // update cell occupancy
+    this->updateCellOccupancy();
+
+    // switch forces back on
+    this->forces().setCalcNonCoupled(true);
+    this->forces().setCalcCoupled(this->solution().coupled());
+}
+
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::info()
+{
+    CloudType::info();
+
+    tmp<volScalarField> alpha = this->theta();
+
+    Info<< "    Min cell volume fraction        = "
+        << gMin(alpha().internalField()) << endl;
+    Info<< "    Max cell volume fraction        = "
+        << gMax(alpha().internalField()) << endl;
+
+    label nOutside = 0;
+
+    forAllIter(typename CloudType, *this, iter)
+    {
+        typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), this->mesh());
+        nOutside += !tetIs.tet(this->mesh()).inside(p.position());
+    }
+
+    reduce(nOutside, plusOp<label>());
+
+    if (nOutside > 0)
+    {
+        Info<< "    Number of parcels outside tets  = " << nOutside << endl;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
new file mode 100644
index 0000000000..2b142fae63
--- /dev/null
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
@@ -0,0 +1,249 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::MPPICCloud
+
+Description
+    Adds collisions to kinematic clouds
+
+SourceFiles
+    MPPICCloudI.H
+    MPPICCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MPPICCloud_H
+#define MPPICCloud_H
+
+#include "particle.H"
+#include "Cloud.H"
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "fvMesh.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+template<class CloudType>
+class PackingModel;
+
+template<class CloudType>
+class DampingModel;
+
+template<class CloudType>
+class IsotropyModel;
+
+/*---------------------------------------------------------------------------*\
+                       Class MPPICCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MPPICCloud
+:
+    public CloudType
+{
+public:
+
+    // Public typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::parcelType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef MPPICCloud<CloudType> MPPICCloudType;
+
+
+private:
+
+    // Private data
+
+        //- Cloud copy pointer
+        autoPtr<MPPICCloud<CloudType> > cloudCopyPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        MPPICCloud(const MPPICCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const MPPICCloud&);
+
+
+protected:
+
+    // Protected data
+
+        // References to the cloud sub-models
+
+            //- Packing model
+            autoPtr<PackingModel<MPPICCloud<CloudType> > > packingModel_;
+
+            //- Damping model
+            autoPtr<DampingModel<MPPICCloud<CloudType> > >
+                dampingModel_;
+
+            //- Exchange model
+            autoPtr<IsotropyModel<MPPICCloud<CloudType> > >
+                isotropyModel_;
+
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        MPPICCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            const dimensionedVector& g,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        MPPICCloud
+        (
+            MPPICCloud<CloudType>& c,
+            const word& name
+        );
+
+        //- Copy constructor with new name - creates bare cloud
+        MPPICCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const MPPICCloud<CloudType>& c
+        );
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new MPPICCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new MPPICCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~MPPICCloud();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the cloud copy
+            inline const MPPICCloud& cloudCopy() const;
+
+            //- Return const access to the packing model
+            inline const PackingModel<MPPICCloud<CloudType> >&
+                packingModel() const;
+
+            //- Return a reference to the packing model
+            inline PackingModel<MPPICCloud<CloudType> >& packingModel();
+
+            //- Return condt access to the damping model
+            inline const DampingModel<MPPICCloud<CloudType> >&
+                dampingModel() const;
+
+            //- Return a reference to the damping model
+            inline DampingModel<MPPICCloud<CloudType> >& dampingModel();
+
+            //- Return condt access to the isotropy model
+            inline const IsotropyModel<MPPICCloud<CloudType> >&
+                isotropyModel() const;
+
+            //- Return a reference to the isotropy model
+            inline IsotropyModel<MPPICCloud<CloudType> >& isotropyModel();
+
+
+        // Cloud evolution functions
+
+            //- Store the current cloud state
+            void storeState();
+
+            //- Reset the current cloud to the previously stored state
+            void restoreState();
+
+            //- Evolve the cloud
+            void evolve();
+
+            //- Particle motion
+            template<class TrackData>
+            void motion(TrackData& td);
+
+
+        //- I-O
+
+            //- Print cloud information
+            void info();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MPPICCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MPPICCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudI.H b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudI.H
new file mode 100644
index 0000000000..5211f87002
--- /dev/null
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudI.H
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::MPPICCloud<CloudType>&
+Foam::MPPICCloud<CloudType>::cloudCopy() const
+{
+    return cloudCopyPtr_();
+}
+
+
+template<class CloudType>
+inline const Foam::PackingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::packingModel() const
+{
+    return packingModel_();
+}
+
+
+template<class CloudType>
+inline Foam::PackingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::packingModel()
+{
+    return packingModel_();
+}
+
+
+template<class CloudType>
+inline const Foam::DampingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::dampingModel() const
+{
+    return dampingModel_();
+}
+
+
+template<class CloudType>
+inline Foam::DampingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::dampingModel()
+{
+    return dampingModel_();
+}
+
+
+template<class CloudType>
+inline const Foam::IsotropyModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::isotropyModel() const
+{
+    return isotropyModel_();
+}
+
+
+template<class CloudType>
+inline Foam::IsotropyModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::isotropyModel()
+{
+    return isotropyModel_();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/basicKinematicMPPICCloud/basicKinematicMPPICCloud.H b/src/lagrangian/intermediate/clouds/derived/basicKinematicMPPICCloud/basicKinematicMPPICCloud.H
new file mode 100644
index 0000000000..e239b0606f
--- /dev/null
+++ b/src/lagrangian/intermediate/clouds/derived/basicKinematicMPPICCloud/basicKinematicMPPICCloud.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::basicKinematicMPPICCloud
+
+Description
+    Cloud class to introduce kinematic MPPIC parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicKinematicMPPICCloud_H
+#define basicKinematicMPPICCloud_H
+
+#include "Cloud.H"
+#include "KinematicCloud.H"
+#include "MPPICCloud.H"
+#include "basicKinematicMPPICParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef MPPICCloud
+    <
+        KinematicCloud
+        <
+            Cloud
+            <
+                basicKinematicMPPICParcel
+            >
+        >
+    > basicKinematicMPPICCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.C b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.C
new file mode 100644
index 0000000000..3844d22f67
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.C
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "MPPICParcel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const MPPICParcel<ParcelType>& p
+)
+:
+    ParcelType(p),
+    UCorrect_(p.UCorrect_)
+{}
+
+
+template<class ParcelType>
+Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const MPPICParcel<ParcelType>& p,
+    const polyMesh& mesh
+)
+:
+    ParcelType(p, mesh),
+    UCorrect_(p.UCorrect_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackData>
+bool Foam::MPPICParcel<ParcelType>::move
+(
+    TrackData& td,
+    const scalar trackTime
+)
+{
+    typename TrackData::cloudType::parcelType& p =
+        static_cast<typename TrackData::cloudType::parcelType&>(*this);
+
+    switch (td.part())
+    {
+        case TrackData::tpLinearTrack:
+        {
+            ParcelType::move(td, trackTime);
+
+            break;
+        }
+        case TrackData::tpDampingNoTrack:
+        {
+            p.UCorrect() =
+                td.cloud().dampingModel().velocityCorrection(p, trackTime);
+
+            break;
+        }
+        case TrackData::tpPackingNoTrack:
+        {
+            p.UCorrect() =
+                td.cloud().packingModel().velocityCorrection(p, trackTime);
+
+            break;
+        }
+        case TrackData::tpCorrectTrack:
+        {
+            vector U = p.U();
+
+            scalar f = p.stepFraction();
+
+            p.U() = (1.0 - f)*p.UCorrect();
+
+            ParcelType::move(td, trackTime);
+
+            p.U() = U + (p.stepFraction() - f)*p.UCorrect();
+
+            break;
+        }
+    }
+
+    return td.keepParticle;
+}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+#include "MPPICParcelIO.C"
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
new file mode 100644
index 0000000000..d0917a06ef
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::MPPICParcel
+
+Description
+    Wrapper around kinematic parcel types to add collision modelling
+
+SourceFiles
+    MPPICParcelI.H
+    MPPICParcel.C
+    MPPICParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MPPICParcel_H
+#define MPPICParcel_H
+
+#include "particle.H"
+#include "labelFieldIOField.H"
+#include "vectorFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of clases
+
+template<class ParcelType>
+class MPPICParcel;
+
+template<class Type>
+class AveragingMethod;
+
+// Forward declaration of friend functions
+
+template<class ParcelType>
+Ostream& operator<<
+(
+    Ostream&,
+    const MPPICParcel<ParcelType>&
+);
+
+/*---------------------------------------------------------------------------*\
+                       Class MPPICParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class MPPICParcel
+:
+    public ParcelType
+{
+public:
+
+    template<class CloudType>
+    class TrackingData
+    :
+        public ParcelType::template TrackingData<CloudType>
+    {
+
+    public:
+
+        enum trackPart
+        {
+            tpLinearTrack,
+            tpDampingNoTrack,
+            tpPackingNoTrack,
+            tpCorrectTrack,
+        };
+
+
+    private:
+
+        // Private data
+
+            // MPPIC Averages
+
+                //- Volume average
+                autoPtr<AveragingMethod<scalar> > volumeAverage_;
+
+                //- Radius average [ volume^(1/3) ]
+                autoPtr<AveragingMethod<scalar> > radiusAverage_;
+
+                //- Density average
+                autoPtr<AveragingMethod<scalar> > rhoAverage_;
+
+                //- Velocity average
+                autoPtr<AveragingMethod<vector> > uAverage_;
+
+                //- Magnitude velocity sqyuared average
+                autoPtr<AveragingMethod<scalar> > uSqrAverage_;
+
+                //- Frequency average
+                autoPtr<AveragingMethod<scalar> > frequencyAverage_;
+
+                //- Mass average
+                autoPtr<AveragingMethod<scalar> > massAverage_;
+
+
+            //- Label specifying the current part of the tracking process
+            trackPart part_;
+
+
+    public:
+
+        //- Constructors
+
+            //- Construct from components
+            inline TrackingData
+            (
+                CloudType& cloud,
+                trackPart part = tpLinearTrack
+            );
+
+
+        //- Update the MPPIC averages
+        inline void updateAverages(CloudType& cloud);
+
+
+        //- Access
+
+            //- Const access to the tracking part label
+            inline trackPart part() const;
+
+            //- Non const access to the tracking part label
+            inline trackPart& part();
+    };
+
+
+protected:
+
+    // Protected data
+
+        //- Velocity correction due to collisions [m/s]
+        vector UCorrect_;
+
+
+public:
+
+    // Static data members
+
+        //- Runtime type information
+        TypeName("MPPICParcel");
+
+        //- String representation of properties
+        AddToPropertyList
+        (
+            ParcelType,
+            "(UCorrectx UCorrecty UCorrectz)"
+        );
+
+
+    // Constructors
+
+        //- Construct from owner, position, and cloud owner
+        //  Other properties initialised as null
+        inline MPPICParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
+        );
+
+        //- Construct from components
+        inline MPPICParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const label typeId,
+            const scalar nParticle0,
+            const scalar d0,
+            const scalar dTarget0,
+            const vector& U0,
+            const vector& UCorrect0,
+            const typename ParcelType::constantProperties& constProps
+        );
+
+        //- Construct from Istream
+        MPPICParcel
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct as a copy
+        MPPICParcel(const MPPICParcel& p);
+
+        //- Construct as a copy
+        MPPICParcel(const MPPICParcel& p, const polyMesh& mesh);
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new MPPICParcel(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>(new MPPICParcel(*this, mesh));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<MPPICParcel<ParcelType> > operator()(Istream& is) const
+            {
+                return autoPtr<MPPICParcel<ParcelType> >
+                (
+                    new MPPICParcel<ParcelType>(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to correction velocity
+            inline const vector& UCorrect() const;
+
+            //- Return access to correction velocity
+            inline vector& UCorrect();
+
+
+        // Tracking
+
+            //- Move the parcel
+            template<class TrackData>
+            bool move(TrackData& td, const scalar trackTime);
+
+
+    // Friend Functions
+
+        // I-O
+
+            //- Read
+            template<class CloudType>
+            static void readFields(CloudType& c);
+
+            //- Write
+            template<class CloudType>
+            static void writeFields(const CloudType& c);
+
+
+        // Ostream operator
+
+            friend Ostream& operator<< <ParcelType>
+            (
+                Ostream&,
+                const MPPICParcel<ParcelType>&
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MPPICParcelI.H"
+#include "MPPICParcelTrackingDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MPPICParcel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
new file mode 100644
index 0000000000..ab01a6421d
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const polyMesh& owner,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
+)
+:
+    ParcelType(owner, position, cellI, tetFaceI, tetPtI),
+    UCorrect_(vector::zero)
+{}
+
+
+template<class ParcelType>
+inline Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const polyMesh& owner,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const label typeId,
+    const scalar nParticle0,
+    const scalar d0,
+    const scalar dTarget0,
+    const vector& U0,
+    const vector& UCorrect0,
+    const typename ParcelType::constantProperties& constProps
+)
+:
+    ParcelType
+    (
+        owner,
+        position,
+        cellI,
+        tetFaceI,
+        tetPtI,
+        typeId,
+        nParticle0,
+        d0,
+        dTarget0,
+        U0,
+        constProps
+    ),
+    UCorrect_(UCorrect0)
+{}
+
+
+// * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline const Foam::vector& Foam::MPPICParcel<ParcelType>::UCorrect() const
+{
+    return UCorrect_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::MPPICParcel<ParcelType>::UCorrect()
+{
+    return UCorrect_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C
new file mode 100644
index 0000000000..f0eff872b0
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "MPPICParcel.H"
+#include "IOstreams.H"
+#include "IOField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::string Foam::MPPICParcel<ParcelType>::propertyList_ =
+    Foam::MPPICParcel<ParcelType>::propertyList();
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    ParcelType(mesh, is, readFields),
+    UCorrect_(vector::zero)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is >> UCorrect_;
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&UCorrect_),
+              + sizeof(UCorrect_)
+            );
+        }
+    }
+
+    is.check
+    (
+        "MPPICParcel<ParcelType>::Collisions"
+        "(const polyMesh&, Istream&, bool)"
+    );
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::MPPICParcel<ParcelType>::readFields(CloudType& c)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    ParcelType::readFields(c);
+
+    IOField<vector> UCorrect(c.fieldIOobject("UCorrect", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, UCorrect);
+
+    label i = 0;
+
+    forAllIter(typename CloudType, c, iter)
+    {
+        MPPICParcel<ParcelType>& p = iter();
+
+        p.UCorrect_ = UCorrect[i];
+
+        i++;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::MPPICParcel<ParcelType>::writeFields(const CloudType& c)
+{
+    ParcelType::writeFields(c);
+
+    label np =  c.size();
+
+    IOField<vector>
+        UCorrect(c.fieldIOobject("UCorrect", IOobject::NO_READ), np);
+
+    label i = 0;
+
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        const MPPICParcel<ParcelType>& p = iter();
+
+        UCorrect[i] = p.UCorrect();
+
+        i++;
+    }
+
+    UCorrect.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const MPPICParcel<ParcelType>& p
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const ParcelType&>(p)
+            << token::SPACE << p.UCorrect();
+    }
+    else
+    {
+        os  << static_cast<const ParcelType&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.UCorrect_),
+          + sizeof(p.UCorrect())
+        );
+    }
+
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const MPPICParcel<ParcelType>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
new file mode 100644
index 0000000000..4f5edcc47e
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
@@ -0,0 +1,269 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "AveragingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class CloudType>
+inline Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::TrackingData
+(
+    CloudType& cloud,
+    trackPart part
+)
+:
+    ParcelType::template TrackingData<CloudType>(cloud),
+    volumeAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":volumeAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    radiusAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":radiusAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    rhoAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":rhoAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uAverage_
+    (
+        AveragingMethod<vector>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uSqrAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uSqrAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    frequencyAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":frequencyAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    massAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":massAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    part_(part)
+{}
+
+
+template<class ParcelType>
+template<class CloudType>
+inline void
+Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::updateAverages
+(
+    CloudType& cloud
+)
+{
+    // zero the sums
+    volumeAverage_() = 0;
+    radiusAverage_() = 0;
+    rhoAverage_() = 0;
+    uAverage_() = vector::zero;
+    uSqrAverage_() = 0;
+    frequencyAverage_() = 0;
+    massAverage_() = 0;
+
+    // temporary weights
+    autoPtr<AveragingMethod<scalar> > weightAveragePtr
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":weightAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    );
+    AveragingMethod<scalar>& weightAverage = weightAveragePtr();
+
+    // averaging sums
+    forAllConstIter(typename CloudType, cloud, iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh());
+
+        const scalar m = p.nParticle()*p.mass();
+
+        volumeAverage_->add(p.position(), tetIs, p.nParticle()*p.volume());
+        rhoAverage_->add(p.position(), tetIs, m*p.rho());
+        uAverage_->add(p.position(), tetIs, m*p.U());
+        massAverage_->add(p.position(), tetIs, m);
+    }
+    volumeAverage_->average();
+    massAverage_->average();
+    rhoAverage_->average(massAverage_);
+    uAverage_->average(massAverage_);
+
+    // squared velocity deviation
+    forAllConstIter(typename CloudType, cloud, iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh());
+
+        const vector u = uAverage_->interpolate(p.position(), tetIs);
+
+        uSqrAverage_->add
+        (
+            p.position(),
+            tetIs,
+            p.nParticle()*p.mass()*magSqr(p.U() - u)
+        );
+    }
+    uSqrAverage_->average(massAverage_);
+
+    // sauter mean radius
+    radiusAverage_() = volumeAverage_();
+    weightAverage = 0;
+    forAllConstIter(typename CloudType, cloud, iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh());
+
+        weightAverage.add
+        (
+            p.position(),
+            tetIs,
+            p.nParticle()*pow(p.volume(), 2.0/3.0)
+        );
+    }
+    weightAverage.average();
+    radiusAverage_->average(weightAverage);
+
+    // collision frequency
+    weightAverage = 0;
+    forAllConstIter(typename CloudType, cloud, iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+        tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), cloud.mesh());
+
+        const scalar a = volumeAverage_->interpolate(p.position(), tetIs);
+        const scalar r = radiusAverage_->interpolate(p.position(), tetIs);
+        const vector u = uAverage_->interpolate(p.position(), tetIs);
+
+        const scalar f = 0.75*a/pow3(r)*sqr(0.5*p.d() + r)*mag(p.U() - u);
+
+        frequencyAverage_->add(p.position(), tetIs, p.nParticle()*f*f);
+
+        weightAverage.add(p.position(), tetIs, p.nParticle()*f);
+    }
+    frequencyAverage_->average(weightAverage);
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+inline typename Foam::MPPICParcel<ParcelType>::template
+TrackingData<CloudType>::trackPart
+Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::part() const
+{
+    return part_;
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+inline typename Foam::MPPICParcel<ParcelType>::template
+TrackingData<CloudType>::trackPart&
+Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::part()
+{
+    return part_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/basicKinematicMPPICParcel.H b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/basicKinematicMPPICParcel.H
new file mode 100644
index 0000000000..ee9b03870e
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/basicKinematicMPPICParcel.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::basicKinematicMPPICParcel
+
+Description
+    Definition of basic kinematic MPPIC parcel
+
+SourceFiles
+    basicKinematicParcel.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicKinematicMPPICParcel_H
+#define basicKinematicMPPICParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "MPPICParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef MPPICParcel<KinematicParcel<particle> > basicKinematicMPPICParcel;
+
+    template<>
+    inline bool contiguous<basicKinematicMPPICParcel>()
+    {
+        return true;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/defineBasicKinematicMPPICParcel.C b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/defineBasicKinematicMPPICParcel.C
new file mode 100644
index 0000000000..2e0050a38f
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/defineBasicKinematicMPPICParcel.C
@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "basicKinematicMPPICParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(basicKinematicMPPICParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicKinematicMPPICParcel>, 0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.C b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.C
new file mode 100644
index 0000000000..5154df3a24
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.C
@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "basicKinematicMPPICCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// kinematic sub-models
+#include "makeParcelForces.H"
+#include "makeParcelDispersionModels.H"
+#include "makeParcelInjectionModels.H"
+#include "makeParcelPatchInteractionModels.H"
+#include "makeParcelStochasticCollisionModels.H"
+#include "makeParcelSurfaceFilmModels.H"
+
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeParcelCloudFunctionObjects(basicKinematicMPPICCloud);
+
+    // kinematic sub-models
+    makeParcelForces(basicKinematicMPPICCloud);
+    makeParcelDispersionModels(basicKinematicMPPICCloud);
+    makeParcelInjectionModels(basicKinematicMPPICCloud);
+    makeParcelPatchInteractionModels(basicKinematicMPPICCloud);
+    makeParcelStochasticCollisionModels(basicKinematicMPPICCloud);
+    makeParcelSurfaceFilmModels(basicKinematicMPPICCloud);
+
+    // MPPIC sub-models
+    makeMPPICParcelDampingModels(basicKinematicMPPICCloud);
+    makeMPPICParcelIsotropyModels(basicKinematicMPPICCloud);
+    makeMPPICParcelPackingModels(basicKinematicMPPICCloud);
+
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/include/makeMPPICParcelDampingModels.H b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelDampingModels.H
new file mode 100644
index 0000000000..9513a59cc3
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelDampingModels.H
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 makeMPPICParcelDampingModels_H
+#define makeMPPICParcelDampingModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DampingModel.H"
+#include "NoDamping.H"
+#include "Relaxation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeMPPICParcelDampingModels(CloudType)                               \
+                                                                              \
+    makeDampingModel(CloudType);                                              \
+                                                                              \
+    namespace DampingModels                                                   \
+    {                                                                         \
+        makeDampingModelType(NoDamping, CloudType);                           \
+        makeDampingModelType(Relaxation, CloudType);                          \
+    }
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/include/makeMPPICParcelIsotropyModels.H b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelIsotropyModels.H
new file mode 100644
index 0000000000..eebab863ea
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelIsotropyModels.H
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 makeMPPICParcelIsotropyModels_H
+#define makeMPPICParcelIsotropyModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoIsotropy.H"
+#include "Stochastic.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeMPPICParcelIsotropyModels(CloudType)                              \
+                                                                              \
+    makeIsotropyModel(CloudType);                                             \
+                                                                              \
+    namespace IsotropyModels                                                  \
+    {                                                                         \
+        makeIsotropyModelType(NoIsotropy, CloudType);                         \
+        makeIsotropyModelType(Stochastic, CloudType);                         \
+    }
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/include/makeMPPICParcelPackingModels.H b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelPackingModels.H
new file mode 100644
index 0000000000..2a57875b49
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/include/makeMPPICParcelPackingModels.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 makeMPPICParcelPackingModels_H
+#define makeMPPICParcelPackingModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoPacking.H"
+#include "Explicit.H"
+#include "Implicit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeMPPICParcelPackingModels(CloudType)                               \
+                                                                              \
+    makePackingModel(CloudType);                                              \
+                                                                              \
+    namespace PackingModels                                                   \
+    {                                                                         \
+        makePackingModelType(NoPacking, CloudType);                           \
+        makePackingModelType(Explicit, CloudType);                            \
+        makePackingModelType(Implicit, CloudType);                            \
+    }
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
index b961f0c4a3..5bd9be8e45 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
@@ -64,7 +64,7 @@ class PlessisMasliyahDragForce
 public:
 
     //- Runtime type information
-    TypeName("PlessisMasliyahDragForce");
+    TypeName("PlessisMasliyahDrag");
 
 
     // Constructors
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
index cd9a477193..46a6ff5e19 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
@@ -64,7 +64,7 @@ class WenYuDragForce
 public:
 
     //- Runtime type information
-    TypeName("WenYuDragForce");
+    TypeName("WenYuDrag");
 
 
     // Constructors
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
new file mode 100644
index 0000000000..21cb38d2ed
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "AveragingMethod.H"
+#include "runTimeSelectionTables.H"
+#include "pointMesh.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethod<Type>::AveragingMethod
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const labelList& size
+)
+:
+    regIOobject(io),
+    FieldField<Field, Type>(),
+    dict_(dict),
+    mesh_(mesh)
+{
+    forAll(size, i)
+    {
+        FieldField<Field, Type>::append
+        (
+            new Field<Type>(size[i], pTraits<Type>::zero)
+        );
+    }
+}
+
+
+template<class Type>
+Foam::AveragingMethod<Type>::AveragingMethod
+(
+    const AveragingMethod<Type>& am
+)
+:
+    regIOobject(am),
+    FieldField<Field, Type>(am),
+    dict_(am.dict_),
+    mesh_(am.mesh_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::autoPtr<Foam::AveragingMethod<Type> >
+Foam::AveragingMethod<Type>::New
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+{
+    word averageType(dict.lookup(typeName));
+
+    //Info<< "Selecting averaging method "
+    //    << averageType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(averageType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "Foam::AveragingMethod<Type>::New"
+            "("
+                "const dictionary&"
+                "const fvMesh&"
+            ")"
+        )   << "Unknown averaging method " << averageType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid averaging methods are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<AveragingMethod<Type> >(cstrIter()(io, dict, mesh));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethod<Type>::~AveragingMethod()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethod<Type>::updateGrad()
+{
+    // do nothing
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethod<Type>::average()
+{
+    updateGrad();
+}
+
+
+template<class Type>
+void Foam::AveragingMethod<Type>::average
+(
+    const AveragingMethod<scalar>& weight
+)
+{
+    updateGrad();
+
+    *this /= max(weight, SMALL);
+}
+
+
+template<class Type>
+bool Foam::AveragingMethod<Type>::writeData(Ostream& os) const
+{
+    return os.good();
+}
+
+
+template<class Type>
+bool Foam::AveragingMethod<Type>::write() const
+{
+    const pointMesh pointMesh_(mesh_);
+
+    // point volumes
+    Field<scalar> pointVolume(mesh_.nPoints(), 0);
+
+    // output fields
+    GeometricField<Type, fvPatchField, volMesh> cellValue
+    (
+        IOobject
+        (
+            this->name() + ":cellValue",
+            this->time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
+    );
+    GeometricField<TypeGrad, fvPatchField, volMesh> cellGrad
+    (
+        IOobject
+        (
+            this->name() + ":cellGrad",
+            this->time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensioned<TypeGrad>("zero", dimless, pTraits<TypeGrad>::zero)
+    );
+    GeometricField<Type, pointPatchField, pointMesh> pointValue
+    (
+        IOobject
+        (
+            this->name() + ":pointValue",
+            this->time().timeName(),
+            mesh_
+        ),
+        pointMesh_,
+        dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
+    );
+    GeometricField<TypeGrad, pointPatchField, pointMesh> pointGrad
+    (
+        IOobject
+        (
+            this->name() + ":pointGrad",
+            this->time().timeName(),
+            mesh_
+        ),
+        pointMesh_,
+        dimensioned<TypeGrad>("zero", dimless, pTraits<TypeGrad>::zero)
+    );
+
+    // tet-volume weighted sums
+    forAll(mesh_.C(), cellI)
+    {
+        const List<tetIndices> cellTets =
+            polyMeshTetDecomposition::cellTetIndices(mesh_, cellI);
+
+        forAll(cellTets, tetI)
+        {
+            const tetIndices& tetIs = cellTets[tetI];
+            const scalar v = tetIs.tet(mesh_).mag();
+
+            cellValue[cellI] += v*interpolate(mesh_.C()[cellI], tetIs);
+            cellGrad[cellI] += v*interpolateGrad(mesh_.C()[cellI], tetIs);
+
+            const face& f = mesh_.faces()[tetIs.face()];
+            labelList vertices(3);
+            vertices[0] = f[tetIs.faceBasePt()];
+            vertices[1] = f[tetIs.facePtA()];
+            vertices[2] = f[tetIs.facePtB()];
+
+            forAll(vertices, vertexI)
+            {
+                const label pointI = vertices[vertexI];
+
+                pointVolume[pointI] += v;
+                pointValue[pointI] +=
+                    v*interpolate(mesh_.points()[pointI], tetIs);
+                pointGrad[pointI] +=
+                    v*interpolateGrad(mesh_.points()[pointI], tetIs);
+            }
+        }
+    }
+
+    // average
+    cellValue.internalField() /= mesh_.V();
+    cellGrad.internalField() /= mesh_.V();
+    pointValue.internalField() /= pointVolume;
+    pointGrad.internalField() /= pointVolume;
+
+    // write
+    if(!cellValue.write()) return false;
+    if(!cellGrad.write()) return false;
+    if(!pointValue.write()) return false;
+    if(!pointGrad.write()) return false;
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H
new file mode 100644
index 0000000000..909459ceb4
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::AveragingMethod
+
+Description
+    Base-class for averaging lagrangian fields
+
+SourceFiles
+    AveragingMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AveragingMethod_H
+#define AveragingMethod_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class AveragingMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class AveragingMethod
+:
+    public regIOobject,
+    public FieldField<Field, Type>
+{
+protected:
+
+    //- Protected typedefs
+
+        //- Gradient type
+        typedef typename outerProduct<vector, Type>::type TypeGrad;
+
+
+    //- Protected data
+
+        //- Dictionary
+        const dictionary& dict_;
+
+        //- The mesh on which the averaging is to be done
+        const fvMesh& mesh_;
+
+
+    //- Protected member functions
+
+        //- Update the gradient calculation
+        virtual void updateGrad();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("averagingMethod");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        AveragingMethod,
+        dictionary,
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const fvMesh& mesh
+        ),
+        (io, dict, mesh)
+    );
+
+
+    //- Constructors
+
+        //- Construct from components
+        AveragingMethod
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const labelList& size
+        );
+
+        //- Construct a copy
+        AveragingMethod(const AveragingMethod<Type>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AveragingMethod<Type> > clone() const = 0;
+
+
+    //- Selector
+    static autoPtr<AveragingMethod<Type> > New
+    (
+        const IOobject& io,
+        const dictionary& dict,
+        const fvMesh& mesh
+    );
+
+
+    //- Destructor
+    virtual ~AveragingMethod();
+
+
+    //- Member Functions
+
+        //- Add point value to interpolation
+        virtual void add
+        (
+            const point position,
+            const tetIndices& tetIs,
+            const Type& value
+        ) = 0;
+
+        //- Interpolate
+        virtual Type interpolate
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const = 0;
+
+        //- Interpolate gradient
+        virtual TypeGrad interpolateGrad
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const = 0;
+
+        //- Calculate the average
+        virtual void average();
+        virtual void average(const AveragingMethod<scalar>& weight);
+
+        //- Dummy write
+        virtual bool writeData(Ostream&) const;
+
+        //- Write using setting from DB
+        virtual bool write() const;
+
+        //- Return an internal field of the average
+        virtual tmp<Field<Type> > internalField() const = 0;
+
+        //- Assign to another average
+        inline void operator=(const AveragingMethod<Type>& x);
+
+        //- Assign to value
+        inline void operator=(const Type& x);
+
+        //- Assign to tmp
+        inline void operator=(tmp<FieldField<Field, Type> > x);
+
+        //- Add-equal tmp
+        inline void operator+=(tmp<FieldField<Field, Type> > x);
+
+        //- Multiply-equal tmp
+        inline void operator*=(tmp<FieldField<Field, Type> > x);
+
+        //- Divide-equal tmp
+        inline void operator/=(tmp<FieldField<Field, scalar> > x);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "AveragingMethodI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "AveragingMethod.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethodI.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethodI.H
new file mode 100644
index 0000000000..87bc4ac904
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethodI.H
@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator=
+(
+    const AveragingMethod<Type>& x
+)
+{
+    FieldField<Field, Type>::operator=(x);
+    updateGrad();
+}
+
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator=
+(
+    const Type& x
+)
+{
+    FieldField<Field, Type>::operator=(x);
+    updateGrad();
+}
+
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator=
+(
+    tmp<FieldField<Field, Type> > x
+)
+{
+    FieldField<Field, Type>::operator=(x());
+    updateGrad();
+}
+
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator+=
+(
+    tmp<FieldField<Field, Type> > x
+)
+{
+    FieldField<Field, Type>::operator+=(x());
+    updateGrad();
+}
+
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator*=
+(
+    tmp<FieldField<Field, Type> > x
+)
+{
+    FieldField<Field, Type>::operator*=(x());
+    updateGrad();
+}
+
+
+template<class Type>
+inline void Foam::AveragingMethod<Type>::operator/=
+(
+    tmp<FieldField<Field, scalar> > x
+)
+{
+    FieldField<Field, Type>::operator/=(x());
+    updateGrad();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
new file mode 100644
index 0000000000..576153d96d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Basic.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Basic<Type>::Basic
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    AveragingMethod<Type>(io, dict, mesh, labelList(1, mesh.nCells())),
+    data_(FieldField<Field, Type>::operator[](0)),
+    dataGrad_(mesh.nCells())
+{}
+
+
+template<class Type>
+Foam::AveragingMethods::Basic<Type>::Basic
+(
+    const Basic<Type>& am
+)
+:
+    AveragingMethod<Type>(am),
+    data_(FieldField<Field, Type>::operator[](0)),
+    dataGrad_(am.dataGrad_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Basic<Type>::~Basic()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Basic<Type>::updateGrad()
+{
+    GeometricField<Type, fvPatchField, volMesh> tempData
+    (
+        IOobject
+        (
+            "BasicAverage::Data",
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh_,
+        dimensioned<Type>("zero", dimless, pTraits<Type>::zero),
+        zeroGradientFvPatchField<Type>::typeName
+    );
+    tempData.internalField() = data_;
+    tempData.correctBoundaryConditions();
+    dataGrad_ = fvc::grad(tempData)->internalField();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Basic<Type>::add
+(
+    const point position,
+    const tetIndices& tetIs,
+    const Type& value
+)
+{
+    data_[tetIs.cell()] += value/this->mesh_.V()[tetIs.cell()];
+}
+
+
+template<class Type>
+Type Foam::AveragingMethods::Basic<Type>::interpolate
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    return data_[tetIs.cell()];
+}
+
+
+template<class Type>
+typename Foam::AveragingMethods::Basic<Type>::TypeGrad
+Foam::AveragingMethods::Basic<Type>::interpolateGrad
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    return dataGrad_[tetIs.cell()];
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::AveragingMethods::Basic<Type>::internalField() const
+{
+    return tmp<Field<Type> >(data_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H
new file mode 100644
index 0000000000..d0b32e36d2
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::AveragingMethods::Basic
+
+SourceFiles
+    Basic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Basic_H
+#define Basic_H
+
+#include "AveragingMethod.H"
+#include "pointMesh.H"
+#include "tetIndices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace AveragingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class Basic Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Basic
+:
+    public AveragingMethod<Type>
+{
+public:
+
+    //- Public typedefs
+
+        //- Gradient type
+        typedef typename AveragingMethod<Type>::TypeGrad TypeGrad;
+
+
+private:
+
+    //- Private data
+
+        //- Cell average field
+        Field<Type>& data_;
+
+        //- Gradient field
+        Field<TypeGrad> dataGrad_;
+
+
+    //- Private member functions
+
+        //- Re-calculate gradient
+        virtual void updateGrad();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basic");
+
+
+    //- Constructors
+
+        //- Construct from components
+        Basic
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const fvMesh &mesh
+        );
+
+        //- Construct a copy
+        Basic(const Basic<Type>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AveragingMethod<Type> > clone() const
+        {
+            return autoPtr<AveragingMethod<Type> >
+            (
+                new Basic<Type>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Basic();
+
+
+    //- Member Functions
+
+        //- Add point value to interpolation
+        void add
+        (
+            const point position,
+            const tetIndices& tetIs,
+            const Type& value
+        );
+
+        //- Interpolate
+        Type interpolate
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Interpolate gradient
+        TypeGrad interpolateGrad
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Return an internal field of the average
+        tmp<Field<Type> > internalField() const;
+
+        //- Return an internal field of the gradient
+        tmp<Field<TypeGrad> > internalFieldGrad() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace AveragingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Basic.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
new file mode 100644
index 0000000000..51d9b5b6e5
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Dual.H"
+#include "coupledPointPatchField.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+template<class Type>
+Foam::autoPtr<labelList> Foam::AveragingMethods::Dual<Type>::size
+(
+    const fvMesh& mesh
+)
+{
+    autoPtr<labelList> s(new labelList(2));
+    s()[0] = mesh.nCells();
+    s()[1] = mesh.nPoints();
+    return s;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Dual<Type>::Dual
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    AveragingMethod<Type>(io, dict, mesh, size(mesh)),
+    volumeCell_(mesh.V()),
+    volumeDual_(mesh.nPoints(), 0.0),
+    dataCell_(FieldField<Field, Type>::operator[](0)),
+    dataDual_(FieldField<Field, Type>::operator[](1)),
+    tetVertices_(3),
+    tetCoordinates_(4)
+{
+    forAll(this->mesh_.C(), cellI)
+    {
+        List<tetIndices> cellTets =
+            polyMeshTetDecomposition::cellTetIndices(this->mesh_, cellI);
+        forAll(cellTets, tetI)
+        {
+            const tetIndices& tetIs = cellTets[tetI];
+            const face& f = this->mesh_.faces()[tetIs.face()];
+            const scalar v = tetIs.tet(this->mesh_).mag();
+            volumeDual_[f[tetIs.faceBasePt()]] += v;
+            volumeDual_[f[tetIs.facePtA()]] += v;
+            volumeDual_[f[tetIs.facePtB()]] += v;
+        }
+    }
+
+    mesh.globalData().syncPointData
+    (
+        volumeDual_,
+        plusEqOp<scalar>(),
+        mapDistribute::transform()
+    );
+}
+
+
+template<class Type>
+Foam::AveragingMethods::Dual<Type>::Dual
+(
+    const Dual<Type>& am
+)
+:
+    AveragingMethod<Type>(am),
+    volumeCell_(am.volumeCell_),
+    volumeDual_(am.volumeDual_),
+    dataCell_(FieldField<Field, Type>::operator[](0)),
+    dataDual_(FieldField<Field, Type>::operator[](1)),
+    tetVertices_(am.tetVertices_),
+    tetCoordinates_(am.tetCoordinates_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Dual<Type>::~Dual()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Dual<Type>::tetGeometry
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    const face& f = this->mesh_.faces()[tetIs.face()];
+
+    tetVertices_[0] = f[tetIs.faceBasePt()];
+    tetVertices_[1] = f[tetIs.facePtA()];
+    tetVertices_[2] = f[tetIs.facePtB()];
+
+    tetIs.tet(this->mesh_).barycentric(position, tetCoordinates_);
+
+    tetCoordinates_ = max(tetCoordinates_, 0.0);
+}
+
+
+template<class Type>
+void Foam::AveragingMethods::Dual<Type>::syncDualData()
+{
+    this->mesh_.globalData().syncPointData
+    (
+        dataDual_,
+        plusEqOp<Type>(),
+        mapDistribute::transform()
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Dual<Type>::add
+(
+    const point position,
+    const tetIndices& tetIs,
+    const Type& value
+)
+{
+    tetGeometry(position, tetIs);
+
+    dataCell_[tetIs.cell()] +=
+        tetCoordinates_[0]*value
+      / (0.25*volumeCell_[tetIs.cell()]);
+
+    for(label i = 0; i < 3; i ++)
+    {
+        dataDual_[tetVertices_[i]] +=
+            tetCoordinates_[i+1]*value
+          / (0.25*volumeDual_[tetVertices_[i]]);
+    }
+}
+
+
+template<class Type>
+Type Foam::AveragingMethods::Dual<Type>::interpolate
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    tetGeometry(position, tetIs);
+
+    return
+        tetCoordinates_[0]*dataCell_[tetIs.cell()]
+      + tetCoordinates_[1]*dataDual_[tetVertices_[0]]
+      + tetCoordinates_[2]*dataDual_[tetVertices_[1]]
+      + tetCoordinates_[3]*dataDual_[tetVertices_[2]];
+}
+
+
+template<class Type>
+typename Foam::AveragingMethods::Dual<Type>::TypeGrad
+Foam::AveragingMethods::Dual<Type>::interpolateGrad
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    tetGeometry(position, tetIs);
+
+    const label cellI(tetIs.cell());
+
+    const tensor T
+    (
+        inv
+        (
+            tensor
+            (
+                this->mesh_.points()[tetVertices_[0]] - this->mesh_.C()[cellI],
+                this->mesh_.points()[tetVertices_[1]] - this->mesh_.C()[cellI],
+                this->mesh_.points()[tetVertices_[2]] - this->mesh_.C()[cellI]
+            )
+        )
+    );
+
+    const vector t( - T.T().x() - T.T().y() - T.T().z());
+
+    const TypeGrad S
+    (
+        dataDual_[tetVertices_[0]],
+        dataDual_[tetVertices_[1]],
+        dataDual_[tetVertices_[2]]
+    );
+
+    const Type s(dataCell_[cellI]);
+
+    return (T & S) + (t*s);
+}
+
+
+template<class Type>
+void Foam::AveragingMethods::Dual<Type>::average()
+{
+    syncDualData();
+
+    AveragingMethod<Type>::average();
+}
+
+
+template<class Type>
+void Foam::AveragingMethods::Dual<Type>::average
+(
+    const AveragingMethod<scalar>& weight
+)
+{
+    syncDualData();
+
+    AveragingMethod<Type>::average(weight);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::AveragingMethods::Dual<Type>::internalField() const
+{
+    return tmp<Field<Type> >(dataCell_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H
new file mode 100644
index 0000000000..d86436fa95
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H
@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::AveragingMethods::Dual
+
+SourceFiles
+    Dual.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Dual_H
+#define Dual_H
+
+#include "AveragingMethod.H"
+#include "pointMesh.H"
+#include "tetIndices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace AveragingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class Dual Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Dual
+:
+    public AveragingMethod<Type>
+{
+public:
+
+    //- Public typedefs
+
+        //- Gradient type
+        typedef typename AveragingMethod<Type>::TypeGrad TypeGrad;
+
+
+private:
+
+    //- Private data
+
+        //- Volume of the cell-centered regions
+        const Field<scalar>& volumeCell_;
+
+        //- Volume of the point-centered regions
+        Field<scalar> volumeDual_;
+
+        //- Data on the cells
+        Field<Type>& dataCell_;
+
+        //- Data on the points
+        Field<Type>& dataDual_;
+
+        //- Tet vertex labels
+        mutable List<label> tetVertices_;
+
+        //- Tet barycentric coordinates
+        mutable List<scalar> tetCoordinates_;
+
+
+    //- Private static member functions
+    
+        //- Return the size of the FieldField parts
+        static autoPtr<labelList> size(const fvMesh &mesh);
+
+
+    //- Private member functions
+
+        //- Calculate indices and barycentric coordinates within a tetrahedron
+        void tetGeometry
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Sync point data over processor boundaries
+        void syncDualData();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dual");
+
+
+    //- Constructors
+
+        //- Construct from components
+        Dual
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const fvMesh &mesh
+        );
+
+        //- Construct a copy
+        Dual(const Dual<Type>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AveragingMethod<Type> > clone() const
+        {
+            return autoPtr<AveragingMethod<Type> >
+            (
+                new Dual<Type>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Dual();
+
+
+    //- Member Functions
+
+        //- Add point value to interpolation
+        void add
+        (
+            const point position,
+            const tetIndices& tetIs,
+            const Type& value
+        );
+
+        //- Interpolate
+        Type interpolate
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Interpolate gradient
+        TypeGrad interpolateGrad
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Calculate the average
+        void average();
+        void average(const AveragingMethod<scalar>& weight);
+
+        //- Return an internal field of the average
+        tmp<Field<Type> > internalField() const;
+
+        //- Return an internal field of the gradient
+        tmp<Field<TypeGrad> > internalFieldGrad() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace AveragingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Dual.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
new file mode 100644
index 0000000000..a71e746a1d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Moment.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Moment<Type>::Moment
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    AveragingMethod<Type>(io, dict, mesh, labelList(4, mesh.nCells())),
+    data_(FieldField<Field, Type>::operator[](0)),
+    dataX_(FieldField<Field, Type>::operator[](1)),
+    dataY_(FieldField<Field, Type>::operator[](2)),
+    dataZ_(FieldField<Field, Type>::operator[](3)),
+    transform_(mesh.nCells(), symmTensor::zero),
+    scale_(0.5*pow(mesh.V(), 1.0/3.0))
+{
+    scalar a = 1.0/24.0;
+    scalar b = 0.5854101966249685;
+    scalar c = 0.1381966011250105;
+
+    scalarField wQ(4);
+    wQ[0] = a;
+    wQ[1] = a;
+    wQ[2] = a;
+    wQ[3] = a;
+
+    vectorField xQ(4);
+    xQ[0] = vector(b, c, c);
+    xQ[1] = vector(c, b, c);
+    xQ[2] = vector(c, c, b);
+    xQ[3] = vector(c, c, c);
+
+    forAll(mesh.C(), cellI)
+    {
+        const List<tetIndices> cellTets =
+            polyMeshTetDecomposition::cellTetIndices(mesh, cellI);
+
+        symmTensor A(symmTensor::zero);
+
+        forAll(cellTets, tetI)
+        {
+            const tetIndices& tetIs = cellTets[tetI];
+            const label faceI = tetIs.face();
+            const face& f = mesh.faces()[faceI];
+
+            const tensor T
+            (
+                tensor
+                (
+                    mesh.points()[f[tetIs.faceBasePt()]] - mesh.C()[cellI],
+                    mesh.points()[f[tetIs.facePtA()]] - mesh.C()[cellI],
+                    mesh.points()[f[tetIs.facePtB()]] - mesh.C()[cellI]
+                ).T()
+            );
+
+            const vectorField d((T & xQ)/scale_[cellI]);
+
+            const scalar v(6.0*tetIs.tet(mesh).mag()/mesh.V()[cellI]);
+
+            A += v*sum(wQ*sqr(d));
+        }
+
+        transform_[cellI] = inv(A);
+    }
+}
+
+
+template<class Type>
+Foam::AveragingMethods::Moment<Type>::Moment
+(
+    const Moment<Type>& am
+)
+:
+    AveragingMethod<Type>(am),
+    data_(FieldField<Field, Type>::operator[](0)),
+    dataX_(FieldField<Field, Type>::operator[](1)),
+    dataY_(FieldField<Field, Type>::operator[](2)),
+    dataZ_(FieldField<Field, Type>::operator[](3)),
+    transform_(am.transform_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AveragingMethods::Moment<Type>::~Moment()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Moment<Type>::updateGrad()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::AveragingMethods::Moment<Type>::add
+(
+    const point position,
+    const tetIndices& tetIs,
+    const Type& value
+)
+{
+    const label cellI = tetIs.cell();
+
+    const Type v = value/this->mesh_.V()[cellI];
+    const TypeGrad dv =
+        transform_[cellI]
+      & (
+            v
+          * (position - this->mesh_.C()[cellI])
+          / scale_[cellI]
+        );
+
+    data_[cellI] += v;
+    dataX_[cellI] += v + dv.x();
+    dataY_[cellI] += v + dv.y();
+    dataZ_[cellI] += v + dv.z();
+}
+
+
+template<class Type>
+Type Foam::AveragingMethods::Moment<Type>::interpolate
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    const label cellI = tetIs.cell();
+
+    return
+        data_[cellI]
+      + (
+            TypeGrad
+            (
+                dataX_[cellI] - data_[cellI],
+                dataY_[cellI] - data_[cellI],
+                dataZ_[cellI] - data_[cellI]
+            )
+          & (position - this->mesh_.C()[cellI])
+          / scale_[cellI]
+        );
+}
+
+
+template<class Type>
+typename Foam::AveragingMethods::Moment<Type>::TypeGrad
+Foam::AveragingMethods::Moment<Type>::interpolateGrad
+(
+    const point position,
+    const tetIndices& tetIs
+) const
+{
+    const label cellI(tetIs.cell());
+
+    return
+        TypeGrad
+        (
+            dataX_[cellI] - data_[cellI],
+            dataY_[cellI] - data_[cellI],
+            dataZ_[cellI] - data_[cellI]
+        )/scale_[cellI];
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::AveragingMethods::Moment<Type>::internalField() const
+{
+    return tmp<Field<Type> >(data_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H
new file mode 100644
index 0000000000..1f8df2cf4d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::AveragingMethods::Moment
+
+SourceFiles
+    Moment.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Moment_H
+#define Moment_H
+
+#include "AveragingMethod.H"
+#include "pointMesh.H"
+#include "tetIndices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace AveragingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class Moment Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Moment
+:
+    public AveragingMethod<Type>
+{
+public:
+
+    //- Public typedefs
+
+        //- Gradient type
+        typedef typename AveragingMethod<Type>::TypeGrad TypeGrad;
+
+
+private:
+
+    //- Private data
+
+        //- Data mean
+        Field<Type>& data_;
+
+        //- X-data moment
+        Field<Type>& dataX_;
+
+        //- Y-data moment
+        Field<Type>& dataY_;
+
+        //- Z-data moment
+        Field<Type>& dataZ_;
+
+        //- Transform tensor from moment to gradient
+        Field<symmTensor> transform_;
+
+        //- Length scale for moment values
+        Field<scalar> scale_;
+
+
+    //- Private member functions
+
+        //- Re-calculate gradient
+        virtual void updateGrad();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("moment");
+
+
+    //- Constructors
+
+        //- Construct from components
+        Moment
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const fvMesh &mesh
+        );
+
+        //- Construct a copy
+        Moment(const Moment<Type>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AveragingMethod<Type> > clone() const
+        {
+            return autoPtr<AveragingMethod<Type> >
+            (
+                new Moment<Type>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Moment();
+
+
+    //- Member Functions
+
+        //- Add point value to interpolation
+        void add
+        (
+            const point position,
+            const tetIndices& tetIs,
+            const Type& value
+        );
+
+        //- Interpolate
+        Type interpolate
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Interpolate gradient
+        TypeGrad interpolateGrad
+        (
+            const point position,
+            const tetIndices& tetIs
+        ) const;
+
+        //- Return an internal field of the average
+        tmp<Field<Type> > internalField() const;
+
+        //- Return an internal field of the gradient
+        tmp<Field<TypeGrad> > internalFieldGrad() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace AveragingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Moment.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/makeAveragingMethods.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/makeAveragingMethods.C
new file mode 100644
index 0000000000..c421ab9cb8
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/makeAveragingMethods.C
@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "fvCFD.H"
+#include "polyMeshTetDecomposition.H"
+
+#include "Basic.H"
+#include "Dual.H"
+#include "Moment.H"
+
+namespace Foam
+{
+    // Scalar interpolation
+    defineNamedTemplateTypeNameAndDebug(AveragingMethod<scalar>, 0);
+    defineTemplateRunTimeSelectionTable
+    (
+        AveragingMethod<scalar>,
+        dictionary
+    );
+
+    // Vector interpolation
+    defineNamedTemplateTypeNameAndDebug(AveragingMethod<vector>, 0);
+    defineTemplateRunTimeSelectionTable
+    (
+        AveragingMethod<vector>,
+        dictionary
+    );
+
+    namespace AveragingMethods
+    {
+        // Basic interpolation
+        defineNamedTemplateTypeNameAndDebug(Basic<scalar>, 0);
+        AveragingMethod<scalar>::
+            adddictionaryConstructorToTable<Basic<scalar> >
+            addBasicscalarConstructorToTable_;
+
+        defineNamedTemplateTypeNameAndDebug(Basic<vector>, 0);
+        AveragingMethod<vector>::
+            adddictionaryConstructorToTable<Basic<vector> >
+            addBasicvectorConstructorToTable_;
+
+        // Dual interpolation
+        defineNamedTemplateTypeNameAndDebug(Dual<scalar>, 0);
+        AveragingMethod<scalar>::
+            adddictionaryConstructorToTable<Dual<scalar> >
+            addDualscalarConstructorToTable_;
+
+        defineNamedTemplateTypeNameAndDebug(Dual<vector>, 0);
+        AveragingMethod<vector>::
+            adddictionaryConstructorToTable<Dual<vector> >
+            addDualvectorConstructorToTable_;
+
+        // Moment interpolation
+        defineNamedTemplateTypeNameAndDebug(Moment<scalar>, 0);
+        AveragingMethod<scalar>::
+            adddictionaryConstructorToTable<Moment<scalar> >
+            addMomentscalarConstructorToTable_;
+
+        defineNamedTemplateTypeNameAndDebug(Moment<vector>, 0);
+        AveragingMethod<vector>::
+            adddictionaryConstructorToTable<Moment<vector> >
+            addMomentvectorConstructorToTable_;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.C
new file mode 100644
index 0000000000..84545e7640
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "CorrectionLimitingMethod.H"
+
+namespace Foam
+{
+    defineTypeNameAndDebug(CorrectionLimitingMethod, 0);
+    defineRunTimeSelectionTable(CorrectionLimitingMethod, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethod::CorrectionLimitingMethod
+(
+    const dictionary& dict
+)
+{}
+
+
+Foam::CorrectionLimitingMethod::CorrectionLimitingMethod
+(
+    const CorrectionLimitingMethod& cl
+)
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::CorrectionLimitingMethod>
+Foam::CorrectionLimitingMethod::New
+(
+    const dictionary& dict
+)
+{
+    word modelType(dict.lookup("type"));
+
+    Info<< "Selecting correction limiter " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "CorrectionLimitingMethod::New"
+            "("
+                "const word&, "
+                "const dictionary&"
+            ")"
+        )   << "Unknown correction limiter type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid correction limiter types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<CorrectionLimitingMethod>(cstrIter()(dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethod::~CorrectionLimitingMethod()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.H b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.H
new file mode 100644
index 0000000000..b5ecba77fb
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/CorrectionLimitingMethod/CorrectionLimitingMethod.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::CorrectionLimitingMethod
+
+Description
+
+SourceFiles
+    CorrectionLimitingMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CorrectionLimitingMethod_H
+#define CorrectionLimitingMethod_H
+
+#include "fvCFD.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class CorrectionLimitingMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class CorrectionLimitingMethod
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const CorrectionLimitingMethod&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("correctionLimitingMethod");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        CorrectionLimitingMethod,
+        dictionary,
+        (const dictionary& dict),
+        (dict)
+    );
+
+
+    // Constructors
+
+        //- Construct from components
+        CorrectionLimitingMethod(const dictionary& dict);
+
+        //- Construct as copy
+        CorrectionLimitingMethod(const CorrectionLimitingMethod& cl);
+
+        //- Construct and return a clone
+        virtual autoPtr<CorrectionLimitingMethod> clone() const = 0;
+
+
+    //- Selector
+    static autoPtr<CorrectionLimitingMethod> New
+    (
+        const dictionary& dict
+    );
+
+
+    //- Destructor
+    virtual ~CorrectionLimitingMethod();
+
+
+    // Member Functions
+
+        //- Return the limited velocity
+        virtual vector limitedVelocity
+        (
+            const vector uP,
+            const vector dU,
+            const vector uMean
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.C b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.C
new file mode 100644
index 0000000000..890e74f720
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "absolute.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+    defineTypeNameAndDebug(absolute, 0);
+
+    addToRunTimeSelectionTable
+    (
+        CorrectionLimitingMethod,
+        absolute,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::absolute::absolute(const dictionary& dict)
+:
+    CorrectionLimitingMethod(dict),
+    e_(readScalar(dict.lookup("e")))
+{}
+
+
+Foam::CorrectionLimitingMethods::absolute::absolute(const absolute& cl)
+:
+    CorrectionLimitingMethod(cl),
+    e_(cl.e_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::absolute::~absolute()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector Foam::CorrectionLimitingMethods::absolute::limitedVelocity
+(
+    const vector uP,
+    const vector dU,
+    const vector uMean
+) const
+{
+    const vector uRelative = uP - uMean;
+
+    return minMod
+    (
+        dU,
+      - (1.0 + this->e_)*uRelative
+       *mag(uP)/max(mag(uRelative), SMALL)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.H b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.H
new file mode 100644
index 0000000000..ef242a4c13
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/absolute/absolute.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::CorrectionLimitingMethods::absolute
+
+Description
+
+SourceFiles
+    absolute.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef absolute_H
+#define absolute_H
+
+#include "CorrectionLimitingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class absolute Declaration
+\*---------------------------------------------------------------------------*/
+
+class absolute
+:
+    public CorrectionLimitingMethod
+{
+protected:
+
+    // Protected data
+
+        //- Coefficient of restitution
+        scalar e_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("absolute");
+
+
+    // Constructors
+
+        //- Construct from components
+        absolute(const dictionary& dict);
+
+        //- Construct as copy
+        absolute(const absolute& cl);
+
+        //- Construct and return a clone
+        virtual autoPtr<CorrectionLimitingMethod> clone() const
+        {
+            return autoPtr<CorrectionLimitingMethod>
+            (
+                new absolute(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~absolute();
+
+
+    // Member Functions
+
+        //- Return the limited velocity
+        virtual vector limitedVelocity
+        (
+            const vector uP,
+            const vector dU,
+            const vector uMean
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CorrectionLimitingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.C b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.C
new file mode 100644
index 0000000000..9e05dce0c0
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.C
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "noCorrectionLimiting.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+    defineTypeNameAndDebug(noCorrectionLimiting, 0);
+
+    addToRunTimeSelectionTable
+    (
+        CorrectionLimitingMethod,
+        noCorrectionLimiting,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::noCorrectionLimiting::noCorrectionLimiting
+(
+    const dictionary& dict
+)
+:
+    CorrectionLimitingMethod(dict)
+{}
+
+
+Foam::CorrectionLimitingMethods::noCorrectionLimiting::noCorrectionLimiting
+(
+    const noCorrectionLimiting& cl
+)
+:
+    CorrectionLimitingMethod(cl)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::noCorrectionLimiting::~noCorrectionLimiting()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector
+Foam::CorrectionLimitingMethods::noCorrectionLimiting::limitedVelocity
+(
+    const vector uP,
+    const vector dU,
+    const vector uMean
+) const
+{
+    return dU;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.H b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.H
new file mode 100644
index 0000000000..b1656e2043
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/noCorrectionLimiting/noCorrectionLimiting.H
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::CorrectionLimitingMethods::noCorrectionLimiting
+
+Description
+
+SourceFiles
+    noCorrectionLimiting.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noCorrectionLimiting_H
+#define noCorrectionLimiting_H
+
+#include "CorrectionLimitingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class noCorrectionLimiting Declaration
+\*---------------------------------------------------------------------------*/
+
+class noCorrectionLimiting
+:
+    public CorrectionLimitingMethod
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from components
+        noCorrectionLimiting(const dictionary& dict);
+
+        //- Construct as copy
+        noCorrectionLimiting(const noCorrectionLimiting& cl);
+
+        //- Construct and return a clone
+        virtual autoPtr<CorrectionLimitingMethod> clone() const
+        {
+            return autoPtr<CorrectionLimitingMethod>
+            (
+                new noCorrectionLimiting(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~noCorrectionLimiting();
+
+
+    // Member Functions
+
+        //- Return the limited velocity
+        virtual vector limitedVelocity
+        (
+            const vector uP,
+            const vector dU,
+            const vector uMean
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CorrectionLimitingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.C b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.C
new file mode 100644
index 0000000000..a863ae5655
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.C
@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "relative.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+    defineTypeNameAndDebug(relative, 0);
+
+    addToRunTimeSelectionTable
+    (
+        CorrectionLimitingMethod,
+        relative,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::relative::relative(const dictionary& dict)
+:
+    CorrectionLimitingMethod(dict),
+    e_(readScalar(dict.lookup("e")))
+{}
+
+
+Foam::CorrectionLimitingMethods::relative::relative(const relative& cl)
+:
+    CorrectionLimitingMethod(cl),
+    e_(cl.e_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CorrectionLimitingMethods::relative::~relative()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector Foam::CorrectionLimitingMethods::relative::limitedVelocity
+(
+    const vector uP,
+    const vector dU,
+    const vector uMean
+) const
+{
+    const vector uRelative = uP - uMean;
+
+    return minMod
+    (
+        dU,
+      - (1.0 + this->e_)*uRelative
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.H b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.H
new file mode 100644
index 0000000000..275ab396a0
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/CorrectionLimitingMethods/relative/relative.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::CorrectionLimitingMethods::relative
+
+Description
+
+SourceFiles
+    relative.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef relative_H
+#define relative_H
+
+#include "CorrectionLimitingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace CorrectionLimitingMethods
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class relative Declaration
+\*---------------------------------------------------------------------------*/
+
+class relative
+:
+    public CorrectionLimitingMethod
+{
+protected:
+
+    // Protected data
+
+        //- Coefficient of restitution
+        scalar e_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relative");
+
+
+    // Constructors
+
+        //- Construct from components
+        relative(const dictionary& dict);
+
+        //- Construct as copy
+        relative(const relative& cl);
+
+        //- Construct and return a clone
+        virtual autoPtr<CorrectionLimitingMethod> clone() const
+        {
+            return autoPtr<CorrectionLimitingMethod>
+            (
+                new relative(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~relative();
+
+
+    // Member Functions
+
+        //- Return the limited velocity
+        virtual vector limitedVelocity
+        (
+            const vector uP,
+            const vector dU,
+            const vector uMean
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CorrectionLimitingMethods
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.C
new file mode 100644
index 0000000000..c91831966e
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.C
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "DampingModel.H"
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModel<CloudType>::DampingModel(CloudType& owner)
+:
+    CloudSubModelBase<CloudType>(owner),
+    timeScaleModel_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::DampingModel<CloudType>::DampingModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    CloudSubModelBase<CloudType>(owner, dict, typeName, type),
+    timeScaleModel_
+    (
+        TimeScaleModel::New
+        (
+            this->coeffDict().subDict(TimeScaleModel::typeName)
+        )
+    )
+{}
+
+
+template<class CloudType>
+Foam::DampingModel<CloudType>::DampingModel(const DampingModel<CloudType>& cm)
+:
+    CloudSubModelBase<CloudType>(cm),
+    timeScaleModel_(cm.timeScaleModel_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModel<CloudType>::~DampingModel()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Selector * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::DampingModel<CloudType> >
+Foam::DampingModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word modelType(dict.lookup(typeName));
+
+    Info<< "Selecting damping model " << modelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "DampingModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown damping model type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid damping model types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return
+        autoPtr<DampingModel<CloudType> >
+        (
+            cstrIter()(dict, owner)
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.H b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.H
new file mode 100644
index 0000000000..356cb72786
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.H
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::DampingModel
+
+Description
+    Templated Damping model class.
+
+SourceFiles
+    DampingModel.C
+    DampingModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DampingModel_H
+#define DampingModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "CloudSubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+class TimeScaleModel;
+
+/*---------------------------------------------------------------------------*\
+                         Class DampingModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class DampingModel
+:
+    public CloudSubModelBase<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Time scale model
+        autoPtr<TimeScaleModel> timeScaleModel_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dampingModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        DampingModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        DampingModel(CloudType& owner);
+
+        //- Construct from components
+        DampingModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        DampingModel(const DampingModel<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<DampingModel<CloudType> > clone() const = 0;
+
+
+    //- Destructor
+    virtual ~DampingModel();
+
+
+    //- Selector
+    static autoPtr<DampingModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Member Functions
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeDampingModel(CloudType)                                            \
+                                                                               \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                          \
+    defineNamedTemplateTypeNameAndDebug                                        \
+    (                                                                          \
+        DampingModel<MPPICCloudType>,                                          \
+        0                                                                      \
+    );                                                                         \
+    defineTemplateRunTimeSelectionTable                                        \
+    (                                                                          \
+        DampingModel<MPPICCloudType>,                                          \
+        dictionary                                                             \
+    );
+
+
+#define makeDampingModelType(SS, CloudType)                                    \
+                                                                               \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                          \
+    defineNamedTemplateTypeNameAndDebug(SS<MPPICCloudType>, 0);                \
+                                                                               \
+    DampingModel<MPPICCloudType>::                                             \
+        adddictionaryConstructorToTable<SS<MPPICCloudType> >                   \
+            add##SS##CloudType##MPPICCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DampingModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.C
new file mode 100644
index 0000000000..c49bf78aa5
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.C
@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "NoDamping.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModels::NoDamping<CloudType>::NoDamping
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    DampingModel<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::DampingModels::NoDamping<CloudType>::NoDamping
+(
+    const NoDamping<CloudType>& cm
+)
+:
+    DampingModel<CloudType>(cm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModels::NoDamping<CloudType>::~NoDamping()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::vector Foam::DampingModels::NoDamping<CloudType>::velocityCorrection
+(
+    typename CloudType::parcelType& p,
+    const scalar deltaT
+) const
+{
+    return vector::zero;
+}
+
+
+template<class CloudType>
+bool Foam::DampingModels::NoDamping<CloudType>::active() const
+{
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.H b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.H
new file mode 100644
index 0000000000..ee8274b783
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/NoDamping/NoDamping.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::NoDamping
+
+Description
+    Templated NoDamping model class.
+
+SourceFiles
+    NoDamping.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoDamping_H
+#define NoDamping_H
+
+#include "DampingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace DampingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class NoDamping Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoDamping
+:
+    public DampingModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+    // Constructors
+
+        //- Construct from components
+        NoDamping(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        NoDamping(const NoDamping<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<DampingModel<CloudType> > clone() const
+        {
+            return autoPtr<DampingModel<CloudType> >
+            (
+                new NoDamping<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoDamping();
+
+
+    //- Member Functions
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const;
+
+        //- Return the model 'active' status
+        virtual bool active() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace DampingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoDamping.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
new file mode 100644
index 0000000000..fadd26c67f
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Relaxation.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModels::Relaxation<CloudType>::Relaxation
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    DampingModel<CloudType>(dict, owner, typeName),
+    uAverage_(NULL),
+    oneByTimeScaleAverage_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::DampingModels::Relaxation<CloudType>::Relaxation
+(
+    const Relaxation<CloudType>& cm
+)
+:
+    DampingModel<CloudType>(cm),
+    uAverage_(NULL),
+    oneByTimeScaleAverage_(cm.oneByTimeScaleAverage_->clone())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DampingModels::Relaxation<CloudType>::
+~Relaxation()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::DampingModels::Relaxation<CloudType>::cacheFields(const bool store)
+{
+    if (store)
+    {
+        const fvMesh& mesh = this->owner().mesh();
+        const word& cloudName = this->owner().name();
+
+        const AveragingMethod<scalar>& volumeAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":volumeAverage"
+            );
+        const AveragingMethod<scalar>& radiusAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":radiusAverage"
+            );
+        const AveragingMethod<vector>& uAverage =
+            mesh.lookupObject<AveragingMethod<vector> >
+            (
+                cloudName + ":uAverage"
+            );
+        const AveragingMethod<scalar>& uSqrAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":uSqrAverage"
+            );
+        const AveragingMethod<scalar>& frequencyAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":frequencyAverage"
+            );
+
+        uAverage_ = &uAverage;
+
+        oneByTimeScaleAverage_.reset
+        (
+            AveragingMethod<scalar>::New
+            (
+                IOobject
+                (
+                    cloudName + ":oneByTimeScaleAverage",
+                    this->owner().db().time().timeName(),
+                    mesh
+                ),
+                this->owner().solution().dict(),
+                mesh
+            ).ptr()
+        );
+
+        oneByTimeScaleAverage_() =
+        (
+            this->timeScaleModel_->oneByTau
+            (
+                volumeAverage,
+                radiusAverage,
+                uSqrAverage,
+                frequencyAverage
+            )
+        )();
+    }
+    else
+    {
+        uAverage_ = NULL;
+        oneByTimeScaleAverage_.clear();
+    }
+}
+
+
+template<class CloudType>
+Foam::vector Foam::DampingModels::Relaxation<CloudType>::velocityCorrection
+(
+    typename CloudType::parcelType& p,
+    const scalar deltaT
+) const
+{
+    const tetIndices
+        tetIs(p.cell(), p.tetFace(), p.tetPt(), this->owner().mesh());
+
+    const scalar x =
+        deltaT*oneByTimeScaleAverage_->interpolate(p.position(), tetIs);
+
+    const vector u = uAverage_->interpolate(p.position(), tetIs);
+
+    return (u - p.U())*x/(x + 2.0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.H b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.H
new file mode 100644
index 0000000000..f1859cd7de
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.H
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::Relaxation
+
+Description
+    Templated Relaxation model class.
+
+SourceFiles
+    Relaxation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Relaxation_H
+#define Relaxation_H
+
+#include "DampingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace DampingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Relaxation Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class Relaxation
+:
+    public DampingModel<CloudType>
+{
+private:
+
+    // Private data
+
+        //- Velocity average
+        const AveragingMethod<vector>* uAverage_;
+
+        //- Reciprocal of the time scale average
+        autoPtr<AveragingMethod<scalar> > oneByTimeScaleAverage_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relaxation");
+
+    // Constructors
+
+        //- Construct from components
+        Relaxation(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        Relaxation(const Relaxation<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<DampingModel<CloudType> > clone() const
+        {
+            return autoPtr<DampingModel<CloudType> >
+            (
+                new Relaxation<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Relaxation();
+
+
+    //- Member Functions
+
+        //- Calculate the damping time scales
+        virtual void cacheFields(const bool store);
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace DampingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Relaxation.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.C
new file mode 100644
index 0000000000..da454553cf
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.C
@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "IsotropyModel.H"
+
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModel<CloudType>::IsotropyModel(CloudType& owner)
+:
+    CloudSubModelBase<CloudType>(owner),
+    timeScaleModel_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::IsotropyModel<CloudType>::IsotropyModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    CloudSubModelBase<CloudType>(owner, dict, typeName, type),
+    timeScaleModel_
+    (
+        TimeScaleModel::New
+        (
+            this->coeffDict().subDict(TimeScaleModel::typeName)
+        )
+    )
+{}
+
+
+template<class CloudType>
+Foam::IsotropyModel<CloudType>::IsotropyModel
+(
+    const IsotropyModel<CloudType>& cm
+)
+:
+    CloudSubModelBase<CloudType>(cm),
+    timeScaleModel_(cm.timeScaleModel_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModel<CloudType>::~IsotropyModel()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Selector * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::IsotropyModel<CloudType> >
+Foam::IsotropyModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word modelType(dict.lookup(typeName));
+
+    Info<< "Selecting isotropy model " << modelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "IsotropyModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown isotropy model type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid isotropy model types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return
+        autoPtr<IsotropyModel<CloudType> >
+        (
+            cstrIter()(dict, owner)
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.H b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.H
new file mode 100644
index 0000000000..ad3b500363
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.H
@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::IsotropyModel
+
+Description
+    Templated Isotropy model class.
+
+SourceFiles
+    IsotropyModel.C
+    IsotropyModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IsotropyModel_H
+#define IsotropyModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "CloudSubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+class TimeScaleModel;
+
+/*---------------------------------------------------------------------------*\
+                         Class IsotropyModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class IsotropyModel
+:
+    public CloudSubModelBase<CloudType>
+{
+protected:
+
+        //- Time scale model
+        autoPtr<TimeScaleModel> timeScaleModel_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("isotropyModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        IsotropyModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        IsotropyModel(CloudType& owner);
+
+        //- Construct from components
+        IsotropyModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct a copy
+        IsotropyModel(const IsotropyModel<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<IsotropyModel<CloudType> > clone() const = 0;
+
+
+    //- Destructor
+    virtual ~IsotropyModel();
+
+
+    //- Selector
+    static autoPtr<IsotropyModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    //- Member Functions
+
+        //- Calculate velocities
+        virtual void calculate() = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeIsotropyModel(CloudType)                                \
+                                                                              \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        IsotropyModel<MPPICCloudType>,                             \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        IsotropyModel<MPPICCloudType>,                             \
+        dictionary                                                            \
+    );
+
+
+#define makeIsotropyModelType(SS, CloudType)                       \
+                                                                              \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<MPPICCloudType>, 0);               \
+                                                                              \
+    IsotropyModel<MPPICCloudType>::                                \
+        adddictionaryConstructorToTable<SS<MPPICCloudType> >                  \
+            add##SS##CloudType##MPPICCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "IsotropyModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.C
new file mode 100644
index 0000000000..d8b0fb075d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "NoIsotropy.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModels::NoIsotropy<CloudType>::NoIsotropy
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    IsotropyModel<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::IsotropyModels::NoIsotropy<CloudType>::NoIsotropy
+(
+    const NoIsotropy<CloudType>& cm
+)
+:
+    IsotropyModel<CloudType>(cm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModels::NoIsotropy<CloudType>::~NoIsotropy()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::IsotropyModels::NoIsotropy<CloudType>::calculate()
+{
+    // do nothing
+}
+
+
+template<class CloudType>
+bool Foam::IsotropyModels::NoIsotropy<CloudType>::active() const
+{
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.H b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.H
new file mode 100644
index 0000000000..071b01d53c
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/NoIsotropy/NoIsotropy.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::NoIsotropy
+
+Description
+    Templated NoIsotropy model class.
+
+SourceFiles
+    NoIsotropy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoIsotropy_H
+#define NoIsotropy_H
+
+#include "IsotropyModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace IsotropyModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class NoIsotropy Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoIsotropy
+:
+    public IsotropyModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+    // Constructors
+
+        //- Construct from components
+        NoIsotropy(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        NoIsotropy(const NoIsotropy<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<IsotropyModel<CloudType> > clone() const
+        {
+            return autoPtr<IsotropyModel<CloudType> >
+            (
+                new NoIsotropy<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoIsotropy();
+
+
+    //- Member Functions
+
+        //- calculate velocities
+        virtual void calculate();
+
+        //- Return the model 'active' status
+        virtual bool active() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IsotropyModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoIsotropy.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
new file mode 100644
index 0000000000..fcf022e5d2
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
@@ -0,0 +1,257 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Stochastic.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModels::Stochastic<CloudType>::Stochastic
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    IsotropyModel<CloudType>(dict, owner, typeName)
+{}
+
+
+template<class CloudType>
+Foam::IsotropyModels::Stochastic<CloudType>::Stochastic
+(
+    const Stochastic<CloudType>& cm
+)
+:
+    IsotropyModel<CloudType>(cm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::IsotropyModels::Stochastic<CloudType>::~Stochastic()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::IsotropyModels::Stochastic<CloudType>::sampleGauss()
+{
+    static bool isCached = true;
+    static scalar xCached;
+
+    if (isCached)
+    {
+        isCached = false;
+
+        return xCached;
+    }
+    else
+    {
+        cachedRandom& rndGen = this->owner().rndGen();
+
+        scalar f, m, x, y;
+
+        do
+        {
+            x = 2.0*rndGen.sample01<scalar>() - 1.0;
+            y = 2.0*rndGen.sample01<scalar>() - 1.0;
+            m = x*x + y*y;
+        } while (m >= 1.0 || m == 0.0);
+
+        f = sqrt(-2.0*log(m)/m);
+        xCached = x*f;
+        isCached = true;
+
+        return y*f;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
+{
+    const fvMesh& mesh = this->owner().mesh();
+    const scalar deltaT(this->owner().db().time().deltaTValue());
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    const scalar oneBySqrtThree = sqrt(1.0/3.0);
+
+    const AveragingMethod<scalar>& volumeAverage =
+        mesh.lookupObject<AveragingMethod<scalar> >
+        (
+            this->owner().name() + ":volumeAverage"
+        );
+    const AveragingMethod<scalar>& radiusAverage =
+        mesh.lookupObject<AveragingMethod<scalar> >
+        (
+            this->owner().name() + ":radiusAverage"
+        );
+    const AveragingMethod<vector>& uAverage =
+        mesh.lookupObject<AveragingMethod<vector> >
+        (
+            this->owner().name() + ":uAverage"
+        );
+    const AveragingMethod<scalar>& uSqrAverage =
+        mesh.lookupObject<AveragingMethod<scalar> >
+        (
+            this->owner().name() + ":uSqrAverage"
+        );
+    const AveragingMethod<scalar>& frequencyAverage =
+        mesh.lookupObject<AveragingMethod<scalar> >
+        (
+            this->owner().name() + ":frequencyAverage"
+        );
+    const AveragingMethod<scalar>& massAverage =
+        mesh.lookupObject<AveragingMethod<scalar> >
+        (
+            this->owner().name() + ":massAverage"
+        );
+
+    // calculate time scales and pdf exponent
+    autoPtr<AveragingMethod<scalar> > exponentAveragePtr
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                this->owner().name() + ":exponentAverage",
+                this->owner().db().time().timeName(),
+                mesh
+            ),
+            this->owner().solution().dict(),
+            mesh
+        )
+    );
+    AveragingMethod<scalar>& exponentAverage = exponentAveragePtr();
+    exponentAverage =
+        exp
+        (
+          - deltaT
+           *this->timeScaleModel_->oneByTau
+            (
+                volumeAverage,
+                radiusAverage,
+                uSqrAverage,
+                frequencyAverage
+            )
+        )();
+
+    // random sampling
+    forAllIter(typename CloudType, this->owner(), iter)
+    {
+        typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
+
+        const scalar x = exponentAverage.interpolate(p.position(), tetIs);
+
+        if(x < rndGen.sample01<scalar>())
+        {
+            const vector r(sampleGauss(), sampleGauss(), sampleGauss());
+
+            const vector u = uAverage.interpolate(p.position(), tetIs);
+            const scalar uRms =
+                sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
+
+            p.U() = u + r*uRms*oneBySqrtThree;
+        }
+    }
+
+    // correction velocity averages
+    autoPtr<AveragingMethod<vector> > uTildeAveragePtr
+    (
+        AveragingMethod<vector>::New
+        (
+            IOobject
+            (
+                this->owner().name() + ":uTildeAverage",
+                this->owner().db().time().timeName(),
+                mesh
+            ),
+            this->owner().solution().dict(),
+            mesh
+        )
+    );
+    AveragingMethod<vector>& uTildeAverage = uTildeAveragePtr();
+    forAllIter(typename CloudType, this->owner(), iter)
+    {
+        typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
+        uTildeAverage.add(p.position(), tetIs, p.nParticle()*p.mass()*p.U());
+    }
+    uTildeAverage.average(massAverage);
+
+    autoPtr<AveragingMethod<scalar> > uTildeSqrAveragePtr
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                this->owner().name() + ":uTildeSqrAverage",
+                this->owner().db().time().timeName(),
+                mesh
+            ),
+            this->owner().solution().dict(),
+            mesh
+        )
+    );
+    AveragingMethod<scalar>& uTildeSqrAverage = uTildeSqrAveragePtr();
+    forAllIter(typename CloudType, this->owner(), iter)
+    {
+        typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
+        const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs);
+        uTildeSqrAverage.add
+        (
+            p.position(),
+            tetIs,
+            p.nParticle()*p.mass()*magSqr(p.U() - uTilde)
+        );
+    }
+    uTildeSqrAverage.average(massAverage);
+
+    // conservation correction
+    forAllIter(typename CloudType, this->owner(), iter)
+    {
+        typename CloudType::parcelType& p = iter();
+        const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
+
+        const vector u = uAverage.interpolate(p.position(), tetIs);
+        const scalar uRms =
+            sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
+
+        const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs);
+        const scalar uTildeRms =
+            sqrt(max(uTildeSqrAverage.interpolate(p.position(), tetIs), 0.0));
+
+        p.U() = u + (p.U() - uTilde)*uRms/max(uTildeRms, SMALL);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.H b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.H
new file mode 100644
index 0000000000..fb3893fa20
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.H
@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::Stochastic
+
+Description
+    Templated Stochastic model class.
+
+SourceFiles
+    Stochastic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Stochastic_H
+#define Stochastic_H
+
+#include "IsotropyModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace IsotropyModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Stochastic Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class Stochastic
+:
+    public IsotropyModel<CloudType>
+{
+private:
+
+    // Private Member Functions
+
+        //- Sample a gaussian distribution using the Box-Muller method
+        scalar sampleGauss();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("stochastic");
+
+    // Constructors
+
+        //- Construct from components
+        Stochastic(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        Stochastic(const Stochastic<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<IsotropyModel<CloudType> > clone() const
+        {
+            return autoPtr<IsotropyModel<CloudType> >
+            (
+                new Stochastic<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Stochastic();
+
+
+    //- Member Functions
+
+        //- Calculate velocities
+        virtual void calculate();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IsotropyModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Stochastic.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
new file mode 100644
index 0000000000..f45df3034d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Explicit.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::Explicit<CloudType>::Explicit
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PackingModel<CloudType>(dict, owner, typeName),
+    stressAverage_(NULL),
+    correctionLimiting_
+    (
+        CorrectionLimitingMethod::New
+        (
+            this->coeffDict().subDict(CorrectionLimitingMethod::typeName)
+        )
+    )
+{}
+
+
+template<class CloudType>
+Foam::PackingModels::Explicit<CloudType>::Explicit
+(
+    const Explicit<CloudType>& cm
+)
+:
+    PackingModel<CloudType>(cm),
+    stressAverage_(cm.stressAverage_->clone()),
+    correctionLimiting_
+    (
+        cm.correctionLimiting_->clone()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::Explicit<CloudType>::~Explicit()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PackingModels::Explicit<CloudType>::cacheFields(const bool store)
+{
+    PackingModel<CloudType>::cacheFields(store);
+
+    if (store)
+    {
+        const fvMesh& mesh = this->owner().mesh();
+        const word& cloudName = this->owner().name();
+
+        const AveragingMethod<scalar>& volumeAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":volumeAverage"
+            );
+        const AveragingMethod<scalar>& rhoAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":rhoAverage"
+            );
+        const AveragingMethod<vector>& uAverage =
+            mesh.lookupObject<AveragingMethod<vector> >
+            (
+                cloudName + ":uAverage"
+            );
+        const AveragingMethod<scalar>& uSqrAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":uSqrAverage"
+            );
+
+        volumeAverage_ = &volumeAverage;
+        uAverage_ = &uAverage;
+
+        stressAverage_.reset
+        (
+            AveragingMethod<scalar>::New
+            (
+                IOobject
+                (
+                    cloudName + ":stressAverage",
+                    this->owner().db().time().timeName(),
+                    mesh
+                ),
+                this->owner().solution().dict(),
+                mesh
+            ).ptr()
+        );
+
+        stressAverage_() =
+            this->particleStressModel_->tau
+            (
+                *volumeAverage_,
+                rhoAverage,
+                uSqrAverage
+            )();
+    }
+    else
+    {
+        volumeAverage_ = NULL;
+        uAverage_ = NULL;
+        stressAverage_.clear();
+    }
+}
+
+
+template<class CloudType>
+Foam::vector Foam::PackingModels::Explicit<CloudType>::velocityCorrection
+(
+    typename CloudType::parcelType& p,
+    const scalar deltaT
+) const
+{
+    const fvMesh& mesh = this->owner().mesh();
+    const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
+
+    // interpolated quantities
+    const scalar alpha =
+        this->volumeAverage_->interpolate(p.position(), tetIs);
+    const vector alphaGrad =
+        this->volumeAverage_->interpolateGrad(p.position(), tetIs);
+    const vector uMean =
+        this->uAverage_->interpolate(p.position(), tetIs);
+
+    // stress gradient
+    const vector tauGrad =
+        stressAverage_->interpolateGrad(p.position(), tetIs);
+
+    // parcel relative velocity
+    const vector uRelative = p.U() - uMean;
+
+    // correction velocity
+    vector dU = vector::zero;
+
+    //// existing forces
+    //const scalar Re = p.Re(p.U(), p.d(), p.rhoc(), p.muc());
+    //const typename CloudType::forceType& forces = this->owner().forces();
+    //const forceSuSp F =
+    //    forces.calcCoupled(p, deltaT, p.mass(), Re, p.muc())
+    //  + forces.calcNonCoupled(p, deltaT, p.mass(), Re, p.muc());
+
+    // correction velocity
+    if((uRelative & alphaGrad) > 0)
+    {
+        dU = - deltaT*tauGrad/(p.rho()*alpha/* + deltaT*F.Sp()*/);
+    }
+
+    // apply the velocity limiters
+    return
+        correctionLimiting_->limitedVelocity
+        (
+            p.U(),
+            dU,
+            uMean
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.H b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.H
new file mode 100644
index 0000000000..b3dd1a5571
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::Explicit
+
+SourceFiles
+    Explicit.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Explicit_H
+#define Explicit_H
+
+#include "PackingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PackingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Explicit Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class Explicit
+:
+    public PackingModel<CloudType>
+{
+private:
+
+    //- Private data
+
+        //- Volume fraction average
+        const AveragingMethod<scalar>* volumeAverage_;
+
+        //- Velocity average
+        const AveragingMethod<vector>* uAverage_;
+
+        //- Stress average field
+        autoPtr<AveragingMethod<scalar> > stressAverage_;
+
+        //- Correction limiter
+        autoPtr<CorrectionLimitingMethod> correctionLimiting_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("explicit");
+
+    // Constructors
+
+        //- Construct from components
+        Explicit(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        Explicit(const Explicit<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PackingModel<CloudType> > clone() const
+        {
+            return autoPtr<PackingModel<CloudType> >
+            (
+                new Explicit<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Explicit();
+
+
+    // Member Functions
+
+        //- Calculate the inter particles stresses
+        virtual void cacheFields(const bool store);
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PackingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Explicit.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
new file mode 100644
index 0000000000..d536dc047f
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -0,0 +1,303 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Implicit.H"
+#include "fixedValueFvsPatchField.H"
+#include "volPointInterpolation.H"                                              
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::Implicit<CloudType>::Implicit
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PackingModel<CloudType>(dict, owner, typeName),
+    alpha_
+    (
+        this->owner().name() + ":alpha",
+        this->owner().theta()
+    ),
+    phiCorrect_(NULL),
+    uCorrect_(NULL),
+    applyGravity_(this->coeffDict().lookup("applyGravity")),
+    alphaMin_(readScalar(this->coeffDict().lookup("alphaMin"))),
+    rhoMin_(readScalar(this->coeffDict().lookup("rhoMin")))
+{
+    alpha_.oldTime();
+}
+
+
+template<class CloudType>
+Foam::PackingModels::Implicit<CloudType>::Implicit
+(
+    const Implicit<CloudType>& cm
+)
+:
+    PackingModel<CloudType>(cm),
+    alpha_(cm.alpha_),
+    phiCorrect_(cm.phiCorrect_()),
+    uCorrect_(cm.uCorrect_()),
+    applyGravity_(cm.applyGravity_),
+    alphaMin_(cm.alphaMin_),
+    rhoMin_(cm.rhoMin_)
+{
+    alpha_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::Implicit<CloudType>::~Implicit()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
+{
+    PackingModel<CloudType>::cacheFields(store);
+
+    if (store)
+    { 
+        const fvMesh& mesh = this->owner().mesh();
+        const dimensionedScalar deltaT = this->owner().db().time().deltaT();
+        const word& cloudName = this->owner().name();
+
+        const dimensionedVector& g = this->owner().g();
+        const volScalarField& rhoc = this->owner().rho();
+
+        const AveragingMethod<scalar>& rhoAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":rhoAverage"
+            );
+        const AveragingMethod<scalar>& uSqrAverage =
+            mesh.lookupObject<AveragingMethod<scalar> >
+            (
+                cloudName + ":uSqrAverage"
+            );
+
+
+        // Property fields
+        // ~~~~~~~~~~~~~~~
+
+        // volume fraction field
+        alpha_ = max(this->owner().theta(), alphaMin_);
+        alpha_.correctBoundaryConditions();
+
+        // average density
+        volScalarField rho
+        (
+            IOobject
+            (
+                cloudName + ":rho",
+                this->owner().db().time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimDensity, 0),
+            zeroGradientFvPatchField<scalar>::typeName
+        );
+        rho.internalField() = max(rhoAverage.internalField(), rhoMin_);
+        rho.correctBoundaryConditions();
+
+        //Info << "       x: " << mesh.C().internalField().component(2) << endl;
+        //Info << "   alpha: " << alpha_.internalField() << endl;
+        //Info << "alphaOld: " << alpha_.oldTime().internalField() << endl;
+        //Info << "     rho: " << rho.internalField() << endl;
+        //Info << endl;
+
+        
+        // Stress field
+        // ~~~~~~~~~~~~
+
+        // stress derivative wrt volume fraction
+        volScalarField tauPrime
+        (
+            IOobject
+            (
+                cloudName + ":tauPrime",
+                this->owner().db().time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimPressure, 0),
+            zeroGradientFvPatchField<scalar>::typeName
+        );
+
+        tauPrime.internalField() =
+            this->particleStressModel_->dTaudTheta
+            (
+                alpha_.internalField(),
+                rho.internalField(),
+                uSqrAverage.internalField()
+            )();
+
+        tauPrime.correctBoundaryConditions();
+
+
+        // Implicit solution for the volume fraction
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        surfaceScalarField
+            tauPrimeByRhoAf
+            (
+                "tauPrimeByRhoAf",
+                fvc::interpolate(deltaT*tauPrime/rho)
+            );
+
+        fvScalarMatrix alphaEqn
+        (
+            fvm::ddt(alpha_)
+          - fvc::ddt(alpha_)
+          - fvm::laplacian(tauPrimeByRhoAf, alpha_)
+        );
+
+        if (applyGravity_)
+        {
+            surfaceScalarField
+                phiGByA
+                (
+                    "phiGByA",
+                    deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
+                );
+
+            alphaEqn += fvm::div(phiGByA, alpha_);
+        }
+
+        alphaEqn.solve();
+
+        //// updated stress
+        //tauPrime.internalField() =
+        //    this->particleStressModel_->tauPrime
+        //    (
+        //        alpha_.internalField(),
+        //        rho.internalField(),
+        //        uSqrAverage.internalField()
+        //    )();
+        //tauPrime.correctBoundaryConditions();
+
+
+        // Generate correction fields
+        // ~~~~~~~~~~~~~~~~~
+
+        // correction volumetric flux
+        phiCorrect_ = tmp<surfaceScalarField>
+        (
+            new surfaceScalarField
+            (
+                cloudName + ":phiCorrect",
+                alphaEqn.flux()/fvc::interpolate(alpha_)
+            )
+        );
+
+        // correction velocity
+        uCorrect_ = tmp<volVectorField>
+        (
+            new volVectorField
+            (
+                cloudName + ":uCorrect",
+                fvc::reconstruct(phiCorrect_())
+            //  - deltaT*fvc::grad(tauPrime)/(rho*alpha_)
+            //  + (applyGravity_ ? deltaT*g*(1.0 - rhoc/rho) : 0.0)
+            )
+
+        );
+        uCorrect_->correctBoundaryConditions();
+
+        //Info << endl;
+        //Info << "     alpha: " << alpha_.internalField() << endl;
+        //Info << "phiCorrect: " << phiCorrect_->internalField() << endl;
+        //Info << "  uCorrect: " << uCorrect_->internalField() << endl;
+        //Info << endl;
+    }
+    else
+    {
+        alpha_.oldTime();
+        phiCorrect_.clear();
+        uCorrect_.clear();
+    }
+}
+
+
+template<class CloudType>
+Foam::vector Foam::PackingModels::Implicit<CloudType>::velocityCorrection
+(
+    typename CloudType::parcelType& p,
+    const scalar deltaT
+) const
+{
+    const fvMesh& mesh = this->owner().mesh();
+
+    // containing tetrahedron and parcel coordinates within
+    const label cellI = p.cell();
+    const label faceI = p.tetFace();
+    const tetIndices tetIs(cellI, faceI, p.tetPt(), mesh);
+    List<scalar> tetCoordinates(4);
+    tetIs.tet(mesh).barycentric(p.position(), tetCoordinates);
+
+    // cell velocity
+    const vector U = uCorrect_()[cellI];
+
+    // face geometry
+    vector nHat = mesh.faces()[faceI].normal(mesh.points());
+    const scalar nMag = mag(nHat);
+    nHat /= nMag;
+
+    // get face flux
+    scalar phi;
+    const label patchI = mesh.boundaryMesh().whichPatch(faceI);
+    if (patchI == -1)
+    {
+        phi = phiCorrect_()[faceI];
+    }
+    else
+    {
+        phi =
+            phiCorrect_().boundaryField()[patchI]
+            [
+                mesh.boundaryMesh()[patchI].whichFace(faceI)
+            ];
+    }
+
+    // interpolant equal to 1 at the cell centre and 0 at the face
+    const scalar t = tetCoordinates[0];
+
+    // the normal component of the velocity correction is interpolated linearly
+    // the tangential component is equal to that at the cell centre
+    return U + (1.0 - t)*nHat*(phi/nMag - (U & nHat));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H
new file mode 100644
index 0000000000..095fed016b
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::Implicit
+
+SourceFiles
+    Implicit.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Implicit_H
+#define Implicit_H
+
+#include "PackingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PackingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Implicit Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class Implicit
+:
+    public PackingModel<CloudType>
+{
+private:
+
+    //- Private data
+            
+        //- Volume fraction field
+        volScalarField alpha_;
+
+        //- Correction flux
+        tmp<surfaceScalarField> phiCorrect_;
+
+        //- Correction cell-centred velocity
+        tmp<volVectorField> uCorrect_;
+
+        //- Flag to indicate whether gravity is applied
+        Switch applyGravity_;
+
+        //- Minimum stable volume fraction
+        scalar alphaMin_;
+
+        //- Minimum stable density
+        scalar rhoMin_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("implicit");
+
+    // Constructors
+
+        //- Construct from components
+        Implicit(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        Implicit(const Implicit<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PackingModel<CloudType> > clone() const
+        {
+            return autoPtr<PackingModel<CloudType> >
+            (
+                new Implicit<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Implicit();
+
+
+    // Member Functions
+
+        //- Calculate the inter particles stresses
+        virtual void cacheFields(const bool store);
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PackingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Implicit.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.C
new file mode 100644
index 0000000000..92adbba777
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.C
@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "NoPacking.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::NoPacking<CloudType>::NoPacking
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PackingModel<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::PackingModels::NoPacking<CloudType>::NoPacking
+(
+    const NoPacking<CloudType>& cm
+)
+:
+    PackingModel<CloudType>(cm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModels::NoPacking<CloudType>::~NoPacking()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::vector Foam::PackingModels::NoPacking<CloudType>::velocityCorrection
+(
+    typename CloudType::parcelType& p,
+    const scalar deltaT
+) const
+{
+    return vector::zero;
+}
+
+
+template<class CloudType>
+bool Foam::PackingModels::NoPacking<CloudType>::active() const
+{
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.H b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.H
new file mode 100644
index 0000000000..a7aa5ce84e
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/NoPacking/NoPacking.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::NoPacking
+
+Description
+    Templated NoPacking model class.
+
+SourceFiles
+    NoPacking.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoPacking_H
+#define NoPacking_H
+
+#include "PackingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PackingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class NoPacking Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoPacking
+:
+    public PackingModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+    // Constructors
+
+        //- Construct from components
+        NoPacking(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        NoPacking(const NoPacking<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PackingModel<CloudType> > clone() const
+        {
+            return autoPtr<PackingModel<CloudType> >
+            (
+                new NoPacking<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoPacking();
+
+
+    // Member Functions
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const;
+
+        //- Return the model 'active' status
+        virtual bool active() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PackingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoPacking.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.C
new file mode 100644
index 0000000000..fd29b22769
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.C
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "PackingModel.H"
+#include "AveragingMethod.H"
+#include "ParticleStressModel.H"
+#include "CorrectionLimitingMethod.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModel<CloudType>::PackingModel(CloudType& owner)
+:
+    CloudSubModelBase<CloudType>(owner),
+    particleStressModel_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::PackingModel<CloudType>::PackingModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    CloudSubModelBase<CloudType>(owner, dict, typeName, type),
+    particleStressModel_
+    (
+        ParticleStressModel::New
+        (
+            this->coeffDict().subDict(ParticleStressModel::typeName)
+        )
+    )
+{}
+
+
+template<class CloudType>
+Foam::PackingModel<CloudType>::PackingModel(const PackingModel<CloudType>& cm)
+:
+    CloudSubModelBase<CloudType>(cm),
+    particleStressModel_(cm.particleStressModel_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PackingModel<CloudType>::~PackingModel()
+{}
+
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::PackingModel<CloudType> >
+Foam::PackingModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word modelType(dict.lookup(typeName));
+
+    Info<< "Selecting packing model " << modelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "PackingModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown packing model type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid packing model types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<PackingModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.H b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.H
new file mode 100644
index 0000000000..d29af9c449
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.H
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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::PackingModel
+
+Description
+    Templated Packing model class.
+
+SourceFiles
+    PackingModel.C
+    PackingModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PackingModel_H
+#define PackingModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "CloudSubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+class ParticleStressModel;
+
+class CorrectionLimitingMethod;
+
+template <class Type>
+class AveragingMethod;
+
+/*---------------------------------------------------------------------------*\
+                         Class PackingModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class PackingModel
+:
+    public CloudSubModelBase<CloudType>
+{
+protected:
+
+    //- Protected data
+
+        //- Particle stress model
+        autoPtr<ParticleStressModel> particleStressModel_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("packingModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        PackingModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        PackingModel(CloudType& owner);
+
+        //- Construct from components
+        PackingModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        PackingModel(const PackingModel<CloudType>& cm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PackingModel<CloudType> > clone() const = 0;
+
+
+    //- Destructor
+    virtual ~PackingModel();
+
+
+    //- Selector
+    static autoPtr<PackingModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Member Functions
+
+        //- Calculate the velocity correction
+        virtual vector velocityCorrection
+        (
+            typename CloudType::parcelType& p,
+            const scalar deltaT
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makePackingModel(CloudType)                                           \
+                                                                              \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        PackingModel<MPPICCloudType>,                                         \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        PackingModel<MPPICCloudType>,                                         \
+        dictionary                                                            \
+    );
+
+
+#define makePackingModelType(SS, CloudType)                                   \
+                                                                              \
+    typedef CloudType::MPPICCloudType MPPICCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<MPPICCloudType>, 0);               \
+                                                                              \
+    PackingModel<MPPICCloudType>::                                            \
+        adddictionaryConstructorToTable<SS<MPPICCloudType> >                  \
+            add##SS##CloudType##MPPICCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PackingModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.C b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.C
new file mode 100644
index 0000000000..ec7715d230
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.C
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "HarrisCrighton.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+    defineTypeNameAndDebug(HarrisCrighton, 0);
+
+    addToRunTimeSelectionTable
+    (
+        ParticleStressModel,
+        HarrisCrighton,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::HarrisCrighton::HarrisCrighton
+(
+    const dictionary& dict
+)
+:
+    ParticleStressModel(dict),
+    pSolid_(readScalar(dict.lookup("pSolid"))),
+    beta_(readScalar(dict.lookup("beta"))),
+    eps_(readScalar(dict.lookup("eps")))
+{}
+
+
+Foam::ParticleStressModels::HarrisCrighton::HarrisCrighton
+(
+    const HarrisCrighton& hc
+)
+:
+    ParticleStressModel(hc),
+    pSolid_(hc.pSolid_),
+    beta_(hc.beta_),
+    eps_(hc.eps_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::HarrisCrighton::~HarrisCrighton()
+{}
+
+
+// * * * * * * * * * * * * * Privare Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::HarrisCrighton::denominator
+(
+    const Field<scalar>& alpha
+) const
+{
+    return
+        max
+        (
+            alphaPacked_ - alpha,
+            max(eps_*(1.0 - alpha), SMALL)
+        );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::HarrisCrighton::tau
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    return
+    (
+        pSolid_
+      * pow(alpha, beta_)
+      / denominator(alpha)
+    );
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::HarrisCrighton::dTaudTheta
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    const Field<scalar> d(denominator(alpha));
+
+    return
+    (
+        pSolid_
+      * pow(alpha, beta_)
+      / d
+      * (beta_/alpha + 1.0/d)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.H b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.H
new file mode 100644
index 0000000000..e9b6f0210c
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/HarrisCrighton/HarrisCrighton.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 HarrisCrighton_H
+#define HarrisCrighton_H
+
+#include "ParticleStressModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class HarrisCrighton Declaration
+\*---------------------------------------------------------------------------*/
+
+class HarrisCrighton
+:
+    public ParticleStressModel
+{
+    // Private data
+
+        //- Solid pressure coefficient
+        scalar pSolid_;
+
+        //- Exponent of the volume fraction
+        scalar beta_;
+
+        //- Smallest allowable difference from the packed volume fraction
+        scalar eps_;
+
+
+    // Private member functions
+
+        //- Return the limited denominator of the radial distribution function
+        tmp<Field<scalar> > denominator(const Field<scalar>& alpha) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("HarrisCrighton");
+
+
+    //- Constructors
+
+        //- Construct from components
+        HarrisCrighton(const dictionary& dict);
+
+        //- Construct copy
+        HarrisCrighton(const HarrisCrighton& hc);
+
+        //- Clone
+        virtual autoPtr<ParticleStressModel> clone() const
+        {
+            return autoPtr<ParticleStressModel>
+            (
+                new HarrisCrighton(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~HarrisCrighton();
+
+
+    //- Member Functions
+
+        //- Collision stress
+        tmp<Field<scalar> > tau
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+
+        //- Collision stress derivaive w.r.t. the volume fraction
+        tmp<Field<scalar> > dTaudTheta
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace ParticleStressModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.C b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.C
new file mode 100644
index 0000000000..9962149fbd
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "Lun.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+    defineTypeNameAndDebug(Lun, 0);
+
+    addToRunTimeSelectionTable
+    (
+        ParticleStressModel,
+        Lun,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::Lun::Lun
+(
+    const dictionary& dict
+)
+:
+    ParticleStressModel(dict),
+    e_(readScalar(dict.lookup("e"))),
+    eps_(readScalar(dict.lookup("eps")))
+{}
+
+
+Foam::ParticleStressModels::Lun::Lun
+(
+    const Lun& ln
+)
+:
+    ParticleStressModel(ln),
+    e_(ln.e_),
+    eps_(ln.eps_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::Lun::~Lun()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::Lun::tau
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    tmp<Field<scalar> > g0
+    (
+        0.6
+      / max
+        (
+            1.0 - pow(alpha/alphaPacked_,1.0/3.0),
+            max(eps_*(1.0 - alpha), SMALL)
+        )
+    );
+
+    tmp<Field<scalar> > gT(uSqr/3.0);
+
+    return alpha*rho*(1.0 + alpha*(1.0 + e_)*g0)*gT;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::Lun::dTaudTheta
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    notImplemented
+    (
+        "Foam::scalar Foam::ParticleStressModels::Lun::dTau_dTheta"
+        "(const Field<scalar>&, const Field<scalar>&, const Field<scalar>&) "
+        "const"
+    );
+
+    return tmp<Field<scalar> >(NULL);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.H b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.H
new file mode 100644
index 0000000000..ace89618fb
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/Lun/Lun.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 Lun_H
+#define Lun_H
+
+#include "ParticleStressModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class Lun Declaration
+\*---------------------------------------------------------------------------*/
+
+class Lun
+:
+    public ParticleStressModel
+{
+    // Private data
+
+        //- Coefficient of restitution
+        scalar e_;
+
+        //- Smallest allowable difference from the packed volume fraction
+        scalar eps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Lun");
+
+
+    //- Constructors
+
+        //- Construct from components
+        Lun(const dictionary& dict);
+
+        //- Construct copy
+        Lun(const Lun& hc);
+
+        //- Clone
+        virtual autoPtr<ParticleStressModel> clone() const
+        {
+            return autoPtr<ParticleStressModel>
+            (
+                new Lun(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~Lun();
+
+
+    //- Member Functions
+
+        //- Collision stress
+        tmp<Field<scalar> > tau
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+
+        //- Collision stress derivaive w.r.t. the volume fraction
+        tmp<Field<scalar> > dTaudTheta
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace ParticleStressModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.C b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.C
new file mode 100644
index 0000000000..d95d42c4a1
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "ParticleStressModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(ParticleStressModel, 0);
+    defineRunTimeSelectionTable(ParticleStressModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModel::ParticleStressModel
+(
+    const dictionary& dict
+)
+:
+    alphaPacked_(readScalar(dict.lookup("alphaPacked")))
+{
+}
+
+
+Foam::ParticleStressModel::ParticleStressModel
+(
+    const ParticleStressModel& cm
+)
+:
+    alphaPacked_(cm.alphaPacked_)
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::ParticleStressModel> Foam::ParticleStressModel::New
+(
+    const dictionary& dict
+)
+{
+    word modelType(dict.lookup("type"));
+
+    Info<< "Selecting particle stress model " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "ParticleStressModel::New"
+            "("
+                "const dictionary&"
+            ")"
+        )   << "Unknown particle stress model type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid particle stress model types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<ParticleStressModel>(cstrIter()(dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModel::~ParticleStressModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::ParticleStressModel::alphaPacked() const
+{
+    return alphaPacked_;
+}
+
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar> >
+Foam::ParticleStressModel::tau
+(
+    const FieldField<Field, scalar>& alpha,
+    const FieldField<Field, scalar>& rho,
+    const FieldField<Field, scalar>& uRms
+) const
+{
+    tmp<FieldField<Field, scalar> > value
+    (
+        new FieldField<Field, scalar>(alpha.size())
+    );
+
+    forAll(alpha, i)
+    {
+        value->set(i, tau(alpha[i], rho[i], uRms[i]));
+    }
+
+    return value;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.H b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.H
new file mode 100644
index 0000000000..bea5dade31
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/ParticleStressModel/ParticleStressModel.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 ParticleStressModel_H
+#define ParticleStressModel_H
+
+#include "fvCFD.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class ParticleStressModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class ParticleStressModel
+{
+private:
+
+    //- Private member functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const ParticleStressModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Close pack volume fraction
+        scalar alphaPacked_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("particleStressModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        ParticleStressModel,
+        dictionary,
+        (const dictionary& dict),
+        (dict)
+    );
+
+
+    //- Constructors
+
+        //- Construct from components
+        ParticleStressModel(const dictionary& dict);
+
+        //- Construct a copy
+        ParticleStressModel(const ParticleStressModel& sm);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleStressModel> clone() const = 0;
+
+
+    //- Selector
+    static autoPtr<ParticleStressModel> New
+    (
+        const dictionary& dict
+    );
+
+
+    //- Destructor
+    virtual ~ParticleStressModel();
+
+
+    //- Member Functions
+
+        //- Access max volume fraction
+        scalar alphaPacked() const;
+
+        //- Collision stress
+        virtual tmp<Field<scalar> > tau
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const = 0;
+
+        //- Collision stress derivaive w.r.t. the volume fraction
+        virtual tmp<Field<scalar> > dTaudTheta
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const = 0;
+
+        //- Collision stress using FieldFields
+        tmp<FieldField<Field, scalar> > tau
+        (
+            const FieldField<Field, scalar>& alpha,
+            const FieldField<Field, scalar>& rho,
+            const FieldField<Field, scalar>& uRms
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.C b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.C
new file mode 100644
index 0000000000..2ac8811fca
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.C
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "exponential.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+    defineTypeNameAndDebug(exponential, 0);
+
+    addToRunTimeSelectionTable
+    (
+        ParticleStressModel,
+        exponential,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::exponential::exponential
+(
+    const dictionary& dict
+)
+:
+    ParticleStressModel(dict),
+    preExp_(readScalar(dict.lookup("preExp"))),
+    expMax_(readScalar(dict.lookup("expMax"))),
+    g0_(readScalar(dict.lookup("g0")))
+{}
+
+
+Foam::ParticleStressModels::exponential::exponential
+(
+    const exponential& hc
+)
+:
+    ParticleStressModel(hc),
+    preExp_(hc.preExp_),
+    expMax_(hc.expMax_),
+    g0_(hc.g0_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::ParticleStressModels::exponential::~exponential()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::exponential::tau
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    return dTaudTheta(alpha, rho, uSqr)/preExp_;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::ParticleStressModels::exponential::dTaudTheta
+(
+    const Field<scalar>& alpha,
+    const Field<scalar>& rho,
+    const Field<scalar>& uSqr
+) const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preExp_*(alpha - alphaPacked_)),
+            expMax_
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.H b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.H
new file mode 100644
index 0000000000..8eb29a0cb6
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/ParticleStressModels/exponential/exponential.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 exponential_H
+#define exponential_H
+
+#include "ParticleStressModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ParticleStressModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class exponential Declaration
+\*---------------------------------------------------------------------------*/
+
+class exponential
+:
+    public ParticleStressModel
+{
+    // Private data
+
+        //- Pre-exponential factor
+        scalar preExp_;
+
+        //- Maximum limit of the exponential
+        scalar expMax_;
+
+        //- Front coefficient
+        scalar g0_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("exponential");
+
+
+    //- Constructors
+
+        //- Construct from components
+        exponential(const dictionary& dict);
+
+        //- Construct copy
+        exponential(const exponential& hc);
+
+        //- Clone
+        virtual autoPtr<ParticleStressModel> clone() const
+        {
+            return autoPtr<ParticleStressModel>
+            (
+                new exponential(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~exponential();
+
+
+    //- Member Functions
+
+        //- Collision stress
+        tmp<Field<scalar> > tau
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+
+        //- Collision stress derivaive w.r.t. the volume fraction
+        tmp<Field<scalar> > dTaudTheta
+        (
+            const Field<scalar>& alpha,
+            const Field<scalar>& rho,
+            const Field<scalar>& uRms
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace ParticleStressModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.C b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.C
new file mode 100644
index 0000000000..bc4c7d5397
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.C
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(TimeScaleModel, 0);
+    defineRunTimeSelectionTable(TimeScaleModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModel::TimeScaleModel
+(
+    const dictionary& dict
+)
+:
+    alphaPacked_(readScalar(dict.lookup("alphaPacked"))),
+    e_(readScalar(dict.lookup("e")))
+{
+}
+
+
+Foam::TimeScaleModel::TimeScaleModel
+(
+    const TimeScaleModel& cm
+)
+:
+    alphaPacked_(cm.alphaPacked_),
+    e_(cm.e_)
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::TimeScaleModel> Foam::TimeScaleModel::New
+(
+    const dictionary& dict
+)
+{
+    word modelType(dict.lookup("type"));
+
+    Info<< "Selecting time scale model " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "TimeScaleModel::New"
+            "("
+                "const dictionary&"
+            ")"
+        )   << "Unknown time scale model type " << modelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid time scale model types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<TimeScaleModel>(cstrIter()(dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModel::~TimeScaleModel()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.H b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.H
new file mode 100644
index 0000000000..b103d528cf
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/TimeScaleModel/TimeScaleModel.H
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 TimeScaleModel_H
+#define TimeScaleModel_H
+
+#include "fvCFD.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class TimeScaleModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class TimeScaleModel
+{
+private:
+
+    //- Private member functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const TimeScaleModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Close pack volume fraction
+        scalar alphaPacked_;
+
+        //- Coefficient of restitution
+        scalar e_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("timeScaleModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        TimeScaleModel,
+        dictionary,
+        (const dictionary& dict),
+        (dict)
+    );
+
+
+    //- Constructors
+
+        //- Construct from components
+        TimeScaleModel(const dictionary& dict);
+
+        //- Construct a copy
+        TimeScaleModel(const TimeScaleModel& sm);
+
+        //- Construct and return a clone
+        virtual autoPtr<TimeScaleModel> clone() const = 0;
+
+
+    //- Selector
+    static autoPtr<TimeScaleModel> New
+    (
+        const dictionary& dict
+    );
+
+
+    //- Destructor
+    virtual ~TimeScaleModel();
+
+
+    //- Member Functions
+
+        //- Time scale
+        virtual tmp<FieldField<Field, scalar> > oneByTau
+        (
+            const FieldField<Field, scalar>& alpha,
+            const FieldField<Field, scalar>& r32,
+            const FieldField<Field, scalar>& uSqr,
+            const FieldField<Field, scalar>& f
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.C b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.C
new file mode 100644
index 0000000000..c3d35c37fe
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "equilibrium.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+    defineTypeNameAndDebug(equilibrium, 0);
+
+    addToRunTimeSelectionTable
+    (
+        TimeScaleModel,
+        equilibrium,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::equilibrium::equilibrium
+(
+    const dictionary& dict
+)
+:
+    TimeScaleModel(dict)
+{}
+
+
+Foam::TimeScaleModels::equilibrium::equilibrium
+(
+    const equilibrium& hc
+)
+:
+    TimeScaleModel(hc)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::equilibrium::~equilibrium()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar> >
+Foam::TimeScaleModels::equilibrium::oneByTau
+(
+    const FieldField<Field, scalar>& alpha,
+    const FieldField<Field, scalar>& r32,
+    const FieldField<Field, scalar>& uSqr,
+    const FieldField<Field, scalar>& f
+) const
+{
+    static const scalar a =
+        16.0/sqrt(3.0*constant::mathematical::pi)
+       *0.25*(1.0 - e_*e_);
+
+    return 
+        a
+       *alpha*sqrt(max(uSqr, 0.0))/max(r32, SMALL)
+       *alphaPacked_/max(alphaPacked_ - alpha, SMALL);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.H b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.H
new file mode 100644
index 0000000000..0fa635f935
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/equilibrium/equilibrium.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 equilibrium_H
+#define equilibrium_H
+
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class equilibrium Declaration
+\*---------------------------------------------------------------------------*/
+
+class equilibrium
+:
+    public TimeScaleModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("equilibrium");
+
+
+    //- Constructors
+
+        //- Construct from components
+        equilibrium(const dictionary& dict);
+
+        //- Construct a copy
+        equilibrium(const equilibrium& hc);
+
+        //- Construct and return a clone
+        virtual autoPtr<TimeScaleModel> clone() const
+        {
+            return autoPtr<TimeScaleModel>
+            (
+                new equilibrium(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~equilibrium();
+
+
+    //- Member Functions
+
+        //- Time scale
+        tmp<FieldField<Field, scalar> > oneByTau
+        (
+            const FieldField<Field, scalar>& alpha,
+            const FieldField<Field, scalar>& r32,
+            const FieldField<Field, scalar>& uSqr,
+            const FieldField<Field, scalar>& f
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace TimeScaleModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.C b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.C
new file mode 100644
index 0000000000..bb8219e2ae
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "isotropic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+    defineTypeNameAndDebug(isotropic, 0);
+
+    addToRunTimeSelectionTable
+    (
+        TimeScaleModel,
+        isotropic,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::isotropic::isotropic
+(
+    const dictionary& dict
+)
+:
+    TimeScaleModel(dict)
+{}
+
+
+Foam::TimeScaleModels::isotropic::isotropic
+(
+    const isotropic& hc
+)
+:
+    TimeScaleModel(hc)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::isotropic::~isotropic()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar> >
+Foam::TimeScaleModels::isotropic::oneByTau
+(
+    const FieldField<Field, scalar>& alpha,
+    const FieldField<Field, scalar>& r32,
+    const FieldField<Field, scalar>& uSqr,
+    const FieldField<Field, scalar>& f
+) const
+{
+    static const scalar a =
+        8.0*sqrt(2.0)/(5.0*constant::mathematical::pi)
+       *0.25*(3.0 - e_)*(1.0 + e_);
+    
+    return a*f*alphaPacked_/max(alphaPacked_ - alpha, SMALL);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.H b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.H
new file mode 100644
index 0000000000..6eb3cef0fc
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/isotropic/isotropic.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 isotropic_H
+#define isotropic_H
+
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class isotropic Declaration
+\*---------------------------------------------------------------------------*/
+
+class isotropic
+:
+    public TimeScaleModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("isotropic");
+
+
+    //- Constructors
+
+        //- Construct from components
+        isotropic(const dictionary& dict);
+
+        //- Construct a copy
+        isotropic(const isotropic& hc);
+
+        //- Construct and return a clone
+        virtual autoPtr<TimeScaleModel> clone() const
+        {
+            return autoPtr<TimeScaleModel>
+            (
+                new isotropic(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~isotropic();
+
+
+    //- Member Functions
+
+        //- Time scale
+        tmp<FieldField<Field, scalar> > oneByTau
+        (
+            const FieldField<Field, scalar>& alpha,
+            const FieldField<Field, scalar>& r32,
+            const FieldField<Field, scalar>& uSqr,
+            const FieldField<Field, scalar>& f
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace TimeScaleModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.C b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.C
new file mode 100644
index 0000000000..d7462e813a
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 "nonEquilibrium.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+    defineTypeNameAndDebug(nonEquilibrium, 0);
+
+    addToRunTimeSelectionTable
+    (
+        TimeScaleModel,
+        nonEquilibrium,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::nonEquilibrium::nonEquilibrium
+(
+    const dictionary& dict
+)
+:
+    TimeScaleModel(dict)
+{}
+
+
+Foam::TimeScaleModels::nonEquilibrium::nonEquilibrium
+(
+    const nonEquilibrium& hc
+)
+:
+    TimeScaleModel(hc)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::TimeScaleModels::nonEquilibrium::~nonEquilibrium()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar> >
+Foam::TimeScaleModels::nonEquilibrium::oneByTau
+(
+    const FieldField<Field, scalar>& alpha,
+    const FieldField<Field, scalar>& r32,
+    const FieldField<Field, scalar>& uSqr,
+    const FieldField<Field, scalar>& f
+) const
+{
+    static const scalar a =
+        8.0*sqrt(2.0)/(3.0*constant::mathematical::pi)
+       *0.25*(1.0 - e_*e_);
+
+    return a*f*alphaPacked_/max(alphaPacked_ - alpha, SMALL);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.H b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.H
new file mode 100644
index 0000000000..4cb1143e84
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/MPPIC/TimeScaleModels/nonEquilibrium/nonEquilibrium.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 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 nonEquilibrium_H
+#define nonEquilibrium_H
+
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace TimeScaleModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class nonEquilibrium Declaration
+\*---------------------------------------------------------------------------*/
+
+class nonEquilibrium
+:
+    public TimeScaleModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("nonEquilibrium");
+
+
+    //- Constructors
+
+        //- Construct from components
+        nonEquilibrium(const dictionary& dict);
+
+        //- Construct a copy
+        nonEquilibrium(const nonEquilibrium& hc);
+
+        //- Construct and return a clone
+        virtual autoPtr<TimeScaleModel> clone() const
+        {
+            return autoPtr<TimeScaleModel>
+            (
+                new nonEquilibrium(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~nonEquilibrium();
+
+
+    //- Member Functions
+
+        //- Time scale
+        tmp<FieldField<Field, scalar> > oneByTau
+        (
+            const FieldField<Field, scalar>& alpha,
+            const FieldField<Field, scalar>& r32,
+            const FieldField<Field, scalar>& uSqr,
+            const FieldField<Field, scalar>& f
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace TimeScaleModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //