diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
index fae07312a494f8b85ddf62f0151c6124fa21aa35..a818ee2e9ea7707a7707b25f40026a2bb08ad591 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
@@ -68,11 +68,12 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
 
         if (g0.value() > 0.0)
         {
-            ppMagf = rAU1f*fvc::interpolate
-            (
-                (1.0/(rho1*(alpha1 + scalar(0.0001))))
-               *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-            );
+            ppMagf =
+                fvc::interpolate((1.0/rho1)*rAU1)
+               *fvc::interpolate
+                (
+                    g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+                );
 
             alpha1Eqn -= fvm::laplacian
             (
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
index c773650be287ff98c19d6754a11702b6066742ec..c969d827d5ffd781a6657125c2700021272651c3 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
@@ -333,11 +333,11 @@
         drag1
     );
 
-    surfaceScalarField rAU1f
+    volScalarField rAU1
     (
         IOobject
         (
-            "rAU1f",
+            "rAU1",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/files b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/files
index d8b7d88442a158fc7d19d058cda764a3e1525e39..2b36d0bc319b5a3b3541a93c316a671d508a4285 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/files
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/files
@@ -16,7 +16,6 @@ conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C
 radialModel/radialModel/radialModel.C
 radialModel/radialModel/newRadialModel.C
 radialModel/CarnahanStarling/CarnahanStarlingRadial.C
-radialModel/Gidaspow/GidaspowRadial.C
 radialModel/LunSavage/LunSavageRadial.C
 radialModel/SinclairJackson/SinclairJacksonRadial.C
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C
index 3c147d7d618a33c9b385e1dae6785fc66abe912e..e250a89e79b3a3c3446cf4fdfe9dc1c579a4d738 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,30 +69,29 @@ Foam::kineticTheoryModels::radialModels::CarnahanStarling::~CarnahanStarling()
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
 
     return
-        1.0/(1.0 - alpha1)
-      + 3.0*alpha1/(2.0*sqr(1.0 - alpha1))
-      + sqr(alpha1)/(2.0*pow(1.0 - alpha1, 3));
+        1.0/(1.0 - alpha)
+      + 3.0*alpha/(2.0*sqr(1.0 - alpha))
+      + sqr(alpha)/(2.0*pow(1.0 - alpha, 3));
 }
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0prime
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
     return
-        - alpha1/sqr(1.0 - alpha1)
-        + (3.0*(1.0 - alpha1) + 6.0*sqr(alpha1))/(2.0*(1.0 - alpha1))
-        + (2.0*alpha1*(1.0 - alpha1) + 3.0*pow(alpha1, 3))
-         /(2.0*pow(1.0 - alpha1, 4));
+        2.5/sqr(1.0 - alpha)
+      + 4.0*alpha/pow(1.0 - alpha, 3.0)
+      + 1.5*sqr(alpha)/pow(1.0 - alpha, 4.0);
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.H
index b4b825fe12483ca6b90f433130e66a8280379645..de5c24a1820785794e337500bfcecdcdd4c481b6 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,13 +75,13 @@ public:
 
         tmp<volScalarField> g0
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 };
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.C
deleted file mode 100644
index e1c642b4f9123017533fd53928ba1bd5a5edeed8..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.C
+++ /dev/null
@@ -1,93 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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 "GidaspowRadial.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace kineticTheoryModels
-{
-namespace radialModels
-{
-    defineTypeNameAndDebug(Gidaspow, 0);
-
-    addToRunTimeSelectionTable
-    (
-        radialModel,
-        Gidaspow,
-        dictionary
-    );
-}
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::kineticTheoryModels::radialModels::Gidaspow::Gidaspow
-(
-    const dictionary& dict
-)
-:
-    radialModel(dict)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::kineticTheoryModels::radialModels::Gidaspow::~Gidaspow()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::radialModels::Gidaspow::g0
-(
-    const volScalarField& alpha1,
-    const dimensionedScalar& alphaMax
-) const
-{
-    return 0.6/(1.0 - pow(alpha1/alphaMax, 1.0/3.0));
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::radialModels::Gidaspow::g0prime
-(
-    const volScalarField& alpha1,
-    const dimensionedScalar& alphaMax
-) const
-{
-    return
-        (-1.0/5.0)*pow(alpha1/alphaMax, -2.0/3.0)
-       /(alphaMax*sqr(1.0 - pow(alpha1/alphaMax, 1.0/3.0)));
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.H
deleted file mode 100644
index 60791274b7f3a3d36efe91a746ccd17694efc580..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/Gidaspow/GidaspowRadial.H
+++ /dev/null
@@ -1,99 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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::kineticTheoryModels::radialModels::Gidaspow
-
-Description
-
-SourceFiles
-    GidaspowRadial.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Gidaspow_H
-#define Gidaspow_H
-
-#include "radialModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace kineticTheoryModels
-{
-namespace radialModels
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class Gidaspow Declaration
-\*---------------------------------------------------------------------------*/
-
-class Gidaspow
-:
-    public radialModel
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("Gidaspow");
-
-
-    // Constructors
-
-        //- Construct from components
-        Gidaspow(const dictionary& dict);
-
-
-    //- Destructor
-    virtual ~Gidaspow();
-
-
-    // Member Functions
-
-        tmp<volScalarField> g0
-        (
-            const volScalarField& alpha1,
-            const dimensionedScalar& alphaMax
-        ) const;
-
-        tmp<volScalarField> g0prime
-        (
-            const volScalarField& alpha1,
-            const dimensionedScalar& alphaMax
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace radialModels
-} // End namespace kineticTheoryModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C
index e116a4b35e9cb4482e2dc5a2756f8578e9bba5e8..58ba2a6210e2244251da11a2a00915b339f956a3 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,23 +69,23 @@ Foam::kineticTheoryModels::radialModels::LunSavage::~LunSavage()
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::LunSavage::g0
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
 
-    return pow(1.0 - alpha1/alphaMax, -2.5*alphaMax);
+    return pow(1.0 - alpha/alphaMax, -2.5*alphaMax);
 }
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::LunSavage::g0prime
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
-    return 2.5*alphaMax*alpha1*pow(1.0 - alpha1, -1.0 - 2.5*alphaMax);
+    return 2.5*pow(1.0 - alpha/alphaMax, -1.0 - 2.5*alphaMax);
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.H
index 4977bfce5ed74de3e41a6064ec55809e875e53ab..1e95f838aecdf21c7da2aea675cb51af500de482 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,13 +74,13 @@ public:
 
         tmp<volScalarField> g0
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 };
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C
index cfbdd4baa99ee29a9b8cd5d8def4f1e70bc24753..74e8dac729ef04ab239d5cc73dc5f91e36c42d72 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,24 +69,24 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::~SinclairJackson()
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::SinclairJackson::g0
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
-    return 1.0/(1.0 - pow(alpha1/alphaMax, 1.0/3.0));
+    return 1.0/(1.0 - pow(alpha/alphaMax, 1.0/3.0));
 }
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
 (
-    const volScalarField& alpha1,
+    const volScalarField& alpha,
     const dimensionedScalar& alphaMax
 ) const
 {
     return
-        (-1.0/3.0)*pow(alpha1/alphaMax, -2.0/3.0)
-       /(alphaMax*sqr(1.0 - pow(alpha1/alphaMax, 1.0/3.0)));
+       (1.0/3.0)*pow(max(alpha, 1.0e-6)/alphaMax, -2.0/3.0)
+      /(alphaMax*sqr(1.0 - pow(alpha/alphaMax, 1.0/3.0)));
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.H
index 956cd568ef1a51859db46f6778ed3941dd684259..8df30f194513cc923060f1e5f5896b20f90fb094 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,13 +74,13 @@ public:
 
         tmp<volScalarField> g0
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const;
 };
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/radialModel/radialModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/radialModel/radialModel.H
index a53b8f1a376eaabceb370ae1f6d81d09e4585e53..bad05a99dff880b4adbb5b80c77af6ba03661763 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/radialModel/radialModel.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/radialModel/radialModel/radialModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,14 +107,14 @@ public:
         //- Radial distribution function
         virtual tmp<volScalarField> g0
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const = 0;
 
         //- Derivative of the radial distribution function
         virtual tmp<volScalarField> g0prime
         (
-            const volScalarField& alpha1,
+            const volScalarField& alpha,
             const dimensionedScalar& alphaMax
         ) const = 0;
 };
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index bf8fd458aa6e96df1a556e9c903ad1c121e26650..67717b07383af9ea70236a02140a8f7e62411556 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -10,7 +10,7 @@
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
-    volScalarField rAU1(1.0/U1Eqn.A());
+    rAU1 = 1.0/U1Eqn.A();
     volScalarField rAU2(1.0/U2Eqn.A());
 
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
@@ -30,6 +30,19 @@
       + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
       + rAlphaAU1f*(g & mesh.Sf())
     );
+
+    if (g0.value() > 0.0)
+    {
+        phiHbyA1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
+    }
+
+    if (kineticTheory.on())
+    {
+        phiHbyA1 -=
+            fvc::interpolate((1.0/rho1)*rAU1)
+           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
+    }
+
     mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2