diff --git a/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H b/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H
index d98e9369b8a87828a16dc020e59c12be469ba769..179ac149badfef1c97ec876a581d746f4195ce69 100644
--- a/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H
@@ -1,4 +1,4 @@
-    Info<< "Reading transportProperties\n" << endl;
+    Info<< "\nReading transportProperties\n" << endl;
 
     IOdictionary transportProperties
     (
@@ -31,7 +31,7 @@
         rhoInfValue
     );
 
-    Info<< "\nReading field U\n" << endl;
+    Info<< "Reading field U\n" << endl;
     volVectorField U
     (
         IOobject
@@ -127,7 +127,7 @@
 
     if (HdotGradHheader.headerOk())
     {
-        Info<< "\nReading field HdotGradH\n" << endl;
+        Info<< "Reading field HdotGradH" << endl;
 
         HdotGradHPtr_.reset
         (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
index 4b41508a868caf95ad2b97fe3e9b9bd979dbf328..a949b8e5a31ee3d255f9f14f350e196317181675 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
@@ -110,13 +110,11 @@ class timeVaryingMappedFixedValueFvPatchField
             label& hi
         ) const;
 
+
         //- Read boundary points and determine interpolation weights to patch
         //  faceCentres
         void readSamplePoints();
 
-        //- Find boundary data inbetween current time and interpolate
-        void checkTable();
-
         //- Do actual interpolation using current weights
         tmp<Field<Type> > interpolate(const Field<Type>&) const;
 
@@ -199,6 +197,12 @@ public:
                 return referenceCS_;
             }
 
+            //- Return startSampledValues
+            const Field<Type> startSampledValues()
+            {
+                 return startSampledValues_;
+            }
+
 
         // Mapping functions
 
@@ -216,6 +220,12 @@ public:
             );
 
 
+        // Utility functions
+
+            //- Find boundary data inbetween current time and interpolate
+            void checkTable();
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.H b/src/lagrangian/basic/InteractionLists/InteractionLists.H
index 77da111828664fa21f5b4e461ff20c7171a658b4..8382c77bd7509f900b4f30d0d45433c3b7d2a23f 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.H
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.H
@@ -287,10 +287,13 @@ public:
             wallFaceIndexAndTransformToDistribute() const;
 
             //- Return access to the referred wall faces
-            const List<referredWallFace>& referredWallFaces() const;
+            inline const List<referredWallFace>& referredWallFaces() const;
+
+            //- Return the name of the velocity field
+            inline const word& UName() const;
 
             //- Return access to the referred wall data
-            const List<vector>& referredWallData() const;
+            inline const List<vector>& referredWallData() const;
 
             //- Return access to the referred particle container
             inline const List<IDLList<ParticleType> >&
diff --git a/src/lagrangian/basic/InteractionLists/InteractionListsI.H b/src/lagrangian/basic/InteractionLists/InteractionListsI.H
index ad28f696c5b1df3b6204d092889aefa1e49ccc79..1b61942f5eb671c6f97db05efbf25527b1b7b6bb 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionListsI.H
+++ b/src/lagrangian/basic/InteractionLists/InteractionListsI.H
@@ -128,6 +128,13 @@ Foam::InteractionLists<ParticleType>::referredWallFaces() const
 }
 
 
+template<class ParticleType>
+const Foam::word& Foam::InteractionLists<ParticleType>::UName() const
+{
+    return UName_;
+}
+
+
 template<class ParticleType>
 const Foam::List<Foam::vector>&
 Foam::InteractionLists<ParticleType>::referredWallData() const
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 486c0959dddfb299e64c1c7747ebc2ff3f98a62f..7fc3cf863634c4167d4fc3dc88fd087f5ba1ad4c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -42,7 +42,7 @@ template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::preEvolve()
 {
     this->dispersion().cacheFields(true);
-    forces_.cacheFields(true);
+    forces_.cacheFields(true, interpolationSchemes_);
 }
 
 
@@ -152,7 +152,7 @@ void Foam::KinematicCloud<ParcelType>::postEvolve()
     }
 
     this->dispersion().cacheFields(false);
-    forces_.cacheFields(false);
+    forces_.cacheFields(false, interpolationSchemes_);
 
     this->postProcessing().post();
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 31c274ea38d9fd0372992665f3b7edd61d975890..fc34e3b1614498aa2c19a7f3559d23bfbaf4d943 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -162,6 +162,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Momentum source due to particle forces
     const vector FCoupled = mass*td.cloud().forces().calcCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,
@@ -173,6 +174,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 4dd608c5bb236e4fd28137a281cbb4143e80755b..84123c89f4c202aa4bd926c10e42142d78686a11 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -182,8 +182,8 @@ public:
 
         // Constructors
 
-           //- Construct from components
-           inline trackData
+            //- Construct from components
+            inline trackData
             (
                 KinematicCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 7f531a05cfd7b0598570694b43ea288eb779e676..ceef9ca8b2eb984d29fbafcf0d8b21c4c784f285 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -439,14 +439,23 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::wallImpactDistance
     const vector&
 ) const
 {
-    // Do not use a wall impact distance to allow proper multiple face
-    // collisions.  In a twisted mesh the particle can be within range
-    // of a wall but not in the cell attached to a wall face, hence
-    // miss the interaction.
+    const KinematicCloud<ParcelType>& c =
+        dynamic_cast<const KinematicCloud<ParcelType>&>(this->cloud());
 
-    return 0.0;
+    if (c.collision().controlsWallInteraction())
+    {
+        // Do not use a wall impact distance if the collision model
+        // controls wall interactions to allow proper multiple face
+        // collisions.  In a twisted mesh the particle can be within
+        // range of a wall but not in the cell attached to a wall
+        // face, hence miss the interaction.
 
-    // return 0.5*d_;
+        return 0.0;
+    }
+    else
+    {
+        return 0.5*d_;
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C
index 5d446037fe6b79a31eb33e8cd47ba965d7f68fe6..31db9460f1b98e485d6dfb6e3b58e8392c85546f 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.C
+++ b/src/lagrangian/intermediate/particleForces/particleForces.C
@@ -30,6 +30,24 @@ License
 #include "mathematicalConstants.H"
 #include "electromagneticConstants.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::particleForces::deleteFields()
+{
+    if (gradUPtr_)
+    {
+        delete gradUPtr_;
+        gradUPtr_ = NULL;
+    }
+
+    if (HdotGradHInterPtr_)
+    {
+        delete HdotGradHInterPtr_;
+        HdotGradHInterPtr_ = NULL;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::particleForces::particleForces
@@ -85,7 +103,7 @@ Foam::particleForces::particleForces(const particleForces& f)
 
 Foam::particleForces::~particleForces()
 {
-    cacheFields(false);
+    deleteFields();
 }
 
 
@@ -151,26 +169,46 @@ const Foam::word& Foam::particleForces::HdotGradHName() const
 }
 
 
-void Foam::particleForces::cacheFields(const bool store)
+void Foam::particleForces::cacheFields
+(
+    const bool store,
+    const dictionary& interpolationSchemes
+)
 {
-    if (store && pressureGradient_)
+    if (store)
     {
-        const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
-        gradUPtr_ = fvc::grad(U).ptr();
+        if (pressureGradient_)
+        {
+            const volVectorField& U =
+                mesh_.lookupObject<volVectorField>(UName_);
+
+            gradUPtr_ = fvc::grad(U).ptr();
+        }
+
+        if (paramagnetic_)
+        {
+            const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+            (
+                HdotGradHName_
+            );
+
+            HdotGradHInterPtr_ = interpolation<vector>::New
+            (
+                interpolationSchemes,
+                HdotGradH
+            ).ptr();
+        }
     }
     else
     {
-        if (gradUPtr_)
-        {
-            delete gradUPtr_;
-            gradUPtr_ = NULL;
-        }
+        deleteFields();
     }
 }
 
 
 Foam::vector Foam::particleForces::calcCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -205,6 +243,7 @@ Foam::vector Foam::particleForces::calcCoupled
 
 Foam::vector Foam::particleForces::calcNonCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -226,15 +265,12 @@ Foam::vector Foam::particleForces::calcNonCoupled
 
     if (paramagnetic_)
     {
-        const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
-        (
-            HdotGradHName_
-        );
+        const interpolation<vector>& HdotGradHInter = *HdotGradHInterPtr_;
 
         accelTot +=
             3.0*constant::electromagnetic::mu0.value()/rho
            *magneticSusceptibility_/(magneticSusceptibility_ + 3)
-           *HdotGradH[cellI];
+           *HdotGradHInter.interpolate(position, cellI);
 
         // force is:
 
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index 90c5ecafa43101deb433ddee8a3389fa04e9e41a..8c76b1f773dfbba1d90f146878e649468b36abcd 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.H
+++ b/src/lagrangian/intermediate/particleForces/particleForces.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "Switch.H"
 #include "vector.H"
 #include "volFieldsFwd.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +70,9 @@ class particleForces
         //- Velocity gradient field
         const volTensorField* gradUPtr_;
 
+        //- HdotGradH interpolator
+        const interpolation<vector>* HdotGradHInterPtr_;
+
 
         // Forces to include in particle motion evaluation
 
@@ -101,6 +105,12 @@ class particleForces
             const word HdotGradHName_;
 
 
+    // Private Member Functions
+
+        //- Delete cached carrier fields
+        void deleteFields();
+
+
 public:
 
     // Constructors
@@ -159,11 +169,16 @@ public:
        // Evaluation
 
             //- Cache carrier fields
-            void cacheFields(const bool store);
+            void cacheFields
+            (
+                const bool store,
+                const dictionary& interpolationSchemes
+            );
 
             //- Calculate action/reaction forces between carrier and particles
             vector calcCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,
@@ -176,6 +191,7 @@ public:
             //- Calculate external forces applied to the particles
             vector calcNonCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.H
index 3420f8d9e5b2922170c051e9371db1616aac71b0..cbce0ee43acc85e130ea867c751128a69a7be97f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.H
@@ -141,6 +141,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const = 0;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const = 0;
+
         // Collision function
         virtual void collide() = 0;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.C
index f99edcd47f47be531341441632f4cc6c6b032635..993fe53b7041abab78c8238ea2e013123dfa5bc1 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.C
@@ -61,6 +61,13 @@ bool Foam::NoCollision<CloudType>::active() const
 }
 
 
+template<class CloudType>
+bool Foam::NoCollision<CloudType>::controlsWallInteraction() const
+{
+    return false;
+}
+
+
 template<class CloudType>
 void Foam::NoCollision<CloudType>::collide()
 {}
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H
index 5806806892e1e40186a2c6b4974498df2a59bc79..e0c0ffbc6b5b9400fd4f87aae33804fd47484102 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H
@@ -82,6 +82,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const;
+
         // Collision function
         virtual void collide();
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
index 967acf06401671c9cb4c47dd9e560f9d5ccc01d3..e1e1587c95dc9afae113bd90d7129d9f2eb85111 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
@@ -175,13 +175,20 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
     const labelListList directWallFaces = il_.dwfil();
 
+    const labelList& patchID = mesh.boundaryMesh().patchID();
+
+    const volVectorField& U = mesh.lookupObject<volVectorField>(il_.UName());
+
     // Storage for the wall interaction sites
-    DynamicList<point> flatSites;
+    DynamicList<point> flatSitePoints;
     DynamicList<scalar> flatSiteExclusionDistancesSqr;
-    DynamicList<point> otherSites;
+    DynamicList<WallSiteData<vector> > flatSiteData;
+    DynamicList<point> otherSitePoints;
     DynamicList<scalar> otherSiteDistances;
-    DynamicList<point> sharpSites;
+    DynamicList<WallSiteData<vector> > otherSiteData;
+    DynamicList<point> sharpSitePoints;
     DynamicList<scalar> sharpSiteExclusionDistancesSqr;
+    DynamicList<WallSiteData<vector> > sharpSiteData;
 
     forAll(dil, realCellI)
     {
@@ -191,19 +198,22 @@ void Foam::PairCollision<CloudType>::wallInteraction()
         // Loop over all Parcels in cell
         forAll(cellOccupancy_[realCellI], cellParticleI)
         {
-            flatSites.clear();
+            flatSitePoints.clear();
             flatSiteExclusionDistancesSqr.clear();
-            otherSites.clear();
+            flatSiteData.clear();
+            otherSitePoints.clear();
             otherSiteDistances.clear();
-            sharpSites.clear();
+            otherSiteData.clear();
+            sharpSitePoints.clear();
             sharpSiteExclusionDistancesSqr.clear();
+            sharpSiteData.clear();
 
             typename CloudType::parcelType& p =
-            *cellOccupancy_[realCellI][cellParticleI];
+                *cellOccupancy_[realCellI][cellParticleI];
 
             const point& pos = p.position();
 
-            scalar r = p.d()/2;
+            scalar r = wallModel_->pREff(p);
 
             // real wallFace interactions
 
@@ -229,6 +239,18 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     scalar normalAlignment = normal & pW/mag(pW);
 
+                    // Find the patchIndex and wallData for WallSiteData object
+                    label patchI = patchID[realFaceI - mesh.nInternalFaces()];
+
+                    label patchFaceI =
+                        realFaceI - mesh.boundaryMesh()[patchI].start();
+
+                    WallSiteData<vector> wSD
+                    (
+                        patchI,
+                        U.boundaryField()[patchI][patchFaceI]
+                    );
+
                     if (normalAlignment > cosPhiMinFlatWall)
                     {
                         // Guard against a flat interaction being
@@ -239,25 +261,29 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         (
                             !duplicatePointInList
                             (
-                                flatSites,
+                                flatSitePoints,
                                 nearPt,
                                 sqr(r*flatWallDuplicateExclusion)
                             )
                         )
                         {
-                            flatSites.append(nearPt);
+                            flatSitePoints.append(nearPt);
 
                             flatSiteExclusionDistancesSqr.append
                             (
                                 sqr(r) - sqr(nearest.distance())
                             );
+
+                            flatSiteData.append(wSD);
                         }
                     }
                     else
                     {
-                        otherSites.append(nearPt);
+                        otherSitePoints.append(nearPt);
 
                         otherSiteDistances.append(nearest.distance());
+
+                        otherSiteData.append(wSD);
                     }
                 }
             }
@@ -269,8 +295,10 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
             forAll(cellRefWallFaces, rWFI)
             {
+                label refWallFaceI = cellRefWallFaces[rWFI];
+
                 const referredWallFace& rwf =
-                il_.referredWallFaces()[cellRefWallFaces[rWFI]];
+                    il_.referredWallFaces()[refWallFaceI];
 
                 const pointField& pts = rwf.points();
 
@@ -288,6 +316,14 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     scalar normalAlignment = normal & pW/mag(pW);
 
+                    // Find the patchIndex and wallData for WallSiteData object
+
+                    WallSiteData<vector> wSD
+                    (
+                        rwf.patchIndex(),
+                        il_.referredWallData()[refWallFaceI]
+                    );
+
                     if (normalAlignment > cosPhiMinFlatWall)
                     {
                         // Guard against a flat interaction being
@@ -298,25 +334,29 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         (
                             !duplicatePointInList
                             (
-                                flatSites,
+                                flatSitePoints,
                                 nearPt,
                                 sqr(r*flatWallDuplicateExclusion)
                             )
                         )
                         {
-                            flatSites.append(nearPt);
+                            flatSitePoints.append(nearPt);
 
                             flatSiteExclusionDistancesSqr.append
                             (
                                 sqr(r) - sqr(nearest.distance())
                             );
+
+                            flatSiteData.append(wSD);
                         }
                     }
                     else
                     {
-                        otherSites.append(nearPt);
+                        otherSitePoints.append(nearPt);
 
                         otherSiteDistances.append(nearest.distance());
+
+                        otherSiteData.append(wSD);
                     }
                 }
             }
@@ -338,13 +378,13 @@ void Foam::PairCollision<CloudType>::wallInteraction()
             {
                 label orderedIndex = sortedOtherSiteIndices[siteI];
 
-                const point& otherPt = otherSites[orderedIndex];
+                const point& otherPt = otherSitePoints[orderedIndex];
 
                 if
                 (
                     !duplicatePointInList
                     (
-                        flatSites,
+                        flatSitePoints,
                         otherPt,
                         flatSiteExclusionDistancesSqr
                     )
@@ -357,23 +397,32 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                     (
                         !duplicatePointInList
                         (
-                            sharpSites,
+                            sharpSitePoints,
                             otherPt,
                             sharpSiteExclusionDistancesSqr
                         )
                     )
                     {
-                        sharpSites.append(otherPt);
+                        sharpSitePoints.append(otherPt);
 
                         sharpSiteExclusionDistancesSqr.append
                         (
                             sqr(r) - sqr(otherSiteDistances[orderedIndex])
                         );
+
+                        sharpSiteData.append(otherSiteData[orderedIndex]);
                     }
                 }
             }
 
-            evaluateWall(p, flatSites, sharpSites);
+            evaluateWall
+            (
+                p,
+                flatSitePoints,
+                flatSiteData,
+                sharpSitePoints,
+                sharpSiteData
+            );
         }
     }
 }
@@ -463,11 +512,20 @@ template<class CloudType>
 void Foam::PairCollision<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
-    const List<point>& flatSites,
-    const List<point>& sharpSites
+    const List<point>& flatSitePoints,
+    const List<WallSiteData<vector> >& flatSiteData,
+    const List<point>& sharpSitePoints,
+    const List<WallSiteData<vector> >& sharpSiteData
 ) const
 {
-    wallModel_->evaluateWall(p, flatSites, sharpSites);
+    wallModel_->evaluateWall
+    (
+        p,
+        flatSitePoints,
+        flatSiteData,
+        sharpSitePoints,
+        sharpSiteData
+    );
 }
 
 
@@ -502,7 +560,8 @@ Foam::PairCollision<CloudType>::PairCollision
     (
         owner.mesh(),
         readScalar(this->coeffDict().lookup("maxInteractionDistance")),
-        Switch(this->coeffDict().lookup("writeReferredParticleCloud"))
+        Switch(this->coeffDict().lookup("writeReferredParticleCloud")),
+        this->coeffDict().lookupOrDefault("UName", word("U"))
     )
 {}
 
@@ -552,6 +611,13 @@ bool Foam::PairCollision<CloudType>::active() const
 }
 
 
+template<class CloudType>
+bool Foam::PairCollision<CloudType>::controlsWallInteraction() const
+{
+    return true;
+}
+
+
 template<class CloudType>
 void Foam::PairCollision<CloudType>::collide()
 {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.H
index 8582db464e25d6d2df5395344e85694b5b49895e..8699d477cdd225f52000a58bd283e8a175a4defd 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.H
@@ -36,6 +36,7 @@ SourceFiles
 
 #include "CollisionModel.H"
 #include "InteractionLists.H"
+#include "WallSiteData.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -137,8 +138,10 @@ class PairCollision
         void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
-            const List<point>& sharpSites
+            const List<point>& flatSitePoints,
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const;
 
 
@@ -171,6 +174,10 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const;
 
+        //- Indicates whether model determines wall collisions or not,
+        //  used to determine what value to use for wallImpactDistance
+        virtual bool controlsWallInteraction() const;
+
         // Collision function
         virtual void collide();
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairModel/PairModel.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairModel/PairModel.H
index 1b6d1e56be43bafb0a9b6f195e86e3ce4bde155a..b5b36fdb68b8e5542ea1d93d6c15e48fbf48f09e 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairModel/PairModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairModel/PairModel.H
@@ -54,7 +54,7 @@ class PairModel
 {
     // Private data
 
-        //- The cloud dictionary
+        //- The CollisionModel dictionary
         const dictionary& dict_;
 
         //- Reference to the owner cloud class
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
index bbefa136cacb2661e18efa1a8ac3406c26ecd01c..6c10eff509f8f3428ae62ee86183a6df5339f9a8 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
@@ -45,13 +45,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::findMinMaxProperties
 
         // Finding minimum diameter to avoid excessive arithmetic
 
-        RMin = min(p.d(), RMin);
+        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+
+        RMin = min(dEff, RMin);
 
         rhoMax = max(p.rho(), rhoMax);
 
         UMagMax = max
         (
-            mag(p.U()) + mag(p.omega())*p.d()/2,
+            mag(p.U()) + mag(p.omega())*dEff/2,
             UMagMax
         );
     }
@@ -91,7 +93,8 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
         (
             this->coeffDict().lookup("collisionResolutionSteps")
         )
-    )
+    ),
+    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
 {
     scalar nu = this->owner().constProps().poissonsRatio();
 
@@ -155,7 +158,11 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 {
     vector r_AB = (pA.position() - pB.position());
 
-    scalar normalOverlapMag = 0.5*(pA.d() + pB.d()) - mag(r_AB);
+    scalar dAEff = pA.d()*cbrt(pA.nParticle()*volumeFactor_);
+
+    scalar dBEff = pB.d()*cbrt(pB.nParticle()*volumeFactor_);
+
+    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - mag(r_AB);
 
     if (normalOverlapMag > 0)
     {
@@ -166,7 +173,7 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         vector U_AB = pA.U() - pB.U();
 
         // Effective radius
-        scalar R = 0.5*pA.d()*pB.d()/(pA.d() + pB.d());
+        scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
 
         // Effective mass
         scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
@@ -185,8 +192,8 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 
         vector USlip_AB =
             U_AB - (U_AB & rHat_AB)*rHat_AB
-          + (pA.omega() ^ (pA.d()/2*-rHat_AB))
-          - (pB.omega() ^ (pB.d()/2*rHat_AB));
+          + (pA.omega() ^ (dAEff/2*-rHat_AB))
+          - (pB.omega() ^ (dBEff/2*rHat_AB));
 
         scalar deltaT = this->owner().mesh().time().deltaTValue();
 
@@ -215,7 +222,7 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         {
             scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;
 
-            scalar& etaT = etaN;
+            scalar etaT = etaN;
 
             // Tangential force
             vector fT_AB;
@@ -233,16 +240,16 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             else
             {
                 fT_AB =
-                -kT*tangentialOverlapMag
-               *tangentialOverlap_AB/tangentialOverlapMag
-              - etaT*USlip_AB;
+                    -kT*tangentialOverlapMag
+                   *tangentialOverlap_AB/tangentialOverlapMag
+                  - etaT*USlip_AB;
             }
 
             pA.f() += fT_AB;
             pB.f() += -fT_AB;
 
-            pA.torque() += (pA.d()/2*-rHat_AB) ^ fT_AB;
-            pB.torque() += (pB.d()/2*rHat_AB) ^ -fT_AB;
+            pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
+            pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
         }
     }
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H
index 819e9a835bf72a7f1d4a530d697ad04a5d03b1c8..18b30753373fc338401f5a656ccb776ed40e60ea 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H
@@ -74,6 +74,24 @@ class PairSpringSliderDashpot
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
 
+        //- Volume factor for determining the equivalent size of a
+        //  parcel where nParticles is not 1.  The equivalent size of
+        //  the parcel is
+        //      parcelEquivVolume = volumeFactor*nParticles*p.volume()
+        //  so
+        //      parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
+        //  + When volumeFactor = 1, the particles are compressed
+        //    together so that the equivalent volume of the parcel is
+        //    the sum of the constituent particles
+        //  + When volumeFactor = 3*sqrt(2)/pi, the particles are
+        //    close packed, but uncompressed.
+        //  + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
+        //    grouped.
+        // 3*sqrt(2)/pi = 1.350474 is the volume factor for close
+        // packing, i.e pi/(3*sqrt(2)) is the maximum close packing
+        // factor
+        scalar volumeFactor_;
+
 
     // Private Member Functions
 
@@ -104,6 +122,12 @@ public:
 
     // Member Functions
 
+        //- Return the volumeFactor
+        inline scalar volumeFactor() const
+        {
+            return volumeFactor_;
+        }
+
         //- Whether the PairModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H
index 7088c4df05116381982734d5155c62b37bb78700..31fde4e2991f3293e29e124f3e6f53453a7800da 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H
@@ -54,7 +54,7 @@ class WallModel
 {
     // Private data
 
-        //- The cloud dictionary
+        //- The CollisionModel dictionary
         const dictionary& dict_;
 
         //- Reference to the owner cloud class
@@ -120,6 +120,9 @@ public:
 
     // Member Functions
 
+        //- Return the effective radius for a particle for the model
+        virtual scalar pREff(const typename CloudType::parcelType& p) const = 0;
+
         //- Whether the WallModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const = 0;
@@ -133,8 +136,10 @@ public:
         virtual void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
-            const List<point>& sharpSites
+            const List<point>& flatSitePoints,
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const = 0;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
index ba22ef38429ea87d3e9447423b3d2fb7a7d6dcb6..60c4218a7509d3786aec5f89440fe854d3ff16ea 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
@@ -44,13 +44,16 @@ void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
         const typename CloudType::parcelType& p = iter();
 
         // Finding minimum diameter to avoid excessive arithmetic
-        rMin = min(p.d(), rMin);
+
+        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+
+        rMin = min(dEff, rMin);
 
         rhoMax = max(p.rho(), rhoMax);
 
         UMagMax = max
         (
-            mag(p.U()) + mag(p.omega())*p.d()/2,
+            mag(p.U()) + mag(p.omega())*dEff/2,
             UMagMax
         );
     }
@@ -66,25 +69,75 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
     const point& site,
+    const WallSiteData<vector>& data,
     scalar pNu,
     scalar pE,
+    scalar pREff,
     scalar Estar,
-    scalar kN
+    scalar kN,
+    scalar Gstar
 ) const
 {
     vector r_PW = p.position() - site;
 
-    scalar normalOverlapMag = p.d()/2 - mag(r_PW);
+    vector U_PW = p.U() - data.wallData();
+
+    scalar normalOverlapMag = pREff - mag(r_PW);
 
     vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
 
     scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
 
     vector fN_PW =
-    rHat_PW
-    *(kN*pow(normalOverlapMag, b_) - etaN*(p.U() & rHat_PW));
+        rHat_PW
+       *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
 
     p.f() += fN_PW;
+
+    vector USlip_PW =
+        U_PW - (U_PW & rHat_PW)*rHat_PW
+      + (p.omega() ^ (pREff*-rHat_PW));
+
+    scalar deltaT = this->owner().mesh().time().deltaTValue();
+
+    // For remembering previous overlap
+    // vector deltaTangentialOverlap_PW = USlip_PW * deltaT;
+    // tangentialOverlap_PW += deltaTangentialOverlap_PW;
+
+    vector tangentialOverlap_PW = USlip_PW * deltaT;
+
+    scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
+
+    if (tangentialOverlapMag > VSMALL)
+    {
+        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
+
+        scalar etaT = etaN;
+
+        // Tangential force
+        vector fT_PW;
+
+        if (kT*tangentialOverlapMag > mu_*mag(fN_PW))
+        {
+            // Tangential force greater than sliding friction,
+            // particle slips
+
+            fT_PW = -mu_*mag(fN_PW)*USlip_PW/mag(USlip_PW);
+
+            // tangentialOverlap_PW = vector::zero;
+        }
+        else
+        {
+            fT_PW =
+                -kT*tangentialOverlapMag
+               *tangentialOverlap_PW/tangentialOverlapMag
+              - etaT*USlip_PW;
+        }
+
+        p.f() += fT_PW;
+
+        p.torque() += (pREff*-rHat_PW) ^ fT_PW;
+    }
 }
 
 
@@ -109,7 +162,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
         (
             this->coeffDict().lookup("collisionResolutionSteps")
         )
-    )
+    ),
+    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
 {}
 
 
@@ -122,6 +176,15 @@ Foam::WallSpringSliderDashpot<CloudType>::~WallSpringSliderDashpot()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+Foam::scalar Foam::WallSpringSliderDashpot<CloudType>::pREff
+(
+    const typename CloudType::parcelType& p
+) const
+{
+    return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
+}
+
 template<class CloudType>
 bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
 {
@@ -164,28 +227,56 @@ template<class CloudType>
 void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
-    const List<point>& flatSites,
-    const List<point>& sharpSites
+    const List<point>& flatSitePoints,
+    const List<WallSiteData<vector> >& flatSiteData,
+    const List<point>& sharpSitePoints,
+    const List<WallSiteData<vector> >& sharpSiteData
 ) const
 {
     scalar pNu = this->owner().constProps().poissonsRatio();
 
     scalar pE = this->owner().constProps().youngsModulus();
 
+    scalar pREff = this->pREff(p);
+
     scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
 
-    scalar kN = (4.0/3.0)*sqrt(p.d()/2)*Estar;
+    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
 
-    forAll(flatSites, siteI)
+    scalar GStar = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu_ - sqr(nu_))/E_));
+
+    forAll(flatSitePoints, siteI)
     {
-        evaluateWall(p, flatSites[siteI], pNu, pE, Estar, kN);
+        evaluateWall
+        (
+            p,
+            flatSitePoints[siteI],
+            flatSiteData[siteI],
+            pNu,
+            pE,
+            pREff,
+            Estar,
+            kN,
+            GStar
+        );
     }
 
-    forAll(sharpSites, siteI)
+    forAll(sharpSitePoints, siteI)
     {
         // Treating sharp sites like flat sites
 
-        evaluateWall(p, sharpSites[siteI], pNu, pE, Estar, kN);
+        evaluateWall
+        (
+            p,
+            sharpSitePoints[siteI],
+            sharpSiteData[siteI],
+            pNu,
+            pE,
+            pREff,
+            Estar,
+            kN,
+            GStar
+        );
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
index fcc83e22e6cb6278b13ba01960cfce9883a6fae9..2aefa7d7e7d70accf278ba26b29262146dd8078c 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
@@ -69,6 +69,24 @@ class WallSpringSliderDashpot
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
 
+        //- Volume factor for determining the equivalent size of a
+        //  parcel where nParticles is not 1.  The equivalent size of
+        //  the parcel is
+        //      parcelEquivVolume = volumeFactor*nParticles*p.volume()
+        //  so
+        //      parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
+        //  + When volumeFactor = 1, the particles are compressed
+        //    together so that the equivalent volume of the parcel is
+        //    the sum of the constituent particles
+        //  + When volumeFactor = 3*sqrt(2)/pi, the particles are
+        //    close packed, but uncompressed.
+        //  + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
+        //    grouped.
+        // 3*sqrt(2)/pi = 1.350474 is the volume factor for close
+        // packing, i.e pi/(3*sqrt(2)) is the maximum close packing
+        // factor
+        scalar volumeFactor_;
+
 
     // Private Member Functions
 
@@ -86,10 +104,13 @@ class WallSpringSliderDashpot
         (
             typename CloudType::parcelType& p,
             const point& site,
+            const WallSiteData<vector>& data,
             scalar pNu,
             scalar pE,
+            scalar pREff,
             scalar Estar,
-            scalar kN
+            scalar kN,
+            scalar Gstar
         ) const;
 
 
@@ -111,6 +132,15 @@ public:
 
     // Member Functions
 
+        //- Return the volumeFactor
+        inline scalar volumeFactor() const
+        {
+            return volumeFactor_;
+        }
+
+        //- Return the effective radius for a particle for the model
+        virtual scalar pREff(const typename CloudType::parcelType& p) const;
+
         //- Whether the WallModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;
@@ -124,8 +154,10 @@ public:
         virtual void evaluateWall
         (
             typename CloudType::parcelType& p,
-            const List<point>& flatSites,
-            const List<point>& sharpSites
+            const List<point>& flatSitePoints,
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.C
new file mode 100644
index 0000000000000000000000000000000000000000..a9aa774c647cca91bf61538b33ac8c3c3f6e701a
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.C
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "WallSiteData.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::WallSiteData<Type>::WallSiteData()
+:
+    patchI_(),
+    wallData_()
+{}
+
+
+template<class Type>
+Foam::WallSiteData<Type>::WallSiteData
+(
+    label patchI,
+    const Type& wallData
+)
+:
+    patchI_(patchI),
+    wallData_(wallData)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::WallSiteData<Type>::~WallSiteData()
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::WallSiteData<Type>::operator==
+(
+    const WallSiteData<Type>& rhs
+) const
+{
+    return patchI_ == rhs.patchI_ && wallData_ == rhs.wallData_;
+}
+
+
+template<class Type>
+bool Foam::WallSiteData<Type>::operator!=
+(
+    const WallSiteData<Type>& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    WallSiteData<Type>& wIS
+)
+{
+    is  >> wIS.patchI_ >> wIS.wallData_;
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>"
+        "(Foam::Istream&, Foam::WallSiteData<Type>&)"
+    );
+
+    return is;
+}
+
+
+template<class Type>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const WallSiteData<Type>& wIS
+)
+{
+    os  << wIS.patchI_ << token::SPACE << wIS.wallData_;
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Ostream&, const WallSiteData<Type>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.H
new file mode 100644
index 0000000000000000000000000000000000000000..a179d4925c97c8399851a56039d5538cd3103528
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteData.H
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::WallSiteData
+
+Description
+    Stores the patch ID and templated data to represent a collision
+    with a wall to be passed to the wall model.
+
+SourceFiles
+    WallSiteDataI.H
+    WallSiteData.C
+    WallSiteDataIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallSiteData_H
+#define WallSiteData_H
+
+#include "label.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class WallSiteData;
+
+template<class Type>
+Istream& operator>>(Istream&, WallSiteData<Type>&);
+
+template<class Type>
+Ostream& operator<<(Ostream&, const WallSiteData<Type>&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class WallSiteData Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class WallSiteData
+{
+    // Private data
+
+        //- Index of originating patch
+        label patchI_;
+
+        //- Wall data
+        Type wallData_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        WallSiteData();
+
+        //- Construct from components
+        WallSiteData
+        (
+            label patchI,
+            const Type& wallData
+        );
+
+
+    //- Destructor
+    ~WallSiteData();
+
+
+    // Member Functions
+
+            //- Return access to the patch index
+            inline label patchIndex() const;
+
+            //- Return non-const access to the patch index
+            inline label& patchIndex();
+
+            //- Return access to wall data
+            inline const Type& wallData() const;
+
+            //- Return non-const access to wall data
+            inline Type& wallData();
+
+
+    // Member Operators
+
+        bool operator==(const WallSiteData<Type>&) const;
+        bool operator!=(const WallSiteData<Type>&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>> <Type>
+        (Istream&, WallSiteData<Type>&);
+
+        friend Ostream& operator<< <Type>
+        (Ostream&, const WallSiteData<Type>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "WallSiteDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallSiteData.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteDataI.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteDataI.H
new file mode 100644
index 0000000000000000000000000000000000000000..b28bf3be23edd1ccd892b08f1d25eb93353d48c4
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallSiteData/WallSiteDataI.H
@@ -0,0 +1,60 @@
+#/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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>
+Foam::label Foam::WallSiteData<Type>::patchIndex() const
+{
+    return patchI_;
+}
+
+
+template<class Type>
+Foam::label& Foam::WallSiteData<Type>::patchIndex()
+{
+    return patchI_;
+}
+
+
+template<class Type>
+const Type& Foam::WallSiteData<Type>::wallData() const
+{
+    return wallData_;
+}
+
+
+template<class Type>
+Type& Foam::WallSiteData<Type>::wallData()
+{
+    return wallData_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index a232a5f0609c754717d9d545817911d299b3054f..9b44324461ba87abce6e45c553f9f52fab06b71d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -129,10 +129,11 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 
 
 template<class CloudType>
-void Foam::InjectionModel<CloudType>::findCellAtPosition
+bool Foam::InjectionModel<CloudType>::findCellAtPosition
 (
     label& cellI,
-    vector& position
+    vector& position,
+    bool errorOnNotFound
 )
 {
     const volVectorField& cellCentres = owner_.mesh().C();
@@ -176,17 +177,26 @@ void Foam::InjectionModel<CloudType>::findCellAtPosition
 
     if (procI == -1)
     {
-        FatalErrorIn
-        (
-            "Foam::InjectionModel<CloudType>::findCellAtPosition"
-            "("
-                "label&, "
-                "vector&"
-            ")"
-        )<< "Cannot find parcel injection cell. "
-         << "Parcel position = " << p0 << nl
-         << abort(FatalError);
+        if (errorOnNotFound)
+        {
+            FatalErrorIn
+            (
+                "Foam::InjectionModel<CloudType>::findCellAtPosition"
+                "("
+                    "label&, "
+                    "vector&"
+                ")"
+            )   << "Cannot find parcel injection cell. "
+                << "Parcel position = " << p0 << nl
+                << abort(FatalError);
+        }
+        else
+        {
+            return false;
+        }
     }
+
+    return true;
 }
 
 
@@ -217,7 +227,7 @@ Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
         }
         case pbFixed:
         {
-            nP = nParticlesFixed_;
+            nP = nParticleFixed_;
             break;
         }
         default:
@@ -290,7 +300,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
-    nParticlesFixed_(0.0),
+    nParticleFixed_(0.0),
     time0_(0.0),
     timeStep0_(0.0)
 {
@@ -316,7 +326,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
-    nParticlesFixed_(0.0),
+    nParticleFixed_(0.0),
     time0_(owner.db().time().value()),
     timeStep0_(0.0)
 {
@@ -340,11 +350,11 @@ Foam::InjectionModel<CloudType>::InjectionModel
     {
         parcelBasis_ = pbFixed;
 
-        Info<< "    Choosing nParticles to be a fixed value, massTotal "
+        Info<< "    Choosing nParticle to be a fixed value, massTotal "
             << "variable now does not determine anything."
             << endl;
 
-        nParticlesFixed_ = readScalar(coeffDict_.lookup("nParticles"));
+        nParticleFixed_ = readScalar(coeffDict_.lookup("nParticle"));
     }
     else
     {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index c6fb6e424bc3270da4d41f4f546a493a6a79c1c3..ad3d4b1daeddba1bd4be74d814b4f79022c9eeec 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -138,9 +138,9 @@ protected:
             //- Parcel basis enumeration
             parcelBasis parcelBasis_;
 
-            //- nParticles to assign to parcels when the 'fixed' basis
+            //- nParticle to assign to parcels when the 'fixed' basis
             //  is selected
-            scalar nParticlesFixed_;
+            scalar nParticleFixed_;
 
             //- Continuous phase time at start of injection time step [s]
             scalar time0_;
@@ -180,7 +180,12 @@ protected:
         //- Find the cell that contains the supplied position
         //  Will modify position slightly towards the owner cell centroid to
         //  ensure that it lies in a cell and not edge/face
-        virtual void findCellAtPosition(label& cellI, vector& position);
+        virtual bool findCellAtPosition
+        (
+            label& cellI,
+            vector& position,
+            bool errorOnNotFound = true
+        );
 
         //- Set number of particles to inject given parcel properties
         virtual scalar setNumberOfParticles
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
index 9b9ec19084e474baca12f0d89176795014e18ab3..f64cdc22681c5282cb6547327e116b9781fc8dbf 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
@@ -25,6 +25,8 @@ License
 
 #include "ManualInjection.H"
 #include "mathematicalConstants.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
 
 using namespace Foam::constant::mathematical;
 
@@ -100,6 +102,40 @@ Foam::ManualInjection<CloudType>::ManualInjection
         )
     )
 {
+    Switch checkAndIgnoreOutOfBounds
+    (
+        this->coeffDict().lookupOrDefault("checkAndIgnoreOutOfBounds", false)
+    );
+
+    label nRejected = 0;
+
+    if (checkAndIgnoreOutOfBounds)
+    {
+        // Dummy cell
+        label cellI = -1;
+
+        PackedBoolList keep(positions_.size(), true);
+
+        forAll(positions_, pI)
+        {
+            if (!this->findCellAtPosition(cellI, positions_[pI], false))
+            {
+                keep[pI] = false;
+
+                nRejected++;
+            }
+        }
+
+        if (nRejected > 0)
+        {
+            inplaceSubset(keep, positions_);
+            inplaceSubset(keep, diameters_);
+
+            Info<< "    " << nRejected
+                << " particles ignored, out of bounds." << endl;
+        }
+    }
+
     // Construct parcel diameters
     forAll(diameters_, i)
     {
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
index 3c63ece5dab6fec115b6dbdb473ab8a3bcf54df0..c6cbc3c11fe6665a55a711c91b34005788fe8edb 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
@@ -36,7 +36,7 @@ Description
 
 SourceFiles
     sixDoFRigidBodyMotionConstraint.C
-    dynamicFvMeshNew.C
+    sixDoFRigidBodyMotionConstraintNew.C
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
index bab292bc4ddf0df3835431c1977b5a3109788713..545b17db8dead568bb815fb0eb090b095a51bc38 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H
@@ -36,7 +36,7 @@ Description
 
 SourceFiles
     sixDoFRigidBodyMotionRestraint.C
-    dynamicFvMeshNew.C
+    sixDoFRigidBodyMotionRestraintNew.C
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/turbulenceModels/incompressible/RAS/Make/files b/src/turbulenceModels/incompressible/RAS/Make/files
index f6107822dfbc1cee80f472537f5ef4448d187c6d..9b95469a648584e4b3311e46e01119a2abc09a49 100644
--- a/src/turbulenceModels/incompressible/RAS/Make/files
+++ b/src/turbulenceModels/incompressible/RAS/Make/files
@@ -44,6 +44,8 @@ derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFv
 derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
 derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
+derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C
+derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C
 
 backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..a69a36f140c85afe0aa2fa26a988aa954e9aab86
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ABLInletEpsilonFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Ustar_(0),
+    z_(pTraits<vector>::zero),
+    z0_(0),
+    kappa_(0.41)
+{}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const ABLInletEpsilonFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    Ustar_(ptf.Ustar_),
+    z_(ptf.z_),
+    z0_(ptf.z0_),
+    kappa_(ptf.kappa_)
+{}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Ustar_(readScalar(dict.lookup("Ustar"))),
+    z_(dict.lookup("z")),
+    z0_(readScalar(dict.lookup("z0"))),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+{
+    if (mag(z_) < SMALL)
+    {
+        FatalErrorIn("ABLInletEpsilonFvPatchScalarField(dict)")
+            << "z is not correct"
+            << abort(FatalError);
+    }
+
+    z_ /= mag(z_);
+
+    evaluate();
+}
+
+
+ABLInletEpsilonFvPatchScalarField::ABLInletEpsilonFvPatchScalarField
+(
+    const ABLInletEpsilonFvPatchScalarField& fcvpvf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(fcvpvf, iF),
+    Ustar_(fcvpvf.Ustar_),
+    z_(fcvpvf.z_),
+    z0_(fcvpvf.z0_),
+    kappa_(fcvpvf.kappa_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ABLInletEpsilonFvPatchScalarField::updateCoeffs()
+{
+    const vectorField& c = patch().Cf();
+    scalarField coord = (c & z_);
+    scalarField::operator=(pow(Ustar_, 3.0)/(kappa_*(coord + z0_)));
+}
+
+
+// Write
+void ABLInletEpsilonFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("Ustar")
+        << Ustar_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z")
+        << z_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z0")
+        << z0_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa")
+        << kappa_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, ABLInletEpsilonFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..39560879aa929566299c32045848f7660db21514
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletEpsilon/ABLInletEpsilonFvPatchScalarField.H
@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    ABLInletEpsilonFvPatchScalarField
+
+Description
+    Boundary condition specifies a epsilon inlet for the atmospheric boundary
+    layer (ABL). This boundaty is to be used in conjunction with
+    ABLInletVelocity.
+
+    @verbatim
+        epsilon = Ustar^3 / (K(z + z0))
+
+    where:
+
+        Ustar is the frictional velocity
+        K is karman's constant
+        z is the verical coordinate
+        z0 is the surface roughness lenght
+
+    @endverbatim
+
+    Reference:
+    D.M. Hargreaves and N.G. Wright
+    "On the use of the k-epsilon model in commercial CFD software to model the
+     neutral atmospheric boundary layer"
+    Journal of Wind Engineering and Industrial Aerodynamics 95(2007) 355-369.
+
+SourceFiles
+    ABLInletEpsilonFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ABLInletEpsilonFvPatchScalarField_H
+#define ABLInletEpsilonFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+/*---------------------------------------------------------------------------*\
+              Class ABLInletEpsilonFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class ABLInletEpsilonFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Peak velocity magnitude
+        scalar Ustar_;
+
+        //- Direction of the z-coordinate
+        vector z_;
+
+         //- Surface roughness lenght
+        scalar z0_;
+
+         //- Von Karman constant
+        scalar kappa_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ABLInletEpsilon");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given ABLInletEpsilonFvPatchScalarField
+        //  onto a new patch
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const ABLInletEpsilonFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new ABLInletEpsilonFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        ABLInletEpsilonFvPatchScalarField
+        (
+            const ABLInletEpsilonFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new ABLInletEpsilonFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return max value
+        scalar& Ustar()
+        {
+            return Ustar_;
+        }
+
+        //- Return z direction
+        vector& z()
+        {
+            return z_;
+        }
+
+        //- Update coefficients
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C
new file mode 100644
index 0000000000000000000000000000000000000000..d8e7ae4639bcfa5f896692aa373805a85d480fac
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.C
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ABLInletVelocityFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    Ustar_(0),
+    n_(pTraits<vector>::zero),
+    z_(pTraits<vector>::zero),
+    z0_(0),
+    kappa_(0.41)
+{}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const ABLInletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    Ustar_(ptf.Ustar_),
+    n_(ptf.n_),
+    z_(ptf.z_),
+    z0_(ptf.z0_),
+    kappa_(ptf.kappa_)
+{}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    Ustar_(readScalar(dict.lookup("Ustar"))),
+    n_(dict.lookup("n")),
+    z_(dict.lookup("z")),
+    z0_(readScalar(dict.lookup("z0"))),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+{
+    if (mag(n_) < SMALL || mag(z_) < SMALL)
+    {
+        FatalErrorIn("ABLInletVelocityFvPatchVectorField(dict)")
+            << "n or z given with zero size not correct"
+            << abort(FatalError);
+    }
+
+    n_ /= mag(n_);
+    z_ /= mag(z_);
+
+    evaluate();
+}
+
+
+ABLInletVelocityFvPatchVectorField::ABLInletVelocityFvPatchVectorField
+(
+    const ABLInletVelocityFvPatchVectorField& fcvpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(fcvpvf, iF),
+    Ustar_(fcvpvf.Ustar_),
+    n_(fcvpvf.n_),
+    z_(fcvpvf.z_),
+    z0_(fcvpvf.z0_),
+    kappa_(fcvpvf.kappa_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ABLInletVelocityFvPatchVectorField::updateCoeffs()
+{
+    const vectorField& c = patch().Cf();
+    scalarField coord = (c & z_);
+    vectorField::operator=(n_*(Ustar_/kappa_)*log((coord + z0_)/z0_));
+}
+
+
+// Write
+void ABLInletVelocityFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchVectorField::write(os);
+    os.writeKeyword("Ustar")
+        << Ustar_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z0")
+        << z0_ << token::END_STATEMENT << nl;
+    os.writeKeyword("n")
+        << n_ << token::END_STATEMENT << nl;
+    os.writeKeyword("z")
+        << z_ << token::END_STATEMENT << nl;
+     os.writeKeyword("kappa")
+        << kappa_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchVectorField, ABLInletVelocityFvPatchVectorField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.H
new file mode 100644
index 0000000000000000000000000000000000000000..0e7a9167cdb4395f7346339e59666bf0f1a01243
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/ABLInletVelocity/ABLInletVelocityFvPatchVectorField.H
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    ABLInletVelocityFvPatchVectorField
+
+Description
+    Boundary condition specifies a atmospheric boundary layer (ABL)
+    velocity inlet profile given the friction velocity value,
+    flow direction n and direction of the parabolic coordinate z.
+
+    @verbatim
+        U = (Ustar/K) ln((z + z0)/z0)
+
+    where:
+
+        Ustar is the frictional velocity
+        K is karman's constant
+        z0 is the surface roughness lenght
+        z is the verical coordinate
+
+    and:
+
+        Ustar = K Uref/ln((Zref + z0)/z0)
+
+    where:
+
+        Uref is the reference velocity at Zref
+        Zref is the reference height.
+
+    @endverbatim
+
+    Reference:
+    D.M. Hargreaves and N.G. Wright
+    "On the use of the k-epsilon model in commercial CFD software to model the
+     neutral atmospheric boundary layer"
+    Journal of Wind Engineering and Industrial Aerodynamics 95(2007) 355-369.
+
+NOTE: D.M. Hargreaves and N.G. Wright recommend Gamma epsilon in the k-epsilon
+      model should be changed from 1.3 to 1.11 for consistency.
+      The roughness height (Er) is given by Er = 20 z0 following the same
+      reference
+
+SourceFiles
+    ABLInletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ABLInletVelocityFvPatchVectorField_H
+#define ABLInletVelocityFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+/*---------------------------------------------------------------------------*\
+             Class ABLInletVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class ABLInletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Friction velocity
+        scalar Ustar_;
+
+        //- Flow direction
+        vector n_;
+
+        //- Direction of the z-coordinate
+        vector z_;
+
+        //- Surface roughness lenght
+        scalar z0_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ABLInletVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        ABLInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        ABLInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given ABLInletVelocityFvPatchVectorField
+        //  onto a new patch
+        ABLInletVelocityFvPatchVectorField
+        (
+            const ABLInletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new ABLInletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        ABLInletVelocityFvPatchVectorField
+        (
+            const ABLInletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new ABLInletVelocityFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return max value
+        scalar& Ustar()
+        {
+            return Ustar_;
+        }
+
+        //- Return flow direction
+        vector& n()
+        {
+            return n_;
+        }
+
+        //- Return y direction
+        vector& z()
+        {
+            return z_;
+        }
+
+        //- Update coefficients
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
index 0397d52fae48f2550a1dde5e662fa914e4b4ac2e..1fd82c9b3c271d744bcc0036ab3ef07183c1e54f 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -10,7 +10,7 @@ License
 
     OpenFOAM is free software; you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
+    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
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H
index 25e65f4bbf2c54efd64f3509be1ee3ffa0f9575a..5597f01bc21dac55d027d3781ce405960056c1f3 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -10,7 +10,7 @@ License
 
     OpenFOAM is free software; you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
+    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
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H
index 3e6be442b59f58f03528a1273a2757d91745c004..397b66806fc5995aba08be194a66a638e9ad1846 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H
@@ -152,6 +152,21 @@ public:
 
     // Member functions
 
+        // Acces functions
+
+            // Return Ks
+             scalarField& Ks()
+            {
+                return Ks_;
+            }
+
+            // Return Cs
+            scalarField& Cs()
+            {
+                return Cs_;
+            }
+
+
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object