diff --git a/applications/solvers/modules/compressibleVoF/alphaPredictor.C b/applications/solvers/modules/compressibleVoF/alphaPredictor.C
index 13b6562868..335182aa60 100644
--- a/applications/solvers/modules/compressibleVoF/alphaPredictor.C
+++ b/applications/solvers/modules/compressibleVoF/alphaPredictor.C
@@ -36,8 +36,6 @@ void Foam::solvers::compressibleVoF::alphaPredictor()
 {
     #include "alphaControls.H"
 
-    volScalarField& alpha2(mixture.alpha2());
-
     const volScalarField& rho1 = mixture.thermo1().rho();
     const volScalarField& rho2 = mixture.thermo2().rho();
 
@@ -47,30 +45,6 @@ void Foam::solvers::compressibleVoF::alphaPredictor()
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();
 
-        talphaPhi1 = new surfaceScalarField
-        (
-            IOobject
-            (
-                "alphaPhi1",
-                runTime.name(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar(alphaPhi1.dimensions(), 0)
-        );
-
-        surfaceScalarField rhoPhiSum
-        (
-            IOobject
-            (
-                "rhoPhiSum",
-                runTime.name(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar(rhoPhi.dimensions(), 0)
-        );
-
         tmp<volScalarField> trSubDeltaT;
 
         if (LTS)
@@ -79,31 +53,59 @@ void Foam::solvers::compressibleVoF::alphaPredictor()
                 fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
         }
 
+        // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
+        tmp<surfaceScalarField> talphaPhi1
+        (
+            surfaceScalarField::New
+            (
+                "alphaPhi1",
+                mesh,
+                dimensionedScalar(alphaPhi1.dimensions(), 0)
+            )
+        );
+
+        // Sub-cycle on both alpha1 and alpha2
+        List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
+
         for
         (
-            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            subCycle<volScalarField, subCycleFields> alphaSubCycle
+            (
+                alphaPtrs,
+                nAlphaSubCycles
+            );
             !(++alphaSubCycle).end();
         )
         {
             #include "alphaEqn.H"
             talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
-            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
         }
 
         alphaPhi1 = talphaPhi1();
-        rhoPhi = rhoPhiSum;
     }
     else
     {
         #include "alphaEqn.H"
     }
 
-    contErr =
-        (
-            fvc::ddt(rho)()() + fvc::div(rhoPhi)()()
-          - (fvModels().source(alpha1, rho1)&rho1)()
-          - (fvModels().source(alpha2, rho2)&rho2)()
-        );
+    mixture.correct();
+
+    alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi1;
+    alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi1);
+
+    rhoPhi = alphaRhoPhi1 + alphaRhoPhi2;
+
+    contErr1 =
+    (
+        fvc::ddt(alpha1, rho1)()() + fvc::div(alphaRhoPhi1)()()
+      - (fvModels().source(alpha1, rho1)&rho1)()
+    );
+
+    contErr2 =
+    (
+        fvc::ddt(alpha2, rho2)()() + fvc::div(alphaRhoPhi2)()()
+      - (fvModels().source(alpha2, rho2)&rho2)()
+    );
 }
 
 
diff --git a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C
index c81bdefc3e..8c7eb7ba27 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C
+++ b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C
@@ -34,13 +34,17 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
     const surfaceScalarField& phi,
     const surfaceScalarField& rhoPhi,
     const surfaceScalarField& alphaPhi1,
+    const surfaceScalarField& alphaRhoPhi1,
+    const surfaceScalarField& alphaRhoPhi2,
     const compressibleTwoPhaseMixture& mixture
 )
 :
     twoPhaseTransport_(false),
     mixture_(mixture),
     phi_(phi),
-    alphaPhi1_(alphaPhi1)
+    alphaPhi1_(alphaPhi1),
+    alphaRhoPhi1_(alphaRhoPhi1),
+    alphaRhoPhi2_(alphaRhoPhi2)
 {
     {
         IOdictionary momentumTransport
@@ -68,38 +72,14 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
 
     if (twoPhaseTransport_)
     {
-        const volScalarField& alpha1(mixture_.alpha1());
-        const volScalarField& alpha2(mixture_.alpha2());
-
-        const volScalarField& rho1 = mixture_.thermo1().rho();
-        const volScalarField& rho2 = mixture_.thermo2().rho();
-
-        alphaRhoPhi1_ =
-        (
-            new surfaceScalarField
-            (
-                IOobject::groupName("alphaRhoPhi", alpha1.group()),
-                fvc::interpolate(rho1)*alphaPhi1_
-            )
-        );
-
-        alphaRhoPhi2_ =
-        (
-            new surfaceScalarField
-            (
-                IOobject::groupName("alphaRhoPhi", alpha2.group()),
-                fvc::interpolate(rho2)*(phi_ - alphaPhi1_)
-            )
-        );
-
         momentumTransport1_ =
         (
             phaseCompressible::momentumTransportModel::New
             (
-                alpha1,
-                rho1,
+                mixture_.alpha1(),
+                mixture_.thermo1().rho(),
                 U,
-                alphaRhoPhi1_(),
+                alphaRhoPhi1_,
                 phi,
                 mixture.thermo1()
             )
@@ -109,10 +89,10 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
         (
             phaseCompressible::momentumTransportModel::New
             (
-                alpha2,
-                rho2,
+                mixture_.alpha2(),
+                mixture_.thermo2().rho(),
                 U,
-                alphaRhoPhi2_(),
+                alphaRhoPhi2_,
                 phi,
                 mixture.thermo2()
             )
@@ -154,19 +134,6 @@ Foam::compressibleInterPhaseTransportModel::divDevTau
 }
 
 
-void Foam::compressibleInterPhaseTransportModel::correctPhasePhi()
-{
-    if (twoPhaseTransport_)
-    {
-        const volScalarField& rho1 = mixture_.thermo1().rho();
-        const volScalarField& rho2 = mixture_.thermo2().rho();
-
-        alphaRhoPhi1_.ref() = fvc::interpolate(rho1)*alphaPhi1_;
-        alphaRhoPhi2_.ref() = fvc::interpolate(rho2)*(phi_ - alphaPhi1_);
-    }
-}
-
-
 void Foam::compressibleInterPhaseTransportModel::predict()
 {
     if (twoPhaseTransport_)
diff --git a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H
index bc3d779386..7dc4188e0b 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H
+++ b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H
@@ -77,10 +77,10 @@ class compressibleInterPhaseTransportModel
         const surfaceScalarField& alphaPhi1_;
 
         //- Phase-1 mass-flux (constructed for two-phase transport)
-        tmp<surfaceScalarField> alphaRhoPhi1_;
+        const surfaceScalarField& alphaRhoPhi1_;
 
         //- Phase-2 mass-flux (constructed for two-phase transport)
-        tmp<surfaceScalarField> alphaRhoPhi2_;
+        const surfaceScalarField& alphaRhoPhi2_;
 
         //- Mixture transport model (constructed for mixture transport)
         autoPtr<compressible::momentumTransportModel> mixtureMomentumTransport_;
@@ -106,6 +106,8 @@ public:
             const surfaceScalarField& phi,
             const surfaceScalarField& rhoPhi,
             const surfaceScalarField& alphaPhi1,
+            const surfaceScalarField& alphaRhoPhi1,
+            const surfaceScalarField& alphaRhoPhi2,
             const compressibleTwoPhaseMixture& mixture
         );
 
@@ -121,10 +123,6 @@ public:
         //- Return the effective momentum stress divergence
         tmp<fvVectorMatrix> divDevTau(volVectorField& U) const;
 
-        //- Correct the phase mass-fluxes
-        // (required for the two-phase transport option)
-        void correctPhasePhi();
-
         //- Predict the phase or mixture transport models
         void predict();
 
diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.C b/applications/solvers/modules/compressibleVoF/compressibleVoF.C
index 951e952989..b75c780762 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleVoF.C
+++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.C
@@ -96,6 +96,7 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
     mixture(U),
 
     alpha1(mixture.alpha1()),
+    alpha2(mixture.alpha2()),
 
     alphaRestart
     (
@@ -169,6 +170,18 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         phi*fvc::interpolate(alpha1)
     ),
 
+    alphaRhoPhi1
+    (
+        IOobject::groupName("alphaRhoPhi", alpha1.group()),
+        fvc::interpolate(mixture.thermo1().rho())*alphaPhi1
+    ),
+
+    alphaRhoPhi2
+    (
+        IOobject::groupName("alphaRhoPhi", alpha2.group()),
+        fvc::interpolate(mixture.thermo2().rho())*(phi - alphaPhi1)
+    ),
+
     K("K", 0.5*magSqr(U)),
 
     momentumTransport
@@ -178,6 +191,8 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         phi,
         rhoPhi,
         alphaPhi1,
+        alphaRhoPhi1,
+        alphaRhoPhi2,
         mixture
     ),
 
@@ -330,7 +345,6 @@ void Foam::solvers::compressibleVoF::prePredictor()
 {
     fvModels().correct();
     alphaPredictor();
-    momentumTransport.correctPhasePhi();
 
     if (pimple.predictTransport())
     {
diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.H b/applications/solvers/modules/compressibleVoF/compressibleVoF.H
index 4f340c6e30..5eea866d6b 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleVoF.H
+++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.H
@@ -95,9 +95,12 @@ protected:
         //- The compressible two-phase mixture
         immiscibleCompressibleTwoPhaseMixture mixture;
 
-        //- Reference to the primary phase-fraction
+        //- Reference to the phase1-fraction
         volScalarField& alpha1;
 
+        //- Reference to the phase2-fraction
+        volScalarField& alpha2;
+
         //- Switch indicating if this is a restart
         bool alphaRestart;
 
@@ -137,6 +140,12 @@ protected:
         // Phase-1 volumetric flux
         surfaceScalarField alphaPhi1;
 
+        // Phase-1 mass-flux
+        surfaceScalarField alphaRhoPhi1;
+
+        // Phase-2 mass-flux
+        surfaceScalarField alphaRhoPhi2;
+
         //- Kinetic energy field
         //  Used in the energy equation
         volScalarField K;
@@ -164,7 +173,9 @@ protected:
 
         tmp<volScalarField> rAU;
 
-        tmp<volScalarField::Internal> contErr;
+        tmp<volScalarField::Internal> contErr1;
+
+        tmp<volScalarField::Internal> contErr2;
 
         //- MULES Correction
         tmp<surfaceScalarField> talphaPhi1Corr0;
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C
index 5f91f13f28..0f4a8991be 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C
@@ -89,7 +89,7 @@ Foam::fv::VoFClouds::VoFClouds
 
 Foam::wordList Foam::fv::VoFClouds::addSupFields() const
 {
-    return wordList({thermo_.rho()().name(), "U", "T"});
+    return wordList({thermo_.rho()().name(), thermo_.he().name(), "U"});
 }
 
 
@@ -108,7 +108,7 @@ void Foam::fv::VoFClouds::correct()
 
 void Foam::fv::VoFClouds::addSup
 (
-    const volScalarField& rho,
+    const volScalarField& alpha,
     fvMatrix<scalar>& eqn,
     const word& fieldName
 ) const
@@ -122,13 +122,31 @@ void Foam::fv::VoFClouds::addSup
     {
         eqn += clouds_.Srho();
     }
-    else if (fieldName == "T")
+    else
     {
-        const volScalarField::Internal Cv(thermo_.Cv());
+        FatalErrorInFunction
+            << "Support for field " << fieldName << " is not implemented"
+            << exit(FatalError);
+    }
+}
 
-        eqn +=
-            clouds_.Sh(eqn.psi())()/Cv
-          + clouds_.Srho()*(eqn.psi() - thermo_.he()/Cv);
+
+void Foam::fv::VoFClouds::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (fieldName == thermo_.he().name())
+    {
+        eqn += clouds_.Sh(eqn.psi());
     }
     else
     {
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
index ed5c571da1..9a9ae34924 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
@@ -126,15 +126,24 @@ public:
 
         // Add explicit and implicit contributions to compressible equation
 
-            //- Add explicit contribution to compressible enthalpy equation
+            //- Add explicit contribution to phase continuity
             virtual void addSup
             (
+                const volScalarField& alpha,
+                fvMatrix<scalar>& eqn,
+                const word& fieldName
+            ) const;
+
+            //- Add explicit contribution to phase energy equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const word& fieldName
             ) const;
 
-            //- Add implicit contribution to compressible momentum equation
+            //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
index 5d0bba9327..834f765fbd 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
@@ -92,6 +92,15 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
     relax_(NaN),
     Cu_(NaN),
     q_(NaN),
+
+    thermo_
+    (
+        mesh().lookupObject<compressibleTwoPhaseMixture>
+        (
+            "phaseProperties"
+        )
+    ),
+
     alphaSolid_
     (
         IOobject
@@ -115,12 +124,13 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
 
 Foam::wordList Foam::fv::VoFSolidificationMeltingSource::addSupFields() const
 {
-    return wordList({"U", "T"});
+    return wordList({"U", thermo_.thermo1().he().name()});
 }
 
 
 void Foam::fv::VoFSolidificationMeltingSource::addSup
 (
+    const volScalarField& alpha,
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const word& fieldName
@@ -131,24 +141,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    const compressibleTwoPhaseMixture& thermo
-    (
-        mesh().lookupObject<compressibleTwoPhaseMixture>
-        (
-            "phaseProperties"
-        )
-    );
-
-    const volScalarField CpVoF(thermo.thermo1().Cp());
-
-    if (eqn.psi().dimensions() == dimTemperature)
-    {
-        eqn += L_/CpVoF*(fvc::ddt(rho, alphaSolid_));
-    }
-    else
-    {
-        eqn += L_*(fvc::ddt(rho, alphaSolid_));
-    }
+    eqn += L_*(fvc::ddt(rho, alphaSolid_));
 }
 
 
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
index 3ad77e8564..deaa07b63f 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
@@ -88,6 +88,9 @@ SourceFiles
 
 namespace Foam
 {
+
+class compressibleTwoPhaseMixture;
+
 namespace fv
 {
 
@@ -119,6 +122,9 @@ class VoFSolidificationMeltingSource
         //- Coefficient used in porosity calc - default = 0.001
         scalar q_;
 
+        //- VoF thermo
+        const compressibleTwoPhaseMixture& thermo_;
+
         //- Solid phase fraction
         mutable volScalarField alphaSolid_;
 
@@ -171,15 +177,16 @@ public:
 
         // Add explicit and implicit contributions to compressible equation
 
-            //- Add explicit contribution to compressible enthalpy equation
+            //- Add explicit contribution to phase energy equation
             virtual void addSup
             (
+                const volScalarField& alpha,
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const word& fieldName
             ) const;
 
-            //- Add implicit contribution to compressible momentum equation
+            //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C
index 4c66ed7596..ccadc1f3e9 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C
@@ -78,7 +78,7 @@ Foam::wordList Foam::fv::VoFSurfaceFilm::addSupFields() const
         {
             surfaceFilm_.rhoPrimary().name(),
             surfaceFilm_.UPrimary().name(),
-            surfaceFilm_.TPrimary().name()
+            surfaceFilm_.primaryThermo().he().name()
         }
     );
 }
@@ -105,7 +105,7 @@ void Foam::fv::VoFSurfaceFilm::correct()
 
 void Foam::fv::VoFSurfaceFilm::addSup
 (
-    const volScalarField& rho,
+    const volScalarField& alpha,
     fvMatrix<scalar>& eqn,
     const word& fieldName
 ) const
@@ -119,14 +119,31 @@ void Foam::fv::VoFSurfaceFilm::addSup
     {
         eqn += surfaceFilm_.Srho();
     }
-    else if (fieldName == surfaceFilm_.TPrimary().name())
+    else
     {
-        const volScalarField::Internal Cv(surfaceFilm_.primaryThermo().Cv());
+        FatalErrorInFunction
+            << "Support for field " << fieldName << " is not implemented"
+            << exit(FatalError);
+    }
+}
 
-        eqn +=
-            surfaceFilm_.Sh()()/Cv
-          + surfaceFilm_.Srho()
-           *(eqn.psi() - surfaceFilm_.primaryThermo().he()/Cv);
+
+void Foam::fv::VoFSurfaceFilm::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (fieldName == surfaceFilm_.primaryThermo().he().name())
+    {
+        eqn += surfaceFilm_.Sh();
     }
     else
     {
diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H
index 8a6d0c9353..33eadabeb5 100644
--- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H
+++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H
@@ -116,15 +116,24 @@ public:
 
         // Add explicit and implicit contributions to compressible equation
 
-            //- Add explicit contribution to compressible enthalpy equation
+            //- Add explicit contribution to phase continuity
             virtual void addSup
             (
+                const volScalarField& alpha,
+                fvMatrix<scalar>& eqn,
+                const word& fieldName
+            ) const;
+
+            //- Add explicit contribution to phase energy equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const word& fieldName
             ) const;
 
-            //- Add implicit contribution to compressible momentum equation
+            //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
diff --git a/applications/solvers/modules/compressibleVoF/momentumPredictor.C b/applications/solvers/modules/compressibleVoF/momentumPredictor.C
index 273895bb81..e221c34d3b 100644
--- a/applications/solvers/modules/compressibleVoF/momentumPredictor.C
+++ b/applications/solvers/modules/compressibleVoF/momentumPredictor.C
@@ -31,7 +31,8 @@ void Foam::solvers::compressibleVoF::momentumPredictor()
 {
     tUEqn =
     (
-        fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::Sp(contErr(), U)
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      - fvm::Sp(contErr1() + contErr2(), U)
       + MRF.DDt(rho, U)
       + momentumTransport.divDevTau(U)
      ==
diff --git a/applications/solvers/modules/compressibleVoF/pressureCorrector.C b/applications/solvers/modules/compressibleVoF/pressureCorrector.C
index 14b3169a84..5295844ab1 100644
--- a/applications/solvers/modules/compressibleVoF/pressureCorrector.C
+++ b/applications/solvers/modules/compressibleVoF/pressureCorrector.C
@@ -34,8 +34,6 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
 {
     volScalarField& p = mixture.p();
 
-    const volScalarField& alpha2(mixture.alpha2());
-
     const volScalarField& rho1 = mixture.rho1();
     const volScalarField& rho2 = mixture.rho2();
 
diff --git a/applications/solvers/modules/compressibleVoF/rhofs.H b/applications/solvers/modules/compressibleVoF/rhofs.H
deleted file mode 100644
index 67dc10f378..0000000000
--- a/applications/solvers/modules/compressibleVoF/rhofs.H
+++ /dev/null
@@ -1,2 +0,0 @@
-surfaceScalarField rho1f(fvc::interpolate(rho1));
-surfaceScalarField rho2f(fvc::interpolate(rho2));
diff --git a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
index 108d7bfe74..169f1598a9 100644
--- a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
+++ b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
@@ -29,29 +29,61 @@ License
 
 void Foam::solvers::compressibleVoF::thermophysicalPredictor()
 {
+    const volScalarField& rho1(mixture.rho1());
+    const volScalarField& rho2(mixture.rho2());
+
+    const volScalarField& e1(mixture.thermo1().he());
+    const volScalarField& e2(mixture.thermo2().he());
+
+    const fvScalarMatrix e1Source(fvModels().source(alpha1, rho1, e1));
+    const fvScalarMatrix e2Source(fvModels().source(alpha2, rho2, e2));
+
     volScalarField& T = mixture.T();
     const volScalarField& p = mixture.p();
-    const volScalarField& alpha2(mixture.alpha2());
 
     fvScalarMatrix TEqn
     (
-        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr(), T)
-      - fvm::laplacian(thermophysicalTransport.alphaEff(), T)
+        correction
+        (
+            mixture.thermo1().Cv()()
+           *(
+                fvm::ddt(alpha1, rho1, T) + fvm::div(alphaRhoPhi1, T)
+              - (
+                    e1Source.hasDiag()
+                  ? fvm::Sp(contErr1(), T) + fvm::Sp(e1Source.A(), T)
+                  : fvm::Sp(contErr1(), T)
+                )
+            )
+          + mixture.thermo2().Cv()()
+           *(
+                fvm::ddt(alpha2, rho2, T) + fvm::div(alphaRhoPhi2, T)
+              - (
+                    e2Source.hasDiag()
+                  ? fvm::Sp(contErr2(), T) + fvm::Sp(e2Source.A(), T)
+                  : fvm::Sp(contErr2(), T)
+                )
+            )
+        )
+
+      + fvc::ddt(alpha1, rho1, e1) + fvc::div(alphaRhoPhi1, e1)
+      - contErr1()*e1
+      + fvc::ddt(alpha2, rho2, e2) + fvc::div(alphaRhoPhi2, e2)
+      - contErr2()*e2
+
+      - fvm::laplacian(thermophysicalTransport.kappaEff(), T)
+
       + (
-             mixture.totalInternalEnergy()
+            mixture.totalInternalEnergy()
           ?
-            fvc::div(fvc::absolute(phi, U), p)()() // - contErr()/rho*p
+            fvc::div(fvc::absolute(phi, U), p)()()
           + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
-          - (U()&(fvModels().source(rho, U)&U)()) - contErr()*K
+          - (U()&(fvModels().source(rho, U)&U)()) - (contErr1() + contErr2())*K
           :
             p*fvc::div(fvc::absolute(phi, U))()()
         )
-       *(
-           alpha1()/mixture.thermo1().Cv()()
-         + alpha2()/mixture.thermo2().Cv()()
-        )
      ==
-        fvModels().source(rho, T)
+        (e1Source&e1)
+      + (e2Source&e2)
     );
 
     TEqn.relax();
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
new file mode 100644
index 0000000000..a5477a0058
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -0,0 +1,50 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
+    // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
+    tmp<surfaceScalarField> talphaPhi1
+    (
+        surfaceScalarField::New
+        (
+            "alphaPhi1",
+            mesh,
+            dimensionedScalar(alphaPhi1.dimensions(), 0)
+        )
+    );
+
+    List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
+
+    for
+    (
+        subCycle<volScalarField, subCycleFields> alphaSubCycle
+        (
+            alphaPtrs,
+            nAlphaSubCycles
+        );
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
+    }
+
+    alphaPhi1 = talphaPhi1();
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+mixture.correct();
+rho == mixture.rho();
+
+// Calculate the mass-flux from the accumulated alphaPhi1
+rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);
diff --git a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C
index 96a2a22ae4..72ade6824b 100644
--- a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C
+++ b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C
@@ -32,14 +32,14 @@ incompressibleInterPhaseTransportModel
 (
     const volVectorField& U,
     const surfaceScalarField& phi,
-    const surfaceScalarField& alphaPhi10,
+    const surfaceScalarField& alphaPhi1,
     const incompressibleTwoPhaseMixture& mixture
 )
 :
     twoPhaseTransport_(false),
     mixture_(mixture),
     phi_(phi),
-    alphaPhi10_(alphaPhi10)
+    alphaPhi1_(alphaPhi1)
 {
     {
         IOdictionary momentumTransport
@@ -75,7 +75,7 @@ incompressibleInterPhaseTransportModel
             new surfaceScalarField
             (
                 IOobject::groupName("alphaPhi", alpha2.group()),
-                (phi_ - alphaPhi10_)
+                (phi_ - alphaPhi1_)
             )
         );
 
@@ -85,7 +85,7 @@ incompressibleInterPhaseTransportModel
             (
                 alpha1,
                 U,
-                alphaPhi10_,
+                alphaPhi1_,
                 phi,
                 mixture.nuModel1()
             )
@@ -143,7 +143,7 @@ void Foam::incompressibleInterPhaseTransportModel::correctPhasePhi()
 {
     if (twoPhaseTransport_)
     {
-        alphaPhi2_.ref() = (phi_ - alphaPhi10_);
+        alphaPhi2_.ref() = (phi_ - alphaPhi1_);
     }
 }
 
diff --git a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H
index e8a7cb09ef..69b9759e4b 100644
--- a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H
+++ b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H
@@ -70,9 +70,9 @@ class incompressibleInterPhaseTransportModel
         const surfaceScalarField& phi_;
 
         //- Phase volumetric flux
-        const surfaceScalarField& alphaPhi10_;
+        const surfaceScalarField& alphaPhi1_;
 
-        //- Phase-2 mass-flux (constructed from alphaPhi10_ and phi_)
+        //- Phase-2 mass-flux (constructed from alphaPhi1_ and phi_)
         tmp<surfaceScalarField> alphaPhi2_;
 
         //- Mixture transport model (constructed for mixture transport)
@@ -94,7 +94,7 @@ public:
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
-            const surfaceScalarField& alphaPhi10,
+            const surfaceScalarField& alphaPhi1,
             const incompressibleTwoPhaseMixture& mixture
         );
 
diff --git a/src/twoPhaseModels/VoF/alphaEqn.H b/src/twoPhaseModels/VoF/alphaEqn.H
index c4ac7d876c..2dbf2ff896 100644
--- a/src/twoPhaseModels/VoF/alphaEqn.H
+++ b/src/twoPhaseModels/VoF/alphaEqn.H
@@ -215,7 +215,8 @@
 
         alpha2 = 1.0 - alpha1;
 
-        mixture.correct();
+        // Correct only the mixture interface for the interface compression flux
+        mixture.interfaceProperties::correct();
     }
 
     if (alphaApplyPrevCorr && MULESCorr)
@@ -228,19 +229,13 @@
         talphaPhi1Corr0.clear();
     }
 
-    #include "rhofs.H"
-
     if
     (
         word(mesh.schemes().ddt("ddt(rho,U)"))
-     == fv::EulerDdtScheme<vector>::typeName
-     || word(mesh.schemes().ddt("ddt(rho,U)"))
-     == fv::localEulerDdtScheme<vector>::typeName
+     != fv::EulerDdtScheme<vector>::typeName
+     && word(mesh.schemes().ddt("ddt(rho,U)"))
+     != fv::localEulerDdtScheme<vector>::typeName
     )
-    {
-        rhoPhi = alphaPhi1*(rho1f - rho2f) + phiCN*rho2f;
-    }
-    else
     {
         if (ocCoeff > 0)
         {
@@ -248,9 +243,6 @@
             alphaPhi1 =
                 (alphaPhi1 - (1.0 - cnCoeff)*alphaPhi1.oldTime())/cnCoeff;
         }
-
-        // Calculate the end-of-time-step mass flux
-        rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "
diff --git a/src/twoPhaseModels/VoF/alphaEqnSubCycle.H b/src/twoPhaseModels/VoF/alphaEqnSubCycle.H
deleted file mode 100644
index bdd4070f87..0000000000
--- a/src/twoPhaseModels/VoF/alphaEqnSubCycle.H
+++ /dev/null
@@ -1,41 +0,0 @@
-if (nAlphaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    surfaceScalarField rhoPhiSum
-    (
-        IOobject
-        (
-            "rhoPhiSum",
-            runTime.name(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(rhoPhi.dimensions(), 0)
-    );
-
-    tmp<volScalarField> trSubDeltaT;
-
-    if (LTS)
-    {
-        trSubDeltaT =
-            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
-    }
-
-    for
-    (
-        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-        !(++alphaSubCycle).end();
-    )
-    {
-        #include "alphaEqn.H"
-        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-    }
-
-    rhoPhi = rhoPhiSum;
-}
-else
-{
-    #include "alphaEqn.H"
-}
-
-rho == mixture.rho();
diff --git a/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C b/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C
index d54d42cbe6..f549ab60bd 100644
--- a/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C
+++ b/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C
@@ -40,7 +40,6 @@ immiscibleCompressibleTwoPhaseMixture(const volVectorField& U)
 void Foam::immiscibleCompressibleTwoPhaseMixture::correct()
 {
     compressibleTwoPhaseMixture::correct();
-    interfaceProperties::correct();
 }
 
 
diff --git a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
index d64172c18a..ce74646565 100644
--- a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
+++ b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
@@ -40,7 +40,6 @@ immiscibleIncompressibleTwoPhaseMixture(const volVectorField& U)
 void Foam::immiscibleIncompressibleTwoPhaseMixture::correct()
 {
     incompressibleTwoPhaseMixture::correct();
-    interfaceProperties::correct();
 }
 
 
diff --git a/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes b/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes
index 8ecc19544b..83cef88b09 100644
--- a/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes
+++ b/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes
@@ -36,7 +36,8 @@ divSchemes
     div(rhoPhi,U)       Gauss linearUpwind limited;
 
     turbulence          Gauss limitedLinear 1;
-    div(rhoPhi,T)       $turbulence;
+    div(alphaRhoPhi,e)  $turbulence;
+    div(alphaRhoPhi,T)  Gauss upwind;
     div(rhoPhi,K)       $turbulence;
     div(phi,p)          $turbulence;
     div(phi,k)          $turbulence;
diff --git a/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes b/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes
index 610de57c32..1bad0a8d9d 100644
--- a/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes
@@ -33,10 +33,11 @@ divSchemes
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
     div(rhoPhi,U)   Gauss limitedLinear 1;
     div(rhoPhi,K)   Gauss limitedLinear 1;
-    div(rhoPhi,T)   Gauss limitedLinear 1;
+    div(phi,p)      Gauss upwind;
+    div(alphaRhoPhi,e) Gauss limitedLinear 1;
+    div(alphaRhoPhi,T) Gauss upwind;
     div(rhoPhi,k)   Gauss upwind;
     div(rhoPhi,epsilon) Gauss upwind;
-    div(phi,p)      Gauss upwind;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes b/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes
index 7ff41630b5..437b790d32 100644
--- a/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes
@@ -26,11 +26,13 @@ gradSchemes
 
 divSchemes
 {
-    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
-    div(rhoPhi,T)   Gauss linearUpwind grad(T);
-
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
+    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
+
+    div(alphaRhoPhi,e) Gauss linearUpwind grad(e);
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(phi,p)      Gauss upwind;
     div(rhoPhi,K)   Gauss upwind;
 
diff --git a/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes b/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes
index 0d7d1c0432..76856d7176 100644
--- a/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes
@@ -26,11 +26,13 @@ gradSchemes
 
 divSchemes
 {
-    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
-    div(rhoPhi,T)   Gauss linearUpwind grad(T);
-
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
+    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
+
+    div(alphaRhoPhi,e) Gauss linearUpwind grad(e);
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(phi,p)      Gauss upwind;
     div(rhoPhi,K)   Gauss upwind;
 
diff --git a/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes b/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes
index 9dbe1506f2..3b8ab02076 100644
--- a/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes
@@ -29,7 +29,10 @@ divSchemes
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
     div(rhoPhi,U)   Gauss linearUpwind grad(U);Gauss upwind;
-    div(rhoPhi,T)   Gauss upwind;
+
+    div(alphaRhoPhi,e) Gauss upwind;
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(rhoPhi,K)   Gauss upwind;
     div(phi,p)      Gauss upwind;
 
diff --git a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes
index 27525c071d..25cf33d862 100644
--- a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes
@@ -29,7 +29,10 @@ divSchemes
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
     div(rhoPhi,U)  Gauss upwind;
-    div(rhoPhi,T)  Gauss upwind;
+
+    div(alphaRhoPhi,e) Gauss upwind;
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(rhoPhi,K)  Gauss upwind;
     div(phi,p)      Gauss upwind;
     div(phi,k)      Gauss upwind;
diff --git a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution
index fa770e4225..99015d85b1 100644
--- a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution
+++ b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution
@@ -97,8 +97,8 @@ PIMPLE
 {
     momentumPredictor no;
     transonic       no;
-    nOuterCorrectors 2;
-    nCorrectors     2;
+    nOuterCorrectors 3;
+    nCorrectors     1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes
index 27525c071d..25cf33d862 100644
--- a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes
@@ -29,7 +29,10 @@ divSchemes
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
     div(rhoPhi,U)  Gauss upwind;
-    div(rhoPhi,T)  Gauss upwind;
+
+    div(alphaRhoPhi,e) Gauss upwind;
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(rhoPhi,K)  Gauss upwind;
     div(phi,p)      Gauss upwind;
     div(phi,k)      Gauss upwind;
diff --git a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution
index 1428d681d5..82b21c5e3c 100644
--- a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution
+++ b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution
@@ -98,7 +98,7 @@ PIMPLE
     momentumPredictor no;
     transonic       no;
     nOuterCorrectors 3;
-    nCorrectors     2;
+    nCorrectors     1;
     nNonOrthogonalCorrectors 0;
 }
 
diff --git a/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes b/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes
index 0d7d1c0432..76856d7176 100644
--- a/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes
@@ -26,11 +26,13 @@ gradSchemes
 
 divSchemes
 {
-    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
-    div(rhoPhi,T)   Gauss linearUpwind grad(T);
-
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
+    div(rhoPhi,U)   Gauss linearUpwindV grad(U);
+
+    div(alphaRhoPhi,e) Gauss linearUpwind grad(e);
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(phi,p)      Gauss upwind;
     div(rhoPhi,K)   Gauss upwind;
 
diff --git a/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes b/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes
index 2271bc4c13..50faa84e82 100644
--- a/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes
+++ b/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes
@@ -29,7 +29,10 @@ divSchemes
     div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
 
     div(rhoPhi,U)  Gauss vanLeerV;
-    div(rhoPhi,T)  Gauss vanLeer;
+
+    div(alphaRhoPhi,e) Gauss vanLeer;
+    div(alphaRhoPhi,T) Gauss upwind;
+
     div(rhoPhi,K)  Gauss linear;
     div(phi,p)     Gauss linear;