diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
new file mode 100644
index 0000000000000000000000000000000000000000..7a3ec75ebd100dfdae47e88d2e0c64303130ecde
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "DistortedSphereDragForce.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::DistortedSphereDragForce<CloudType>::CdRe
+(
+    const scalar Re
+) const
+{
+    if (Re > 1000.0)
+    {
+        return 0.424*Re;
+    }
+    else
+    {
+        return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::DistortedSphereDragForce
+(
+    CloudType& owner,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, false)
+{}
+
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::DistortedSphereDragForce
+(
+    const DistortedSphereDragForce<CloudType>& df
+)
+:
+    ParticleForce<CloudType>(df)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::DistortedSphereDragForce<CloudType>::~DistortedSphereDragForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::forceSuSp Foam::DistortedSphereDragForce<CloudType>::calcCoupled
+(
+    const typename CloudType::parcelType& p,
+    const scalar dt,
+    const scalar mass,
+    const scalar Re,
+    const scalar muc
+) const
+{
+    forceSuSp value(vector::zero, 0.0);
+
+    // Limit the drop distortion to y=0 (sphere) and y=1 (disk)
+    scalar y = min(max(p.y(), 0), 1);
+
+    value.Sp() = mass*0.75*muc*CdRe(Re)*(1 + 2.632*y)/(p.rho()*sqr(p.d()));
+
+    return value;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
new file mode 100644
index 0000000000000000000000000000000000000000..eb9dd62d5da22eaef21e61fe144b1d0f4ad501c8
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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::DistortedSphereDragForce
+
+Description
+    Drag model based on assumption of distorted spheres according to:
+
+    \verbatim
+        "Effects of Drop Drag and Breakup on Fuel Sprays"
+        Liu, A.B., Mather, D., Reitz, R.D.,
+        SAE Paper 930072,
+        SAE Transactions, Vol. 102, Section 3, Journal of Engines, 1993,
+        pp. 63-95
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DistortedSphereDragForce_H
+#define DistortedSphereDragForce_H
+
+#include "ParticleForce.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class DistortedSphereDragForce Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class DistortedSphereDragForce
+:
+    public ParticleForce<CloudType>
+{
+    // Private Member Functions
+
+        //- Drag coefficient multiplied by Reynolds number
+        scalar CdRe(const scalar Re) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distortedSphereDrag");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        DistortedSphereDragForce
+        (
+            CloudType& owner,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct copy
+        DistortedSphereDragForce(const DistortedSphereDragForce<CloudType>& df);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleForce<CloudType> > clone() const
+        {
+            return autoPtr<ParticleForce<CloudType> >
+            (
+                new DistortedSphereDragForce<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~DistortedSphereDragForce();
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Calculate the coupled force
+            virtual forceSuSp calcCoupled
+            (
+                const typename CloudType::parcelType& p,
+                const scalar dt,
+                const scalar mass,
+                const scalar Re,
+                const scalar muc
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DistortedSphereDragForce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index cb284c9b5d55a5964b2fea3b9e0e0321433da6f0..0b1419c5a393a041cfb3cc38af7a79e7b798061c 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -368,7 +368,7 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
 )
 {
     const scalar& TABCmu = td.cloud().breakup().TABCmu();
-    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
+    const scalar& TABtwoWeCrit = td.cloud().breakup().TABtwoWeCrit();
     const scalar& TABComega = td.cloud().breakup().TABComega();
 
     scalar r = 0.5*this->d();
@@ -385,27 +385,19 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
     {
         scalar omega = sqrt(omega2);
         scalar rhoc = this->rhoc();
-        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABWeCrit;
+        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABtwoWeCrit;
 
-        scalar y1 = this->y() - We;
-        scalar y2 = this->yDot()/omega;
+        // Initial values for y and yDot
+        scalar y0 = this->y() - We;
+        scalar yDot0 = this->yDot() + y0*rtd;
 
         // Update distortion parameters
         scalar c = cos(omega*dt);
         scalar s = sin(omega*dt);
         scalar e = exp(-rtd*dt);
-        y2 = (this->yDot() + y1*rtd)/omega;
 
-        this->y() = We + e*(y1*c + y2*s);
-        if (this->y() < 0)
-        {
-            this->y() = 0.0;
-            this->yDot() = 0.0;
-        }
-        else
-        {
-            this->yDot() = (We - this->y())*rtd + e*omega*(y2*c - y1*s);
-        }
+        this->y() = We + e*(y0*c + (yDot0/omega)*s);
+        this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
     }
     else
     {
diff --git a/src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C b/src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C
index 3f92b79a57c1a8115e42219639102678c26f4164..0430e33f91665506076f661bc559a5e8d4e103cd 100644
--- a/src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C
+++ b/src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C
@@ -45,6 +45,7 @@ License
 #include "makeReactingParcelSurfaceFilmModels.H"
 
 // Spray
+#include "DistortedSphereDragForce.H"
 #include "makeSprayParcelAtomizationModels.H"
 #include "makeSprayParcelBreakupModels.H"
 
@@ -72,6 +73,7 @@ namespace Foam
     makeReactingParcelSurfaceFilmModels(basicSprayCloud);
 
     // Spray sub-models
+    makeParticleForceModelType(DistortedSphereDragForce, basicSprayCloud);
     makeSprayParcelAtomizationModels(basicSprayCloud);
     makeSprayParcelBreakupModels(basicSprayCloud);
 };
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
index 2b9efbcea4545deb3778b5f6aa36562259f2a1f1..60047fd32fd0963a8d088b38fd1ead39216312b8 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     yDot0_(0.0),
     TABComega_(0.0),
     TABCmu_(0.0),
-    TABWeCrit_(0.0)
+    TABtwoWeCrit_(0.0)
 {}
 
 
@@ -55,7 +55,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     yDot0_(bum.yDot0_),
     TABComega_(bum.TABComega_),
     TABCmu_(bum.TABCmu_),
-    TABWeCrit_(bum.TABWeCrit_)
+    TABtwoWeCrit_(bum.TABtwoWeCrit_)
 {}
 
 
@@ -70,20 +70,19 @@ Foam::BreakupModel<CloudType>::BreakupModel
 :
     CloudSubModelBase<CloudType>(owner, dict, typeName, type),
     solveOscillationEq_(solveOscillationEq),
-    y0_(0.0),
-    yDot0_(0.0),
-    TABComega_(0.0),
-    TABCmu_(0.0),
-    TABWeCrit_(0.0)
+    y0_(this->coeffDict().template lookupOrDefault<scalar>("y0", 0.0)),
+    yDot0_(this->coeffDict().template lookupOrDefault<scalar>("yDot0", 0.0)),
+    TABComega_(8),
+    TABCmu_(5),
+    TABtwoWeCrit_(12)
 {
-    if (solveOscillationEq_)
+    if (solveOscillationEq_ && dict.found("TABCoeffs"))
     {
         const dictionary coeffs(dict.subDict("TABCoeffs"));
-        y0_ = coeffs.template lookupOrDefault<scalar>("y0", 0.0);
-        yDot0_ = coeffs.template lookupOrDefault<scalar>("yDot0", 0.0);
-        TABComega_ = coeffs.template lookupOrDefault<scalar>("Comega", 8.0);
-        TABCmu_ = coeffs.template lookupOrDefault<scalar>("Cmu", 10.0);
-        TABWeCrit_ = coeffs.template lookupOrDefault<scalar>("WeCrit", 12.0);
+        coeffs.lookup("Comega") >> TABComega_;
+        coeffs.lookup("Cmu") >> TABCmu_;
+        scalar WeCrit(readScalar(coeffs.lookup("WeCrit")));
+        TABtwoWeCrit_ = 2*WeCrit;
     }
 }
 
@@ -160,4 +159,3 @@ bool Foam::BreakupModel<CloudType>::update
 #include "BreakupModelNew.C"
 
 // ************************************************************************* //
-
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
index ec51693b03bcab3b3f125536c846d5da38ae15e6..fdbd849b609b0fa3aa8aa58da7c77b31da06846b 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,9 +63,10 @@ protected:
 
         scalar y0_;
         scalar yDot0_;
+
         scalar TABComega_;
         scalar TABCmu_;
-        scalar TABWeCrit_;
+        scalar TABtwoWeCrit_;
 
 
 public:
@@ -153,9 +154,9 @@ public:
             return TABCmu_;
         }
 
-        inline const scalar& TABWeCrit() const
+        inline const scalar& TABtwoWeCrit() const
         {
-            return TABWeCrit_;
+            return TABtwoWeCrit_;
         }
 
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
index dec1e0bae0fb986b7f5248e250870239e8a6c4b9..c30ded6d7c89553fe6e28a33afbcf575437ca698 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,7 +116,7 @@ bool Foam::ETAB<CloudType>::update
         scalar romega = 1.0/omega;
 
         scalar We = rhoc*sqr(Urmag)*r/sigma;
-        scalar Wetmp = We/this->TABWeCrit_;
+        scalar Wetmp = We/this->TABtwoWeCrit_;
 
         scalar y1 = y - Wetmp;
         scalar y2 = yDot*romega;
diff --git a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
index 817542ad7cf23d354a82c76638b7f9b5c3fdf575..92a92157ad31707993d7eba0200bd5430e1ef53b 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,7 +128,7 @@ bool Foam::TAB<CloudType>::update
     {
         scalar omega = sqrt(omega2);
         scalar We = rhoc*sqr(Urmag)*r/sigma;
-        scalar Wetmp = We/this->TABWeCrit_;
+        scalar Wetmp = We/this->TABtwoWeCrit_;
 
         scalar y1 = y - Wetmp;
         scalar y2 = yDot/omega;