diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H
index c040545617..68744f7043 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H
@@ -1,7 +1,7 @@
 {
     radiation->correct();
     rhoCp =  rho*fluid.Cp();
-    
+
     const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi);
 
     const volScalarField kappaEff
@@ -31,7 +31,7 @@
     fvOptions.correct(T);
 
     fluid.correct();
-    
+
     Info<< "min/max(T) = "
         << min(T).value() << ", " << max(T).value() << endl;
 }
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C
index 190b7fbb1d..9f5a3b03ce 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C
@@ -59,14 +59,14 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Kexp
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
-        
+
         const volScalarField coeff
         (
             C_*from*this->pair().from().rho()*pos(from - alphaMin_)
            *(refValue - Tactivate_)
            /Tactivate_
         );
-        
+
         if (sign(C_.value()) > 0)
         {
             return
@@ -99,7 +99,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSp
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
-        
+
         const volScalarField coeff
         (
             C_*from*this->pair().from().rho()*pos(from - alphaMin_)
@@ -142,7 +142,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSu
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
-        
+
         const volScalarField coeff
         (
             C_*from*this->pair().from().rho()*pos(from - alphaMin_)
@@ -179,7 +179,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Tactivate() const
 
 
 template<class Thermo, class OtherThermo>
-bool 
+bool
 Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::includeDivU()
 {
     return true;
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
index 97059d4f29..e1f6013373 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
@@ -148,14 +148,14 @@ public:
         (
             const volScalarField& field
         );
-        
+
         //- Implicit mass transfer coefficient
         virtual tmp<volScalarField> KSp
         (
             label modelVariable,
             const volScalarField& field
         );
-        
+
         //- Explicit mass transfer coefficient
         virtual tmp<volScalarField> KSu
         (
@@ -165,9 +165,9 @@ public:
 
         //- Return T transition between phases
         virtual const dimensionedScalar& Tactivate() const;
-        
+
         //- Adds and substract alpha*div(U) as a source term
-        //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) 
+        //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
         virtual bool includeDivU();
 };
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
index 498290199b..2260b1db1d 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
@@ -176,25 +176,24 @@ public:
         (
             const volScalarField& field
         ) = 0;
-        
+
         //- Implicit mass transfer
         virtual tmp<volScalarField> KSp
         (
             label modelVariable,
             const volScalarField& field
         ) = 0;
-        
+
         //- Explicit mass transfer
         virtual tmp<volScalarField> KSu
         (
             label modelVariable,
             const volScalarField& field
         ) = 0;
-        
 
         //- Reference value
         virtual const dimensionedScalar& Tactivate() const = 0;
-        
+
         //- Adds and substract alpha*div(U) as a source term
         //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
         virtual bool includeDivU();
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
index 4c85bf2e03..417c51f226 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
@@ -187,14 +187,14 @@ public:
         (
             const volScalarField& field
         );
-        
+
         //- Implicit mass transfer coefficient
         virtual tmp<volScalarField> KSp
         (
             label modelVariable,
             const volScalarField& field
         );
-        
+
         //- Explicit mass transfer coefficient
         virtual tmp<volScalarField> KSu
         (
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
index 463f59da9e..8128eb068a 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
@@ -35,7 +35,7 @@
           - fvm::laplacian(rAUf, p_rgh)
           + fluid.volTransfer(p_rgh)
         );
-          
+
         p_rghEqn.setReference(pRefCell, pRefValue);
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
index 73a763bc3b..99777692bf 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
@@ -219,7 +219,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                     )
                 );
                 volScalarField& dmdtNetki = tdmdtNetki.ref();
-                
+
                 tmp<volScalarField> tSp
                 (
                     new volScalarField
@@ -235,7 +235,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                     )
                 );
                 volScalarField& Sp = tSp.ref();
-                
+
                 tmp<volScalarField> tSu
                 (
                     new volScalarField
@@ -259,23 +259,23 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                         massTransferModels_[keyik];
 
                     dmdtNetki -= *dmdt_[keyik];
-                    
+
                     tmp<volScalarField> KSp =
                         interfacePtr->KSp(interfaceCompositionModel::T, T);
-                    
+
                     if (KSp.valid())
                     {
                         Sp -= KSp.ref();
                     }
-                    
+
                     tmp<volScalarField> KSu =
                         interfacePtr->KSu(interfaceCompositionModel::T, T);
-                    
+
                     if (KSu.valid())
                     {
                         Su -= KSu.ref();
                     }
-                    
+
                     // If linearization is not provided used full explicit
                     if (!KSp.valid() && !KSu.valid())
                     {
@@ -291,23 +291,23 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
 
                     dmdtNetki += *dmdt_[keyki];
 
-                    
+
                     tmp<volScalarField> KSp =
                         interfacePtr->KSp(interfaceCompositionModel::T, T);
-                    
+
                     if (KSp.valid())
                     {
                         Sp += KSp.ref();
                     }
-                    
+
                     tmp<volScalarField> KSu =
                         interfacePtr->KSu(interfaceCompositionModel::T, T);
-                    
+
                     if (KSu.valid())
                     {
                         Su += KSu.ref();
                     }
-                    
+
                     // If linearization is not provided used full explicit
                     if (!KSp.valid() && !KSu.valid())
                     {
@@ -337,9 +337,9 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
     (
         new fvScalarMatrix(p, dimVolume/dimTime)
     );
-    
+
     fvScalarMatrix& eqn = tEqnPtr.ref();
-    
+
     tmp<volScalarField> tSp
     (
         new volScalarField
@@ -355,7 +355,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
         )
     );
     volScalarField& Sp = tSp.ref();
-    
+
     tmp<volScalarField> tSu
     (
         new volScalarField
@@ -371,7 +371,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
         )
     );
     volScalarField& Su = tSu.ref();
-    
+
     forAllConstIters(this->totalPhasePairs(), iter)
     {
         const phasePair& pair = iter()();
@@ -385,72 +385,72 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
             phase2.name(),
             true
         );
-        
+
         if (massTransferModels_.found(key12))
         {
             autoPtr<interfaceCompositionModel>& interfacePtr =
                 massTransferModels_[key12];
-           
+
             tmp<volScalarField> KSp =
                 interfacePtr->KSp(interfaceCompositionModel::P, p);
-            
+
             if (KSp.valid())
             {
                 Sp += KSp.ref();
             }
-            
+
             tmp<volScalarField> KSu =
                 interfacePtr->KSu(interfaceCompositionModel::P, p);
-            
+
             if (KSu.valid())
             {
                 Su += KSu.ref();
             }
-            
+
             // If linearization is not provided used full explicit
             if (!KSp.valid() && !KSu.valid())
             {
-                Su -= 
-                    *dmdt_[key12] 
+                Su -=
+                    *dmdt_[key12]
                     *(
                         - this->coeffs(phase1.name())
                         + this->coeffs(phase2.name())
                     );
             }
         }
-        
+
         const phasePairKey key21
         (
             phase2.name(),
             phase1.name(),
             true
         );
-        
+
         if (massTransferModels_.found(key21))
         {
             autoPtr<interfaceCompositionModel>& interfacePtr =
                 massTransferModels_[key21];
-           
+
             tmp<volScalarField> KSp =
                 interfacePtr->KSp(interfaceCompositionModel::P, p);
-            
+
             if (KSp.valid())
             {
                 Sp += KSp.ref();
             }
-            
+
             tmp<volScalarField> KSu =
                 interfacePtr->KSu(interfaceCompositionModel::P, p);
-            
+
             if (KSu.valid())
             {
                 Su += KSu.ref();
             }
-            
+
             // If linearization is not provided used full explicit
             if (!KSp.valid() && !KSu.valid())
             {
-                Su += 
+                Su +=
                     *dmdt_[key21]
                     *(
                         - this->coeffs(phase1.name())
@@ -460,7 +460,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
         }
 
     }
-    
+
     eqn += fvm::Sp(Sp, p) + Su;
     return tEqnPtr;
 }
@@ -489,18 +489,18 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::correctMassSources
 
                 // Phase k to phase i
                 const phasePairKey keyki(phasek.name(), phasei.name(), true);
-                
+
                 if (massTransferModels_.found(keyik))
                 {
                     autoPtr<interfaceCompositionModel>& interfacePtr =
                         massTransferModels_[keyik];
-                        
+
                     tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
-                    
+
                     *dmdt_[keyik] = Kexp.ref();
-                    
+
                 }
-                
+
                 if (massTransferModels_.found(keyki))
                 {
                     autoPtr<interfaceCompositionModel>& interfacePtr =
@@ -508,7 +508,7 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::correctMassSources
 
                     // Explicit temperature mass transfer rate
                     const tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
-                    
+
                     *dmdt_[keyki] = Kexp.ref();
                 }
             }
@@ -527,7 +527,7 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
     // This term adds and substract alpha*div(U) as a source term
     // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
     bool includeDivU(true);
-    
+
     forAllConstIters(this->totalPhasePairs(), iter)
     {
         const phasePair& pair = iter()();
@@ -551,26 +551,26 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
             phase2.name(),
             true
         );
-        
+
         tmp<volScalarField> tdmdt12(this->dmdt(key12));
         volScalarField& dmdt12 = tdmdt12.ref();
-        
+
         if (massTransferModels_.found(key12))
         {
             autoPtr<interfaceCompositionModel>& interfacePtr =
                 massTransferModels_[key12];
-             
+
             tmp<volScalarField> KSu =
                 interfacePtr->KSu(interfaceCompositionModel::alpha, phase1);
-            
+
             if (KSu.valid())
             {
                 dmdt12 = KSu.ref();
             }
-            
+
             includeDivU = interfacePtr->includeDivU();
         }
-       
+
 
         // Phase 2 to phase 1
         const phasePairKey key21
@@ -582,20 +582,20 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
 
         tmp<volScalarField> tdmdt21(this->dmdt(key21));
         volScalarField& dmdt21 = tdmdt21.ref();
-        
+
         if (massTransferModels_.found(key21))
         {
             autoPtr<interfaceCompositionModel>& interfacePtr =
                 massTransferModels_[key21];
-            
+
             tmp<volScalarField> KSu =
                 interfacePtr->KSu(interfaceCompositionModel::alpha, phase2);
-            
+
             if (KSu.valid())
             {
                 dmdt21 = KSu.ref();
             }
-            
+
             includeDivU = interfacePtr->includeDivU();
         }
 
@@ -612,16 +612,16 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
         const volScalarField coeffs12(coeffs1 - coeffs2);
 
         const surfaceScalarField& phi = this->phi();
-        
+
         if (includeDivU)
         {
             SuPhase1 +=
                 fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1));
-                
+
             SuPhase2 +=
                 fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1));
         }
-        
+
         // NOTE: dmdtNet is distributed in terms =
         //  Source for phase 1 =
         //      dmdtNet/rho1
@@ -715,12 +715,12 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
         }
 
         // Update ddtAlphaMax
-        this->ddtAlphaMax_ =  
+        this->ddtAlphaMax_ =
             max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
     }
 }
 
-        
+
 template<class BasePhaseSystem>
 void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer
 (
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
index 0a6e77ab10..030d02e992 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
@@ -67,8 +67,8 @@ public:
                 phasePairKey::hash
             >
             massTransferModelTable;
-            
-            
+
+
     typedef HashTable<volScalarField::Internal> SuSpTable;
 
 protected:
@@ -128,13 +128,13 @@ public:
 
         //- Return the heat transfer matrix
         virtual  tmp<fvScalarMatrix> heatTransfer(const volScalarField& T);
-        
+
         //- Return the volumetric rate transfer matrix
         virtual  tmp<fvScalarMatrix> volTransfer(const volScalarField& p);
-        
+
         //- Correct/calculates mass sources dmdt for phases
         virtual void correctMassSources(const volScalarField& T);
-        
+
         //- Calculate mass transfer for alpha's
         virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp);
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
index a6dcf73aa7..36066ac389 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
@@ -120,7 +120,7 @@ Foam::multiphaseSystem::multiphaseSystem
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::multiphaseSystem::calculateSuSp()
-{   
+{
     this->alphaTransfer(Su_, Sp_);
 }
 
@@ -129,9 +129,9 @@ void Foam::multiphaseSystem::solve()
 {
     const dictionary& alphaControls = mesh_.solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
-    
+
     volScalarField& alpha = phases_.first();
-    
+
     if (nAlphaSubCycles > 1)
     {
         surfaceScalarField rhoPhiSum
@@ -164,16 +164,16 @@ void Foam::multiphaseSystem::solve()
     {
         solveAlphas();
     }
-    
+
 }
 
 void Foam::multiphaseSystem::solveAlphas()
 {
-    
-    mesh_.solverDict("alpha").readEntry("cAlphas", cAlphas_);
+
     const dictionary& alphaControls = mesh_.solverDict("alpha");
+    alphaControls.readEntry("cAlphas", cAlphas_);
     label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
-    
+
     PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
 
     const surfaceScalarField& phi = this->phi();
@@ -361,7 +361,7 @@ void Foam::multiphaseSystem::solveAlphas()
             alpha1Eqn.solve();
 
             phiAlpha += alpha1Eqn.flux();
-            
+
             MULES::explicitSolve
             (
                 geometricOneField(),
@@ -375,7 +375,7 @@ void Foam::multiphaseSystem::solveAlphas()
             );
 
             phase.alphaPhi() = phiAlpha;
-            
+
             ++phasei;
         }
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.H
index 9b827bad61..e8329d7cf1 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.H
@@ -89,7 +89,7 @@ protected:
 
         //- Calculate Sp and Su
         void calculateSuSp();
-        
+
         //- Solve alphas
         void solveAlphas();
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem_backUp.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem_backUp.C
deleted file mode 100644
index 681de149fe..0000000000
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem_backUp.C
+++ /dev/null
@@ -1,642 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 "multiphaseSystem.H"
-
-#include "fixedValueFvsPatchFields.H"
-#include "Time.H"
-#include "subCycle.H"
-#include "fvcMeshPhi.H"
-
-#include "surfaceInterpolate.H"
-#include "fvcGrad.H"
-#include "fvcSnGrad.H"
-#include "fvcDiv.H"
-#include "fvcDdt.H"
-#include "fvcFlux.H"
-#include "fvmDdt.H"
-#include "fvcAverage.H"
-#include "fvMatrix.H"
-#include "fvmSup.H"
-#include "CMULES.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(multiphaseSystem, 0);
-    defineRunTimeSelectionTable(multiphaseSystem, dictionary);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::multiphaseSystem::multiphaseSystem
-(
-     const fvMesh& mesh
-)
-:
-    phaseSystem(mesh),
-    cAlphas_(),
-    ddtAlphaMax_(0.0),
-    limitedPhiAlphas_(phaseModels_.size()),
-    Su_(phaseModels_.size()),
-    Sp_(phaseModels_.size())
-{
-    label phasei = 0;
-    phases_.setSize(phaseModels_.size());
-    forAllIters(phaseModels_, iter)
-    {
-        phaseModel& pm = iter()();
-        phases_.set(phasei++, &pm);
-    }
-
-    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
-
-    // Initiate Su and Sp
-    forAllConstIters(phaseModels_, iter)
-    {
-        const phaseModel& pm = iter()();
-
-        Su_.insert
-        (
-            pm.name(),
-            volScalarField::Internal
-            (
-                IOobject
-                (
-                    "Su" + pm.name(),
-                    mesh_.time().timeName(),
-                    mesh_
-                ),
-                mesh_,
-                dimensionedScalar(dimless/dimTime, Zero)
-            )
-        );
-
-        Sp_.insert
-        (
-            pm.name(),
-            volScalarField::Internal
-            (
-                IOobject
-                (
-                    "Sp" + pm.name(),
-                    mesh_.time().timeName(),
-                    mesh_
-                ),
-                mesh_,
-                dimensionedScalar(dimless/dimTime, Zero)
-            )
-        );
-    }
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::multiphaseSystem::calculateSuSp()
-{
-    forAllConstIters(totalPhasePairs_, iter)
-    {
-        const phasePair& pair = iter()();
-
-        const phaseModel& phase1 = pair.phase1();
-        const phaseModel& phase2 = pair.phase2();
-
-        const volScalarField& alpha1 = pair.phase1();
-        const volScalarField& alpha2 = pair.phase2();
-
-        tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
-        const volScalarField&  coeffs1 = tCoeffs1();
-
-        tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
-        const volScalarField&  coeffs2 = tCoeffs2();
-
-        // Phase 1 to phase 2
-        const phasePairKey key12
-        (
-            phase1.name(),
-            phase2.name(),
-            true
-        );
-
-
-        tmp<volScalarField> tdmdt12(this->dmdt(key12));
-        const volScalarField& dmdt12 = tdmdt12();
-
-        // Phase 2 to phase 1
-        const phasePairKey key21
-        (
-            phase2.name(),
-            phase1.name(),
-            true
-        );
-
-        tmp<volScalarField> tdmdt21(this->dmdt(key21));
-        const volScalarField& dmdt21 = tdmdt21();
-
-        volScalarField::Internal& SpPhase1 = Sp_[phase1.name()];
-
-        volScalarField::Internal& SuPhase1 = Su_[phase1.name()];
-
-        volScalarField::Internal& SpPhase2 = Sp_[phase2.name()];
-
-        volScalarField::Internal& SuPhase2 = Su_[phase2.name()];
-
-        const volScalarField dmdtNet(dmdt21 - dmdt12);
-
-        const volScalarField coeffs12(coeffs1 - coeffs2);
-
-        // NOTE: dmdtNet is distributed in terms =
-        //  Source for phase 1 =
-        //      dmdtNet/rho1
-        //    - alpha1*dmdtNet(1/rho1 - 1/rho2)
-
-        forAll(dmdtNet, celli)
-        {
-            scalar dmdt21 = dmdtNet[celli];
-            scalar coeffs12Cell = coeffs12[celli];
-
-            scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
-
-            // exp.
-            SuPhase1[celli] += coeffs1[celli]*dmdt21;
-
-            if (dmdt21 > 0)
-            {
-                if (coeffs12Cell > 0)
-                {
-                    // imp
-                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
-                }
-                else if (coeffs12Cell < 0)
-                {
-                    // exp
-                    SuPhase1[celli] -=
-                        dmdt21*coeffs12Cell*alpha1Limited;
-                }
-            }
-            else if (dmdt21 < 0)
-            {
-                if (coeffs12Cell > 0)
-                {
-                    // exp
-                    SuPhase1[celli] -=
-                        dmdt21*coeffs12Cell*alpha1Limited;
-                }
-                else if (coeffs12Cell < 0)
-                {
-                    // imp
-                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
-                }
-            }
-        }
-
-        forAll(dmdtNet, celli)
-        {
-            scalar dmdt12 = -dmdtNet[celli];
-            scalar coeffs21Cell = -coeffs12[celli];
-
-            scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
-
-            // exp
-            SuPhase2[celli] += coeffs2[celli]*dmdt12;
-
-            if (dmdt12 > 0)
-            {
-                if (coeffs21Cell > 0)
-                {
-                    // imp
-                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
-                }
-                else if (coeffs21Cell < 0)
-                {
-                    // exp
-                    SuPhase2[celli] -=
-                        dmdt12*coeffs21Cell*alpha2Limited;
-                }
-            }
-            else if (dmdt12 < 0)
-            {
-                if (coeffs21Cell > 0)
-                {
-                    // exp
-                    SuPhase2[celli] -=
-                        coeffs21Cell*dmdt12*alpha2Limited;
-                }
-                else if (coeffs21Cell < 0)
-                {
-                    // imp
-                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
-                }
-            }
-
-        }
-
-        // Update ddtAlphaMax
-        ddtAlphaMax_ =  
-            max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
-    }
-}
-
-
-void Foam::multiphaseSystem::solve()
-{
-    const fvMesh& mesh = this->mesh();
-
-    const dictionary& alphaControls = mesh.solverDict("alpha");
-    label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
-    label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
-    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
-
-    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
-
-    const surfaceScalarField& phi = this->phi();
-
-    surfaceScalarField phic(mag((phi)/mesh_.magSf()));
-
-    // Do not compress interface at non-coupled boundary faces
-    // (inlets, outlets etc.)
-    surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
-    forAll(phic.boundaryField(), patchi)
-    {
-        fvsPatchScalarField& phicp = phicBf[patchi];
-
-        if (!phicp.coupled())
-        {
-            phicp == 0;
-        }
-    }
-
-    for (int acorr=0; acorr<nAlphaCorr; acorr++)
-    {
-        label phasei = 0;
-        for (phaseModel& phase1 : phases_)
-        {
-            const volScalarField& alpha1 = phase1;
-
-            phiAlphaCorrs.set
-            (
-                phasei,
-                new surfaceScalarField
-                (
-                    "phi" + alpha1.name() + "Corr",
-                    fvc::flux
-                    (
-                        phi,
-                        alpha1,
-                        "div(phi," + alpha1.name() + ')'
-                    )
-                )
-            );
-
-            surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
-
-            for (phaseModel& phase2 : phases_)
-            {
-                const volScalarField& alpha2 = phase2;
-
-                if (&phase2 == &phase1) continue;
-
-                const phasePairKey key12(phase1.name(), phase2.name());
-
-                if (!cAlphas_.found(key12))
-                {
-                    FatalErrorInFunction
-                        << "Phase compression factor (cAlpha) not found for : "
-                        << key12
-                        << exit(FatalError);
-                }
-                scalar cAlpha = cAlphas_.find(key12)();
-
-                phic = min(cAlpha*phic, max(phic));
-
-                surfaceScalarField phir(phic*nHatf(alpha1, alpha2));
-
-                word phirScheme
-                (
-                    "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
-                );
-
-                phiAlphaCorr += fvc::flux
-                (
-                   -fvc::flux(-phir, alpha2, phirScheme),
-                    alpha1,
-                    phirScheme
-                );
-            }
-
-            // Ensure that the flux at inflow BCs is preserved
-            forAll(phiAlphaCorr.boundaryField(), patchi)
-            {
-                fvsPatchScalarField& phiAlphaCorrp =
-                    phiAlphaCorr.boundaryFieldRef()[patchi];
-
-                if (!phiAlphaCorrp.coupled())
-                {
-                    const scalarField& phi1p = phi.boundaryField()[patchi];
-                    const scalarField& alpha1p =
-                        alpha1.boundaryField()[patchi];
-
-                    forAll(phiAlphaCorrp, facei)
-                    {
-                        if (phi1p[facei] < 0)
-                        {
-                            phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
-                        }
-                    }
-                }
-            }
-
-            ++phasei;
-        }
-
-        // Set Su and Sp to zero
-        for (const phaseModel& phase : phases_)
-        {
-            Su_[phase.name()] = dimensionedScalar("Su", dimless/dimTime, Zero);
-            Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
-
-            // Add alpha*div(U)
-            const volScalarField& alpha = phase;
-            Sp_[phase.name()] +=
-                fvc::div(phi);//*min(max(alpha, scalar(0)), scalar(1));
-        }
-
-
-        // Fill Su and Sp
-        calculateSuSp();
-
-        // Limit phiAlphaCorr on each phase
-        phasei = 0;
-        for (phaseModel& phase : phases_)
-        {
-            volScalarField& alpha1 = phase;
-
-            surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
-
-            volScalarField::Internal& Su = Su_[phase.name()];
-            volScalarField::Internal& Sp = Sp_[phase.name()];
-
-            MULES::limit
-            (
-                1.0/mesh_.time().deltaT().value(),
-                geometricOneField(),
-                alpha1,
-                phi,
-                phiAlphaCorr,
-                Sp,
-                Su,
-                oneField(),
-                zeroField(),
-                true
-            );
-            ++phasei;
-        }
-
-        MULES::limitSum(phiAlphaCorrs);
-
-        volScalarField sumAlpha
-        (
-            IOobject
-            (
-                "sumAlpha",
-                mesh_.time().timeName(),
-                mesh_
-            ),
-            mesh_,
-            dimensionedScalar(dimless, Zero)
-        );
-
-        phasei = 0;
-        for (phaseModel& phase : phases_)
-        {
-            volScalarField& alpha1 = phase;
-
-            const volScalarField::Internal& Su = Su_[phase.name()];
-
-            const volScalarField::Internal& Sp = Sp_[phase.name()];
-
-            surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
-
-            // Add a bounded upwind U-mean flux
-            phiAlpha += upwind<scalar>(mesh_, phi).flux(alpha1);
-//             fvScalarMatrix alpha1Eqn
-//             (
-//                 fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
-//               + fv::gaussConvectionScheme<scalar>
-//                 (
-//                     mesh,
-//                     phi,
-//                     upwind<scalar>(mesh, phi)
-//                 ).fvmDiv(phi, alpha1)
-//               ==
-//                  Su + fvm::Sp(Sp, alpha1)
-//             );
-// 
-//             alpha1Eqn.solve();
-
-            //phiAlpha += alpha1Eqn.flux();
-            
-            if (nAlphaSubCycles > 1)
-            {
-                for
-                (
-                    subCycle<volScalarField> alphaSubCycle
-                    (
-                        alpha1,
-                        nAlphaSubCycles
-                    );
-                    !(++alphaSubCycle).end();
-                )
-                {
-                    MULES::explicitSolve
-                    (
-                        geometricOneField(),
-                        alpha1,
-                        phi,
-                        phiAlpha,
-                        (alphaSubCycle.index()*Sp)(),
-                        (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
-                        oneField(),
-                        zeroField()
-                    );
-
-                    if (alphaSubCycle.index() == 1)
-                    {
-                        phase.alphaPhi() = phiAlpha;
-                    }
-                    else
-                    {
-                        phase.alphaPhi() += phiAlpha;
-                    }
-                }
-
-                phase.alphaPhi() /= nAlphaSubCycles;
-            }
-            else
-            {
-                MULES::explicitSolve
-                (
-                    geometricOneField(),
-                    alpha1,
-                    phi,
-                    phiAlpha,
-                    Sp,
-                    Su,
-                    oneField(),
-                    zeroField()
-                );
-
-                phase.alphaPhi() = phiAlpha;
-            }
-
-            ++phasei;
-        }
-
-        if (acorr == nAlphaCorr - 1)
-        {
-            volScalarField sumAlpha
-            (
-                IOobject
-                (
-                    "sumAlpha",
-                    mesh_.time().timeName(),
-                    mesh_
-                ),
-                mesh_,
-                dimensionedScalar(dimless, Zero)
-            );
-
-            // Reset rhoPhi
-            rhoPhi_ = dimensionedScalar("rhoPhi", dimMass/dimTime, Zero);
-
-            for (phaseModel& phase : phases_)
-            {
-                volScalarField& alpha1 = phase;
-                sumAlpha += alpha1;
-
-                // Update rhoPhi
-                rhoPhi_ += fvc::interpolate(phase.rho()) * phase.alphaPhi();
-            }
-
-            Info<< "Phase-sum volume fraction, min, max = "
-                << sumAlpha.weightedAverage(mesh_.V()).value()
-                << ' ' << min(sumAlpha).value()
-                << ' ' << max(sumAlpha).value()
-                << endl;
-
-            volScalarField sumCorr(1.0 - sumAlpha);
-
-            for (phaseModel& phase : phases_)
-            {
-                volScalarField& alpha = phase;
-                //alpha += alpha*sumCorr;
-
-                Info<< alpha.name() << " volume fraction = "
-                    << alpha.weightedAverage(mesh.V()).value()
-                    << "  Min(alpha) = " << min(alpha).value()
-                    << "  Max(alpha) = " << max(alpha).value()
-                    << endl;
-            }
-        }
-    }
-}
-
-
-const Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases() const
-{
-    return phases_;
-}
-
-
-Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases()
-{
-    return phases_;
-}
-
-
-const Foam::phaseModel& Foam::multiphaseSystem::phase(const label i) const
-{
-    return phases_[i];
-}
-
-
-Foam::phaseModel& Foam::multiphaseSystem::phase(const label i)
-{
-    return phases_[i];
-}
-
-
-Foam::dimensionedScalar Foam::multiphaseSystem::ddtAlphaMax() const
-{
-    return ddtAlphaMax_;
-}
-
-
-Foam::scalar Foam::multiphaseSystem::maxDiffNo() const
-{
-    auto iter = phaseModels_.cbegin();
-
-    scalar maxVal = max(iter()->diffNo()).value();
-
-    for (++iter; iter != phaseModels_.cend(); ++iter)
-    {
-        maxVal = max(maxVal, max(iter()->diffNo()).value());
-    }
-
-    return maxVal * mesh_.time().deltaT().value();
-}
-
-
-const Foam::multiphaseSystem::compressionFluxTable&
-Foam::multiphaseSystem::limitedPhiAlphas() const
-{
-    return limitedPhiAlphas_;
-}
-
-
-Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Su()
-{
-    return Su_;
-}
-
-
-Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Sp()
-{
-    return Sp_;
-}
-
-
-bool Foam::multiphaseSystem::read()
-{
-    return true;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
index 242d9cadec..395a4aa648 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
@@ -281,7 +281,7 @@ Foam::phaseSystem::phaseSystem
 
     // Total phase pair
     generatePairsTable();
-    
+
     // Update mu_
     calcMu();
 }
@@ -903,7 +903,7 @@ void Foam::phaseSystem::correct()
     {
         iter()->correct();
     }
-    
+
     calcMu();
 }
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
index d3a17be0ea..a5394939ba 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
@@ -78,7 +78,7 @@ public:
 
 
         typedef HashTable<autoPtr<phaseModel>> phaseModelTable;
-        
+
 
         typedef HashTable<volScalarField::Internal> SuSpTable;
 
@@ -116,7 +116,7 @@ protected:
 
         //- Reference to the mesh
         const fvMesh& mesh_;
-        
+
         //- Dynamic viscocity
         volScalarField mu_;
 
@@ -510,13 +510,13 @@ public:
         (
             const volScalarField& T
         ) = 0;
-        
+
         //- Return the volumetric rate transfer matrix
         virtual tmp<fvScalarMatrix> volTransfer
         (
             const volScalarField& p
         ) = 0;
-        
+
         //- Calculate mass transfer for alpha's
         virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp) = 0;
 
@@ -537,7 +537,7 @@ public:
 
         //- Correct the mixture thermos
         virtual void correct();
-        
+
         //- Correct mass sources
         virtual void correctMassSources(const volScalarField& T) = 0;
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/TEqn.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/TEqn.H
index 66eeb0e86c..23cab6e0fc 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/TEqn.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/TEqn.H
@@ -1,14 +1,14 @@
 {
     tmp<volScalarField> tcp(thermo->Cp());
     const volScalarField& cp = tcp();
-    
+
     const dimensionedScalar Cp1 = thermo->Cp1();
     const dimensionedScalar Cp2 = thermo->Cp2();
 
     rhoCp = rho*cp;
 
     kappaEff = thermo->kappa() + rho*cp*turbulence->nut()/Prt;
-    
+
     const surfaceScalarField rhoCpPhi
     (
         "rhoCpPhi",
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/interCondensatingEvaporatingFoam.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/interCondensatingEvaporatingFoam.C
index 9fad2b220b..2c3d939c7b 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/interCondensatingEvaporatingFoam.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/interCondensatingEvaporatingFoam.C
@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
     #include "createAlphaFluxes.H"
     #include "initCorrectPhi.H"
     #include "createUfIfPresent.H"
-    
+
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -99,13 +99,13 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readDyMControls.H"
-        
+
         // Store divU from the previous mesh so that it can be mapped
         // and used in correctPhi to ensure the corrected phi has the
         // same divergence
 
         volScalarField divU("divU", fvc::div(fvc::absolute(phi, U)));
-        
+
         {
             #include "CourantNo.H"
             #include "alphaCourantNo.H"
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
                     }
                 }
             }
-            
+
             mixture->correct();
 
             #include "alphaControls.H"
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
 
             #include "UEqn.H"
             #include "TEqn.H"
-            
+
             // --- Pressure corrector loop
             while (pimple.correct())
             {
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C
index 6350420fe5..5be98a1344 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C
@@ -58,14 +58,14 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::constant
     temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
     coeffC_
     (
-        "coeffC", 
-        dimless/dimTime/dimTemperature, 
+        "coeffC",
+        dimless/dimTime/dimTemperature,
         optionalSubDict(type() + "Coeffs")
     ),
     coeffE_
     (
-        "coeffE", 
-        dimless/dimTime/dimTemperature, 
+        "coeffE",
+        dimless/dimTime/dimTemperature,
         optionalSubDict(type() + "Coeffs")
     )
 {}
@@ -121,7 +121,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
     const dimensionedScalar& TSat = thermo.TSat();
 
     const dimensionedScalar T0(dimTemperature, Zero);
-    
+
     if (mesh_.time().outputTime())
     {
         volScalarField mDot
@@ -194,9 +194,9 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
         );
 
     const dimensionedScalar& TSat = thermo.TSat();
-    
+
     dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
-    
+
     volScalarField limitedAlpha1
     (
         min(max(mixture_.alpha1(), scalar(0)), scalar(1))
@@ -206,7 +206,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
     (
         min(max(mixture_.alpha2(), scalar(0)), scalar(1))
     );
-   
+
     const volScalarField Vcoeff
     (
         coeffE_*mixture_.rho1()*limitedAlpha1*L
@@ -216,7 +216,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
         coeffC_*mixture_.rho2()*limitedAlpha2*L
     );
 
-    TSource = 
+    TSource =
         fvm::Sp(Vcoeff, T) - Vcoeff*TSat
       - fvm::Sp(Ccoeff, T) + Ccoeff*TSat;
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.H
index 87d226c1f1..6c77da33e7 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.H
@@ -98,7 +98,7 @@ public:
         //  coefficient to multiply (Tsat - T) for the condensation rate
         //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
         virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
-        
+
         //- Source for T equarion
         virtual tmp<fvScalarMatrix> TSource() const;
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
index 468cf05c99..e786bfbb05 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
@@ -64,10 +64,10 @@ interfaceHeatResistance
     temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
     R_
     (
-        "R", 
+        "R",
         dimPower/dimArea/dimTemperature, optionalSubDict(type() + "Coeffs")
     ),
-    
+
     interfaceArea_
     (
         IOobject
@@ -81,7 +81,7 @@ interfaceHeatResistance
         mesh_,
         dimensionedScalar(dimless/dimLength, Zero)
     ),
-    
+
     mDotc_
     (
         IOobject
@@ -95,7 +95,7 @@ interfaceHeatResistance
         mesh_,
         dimensionedScalar(dimDensity/dimTime, Zero)
     ),
-    
+
     mDote_
     (
         IOobject
@@ -109,7 +109,7 @@ interfaceHeatResistance
         mesh_,
         dimensionedScalar(dimDensity/dimTime, Zero)
     ),
-    
+
     spread_
     (
         optionalSubDict(type() + "Coeffs").get<scalar>("spread")
@@ -148,7 +148,7 @@ mDotAlphal() const
     (
         min(max(mixture_.alpha2(), scalar(0)), scalar(1))
     );
-    
+
     return Pair<tmp<volScalarField>>
     (
         (mDotc_/(limitedAlpha2 + SMALL)),
@@ -161,10 +161,10 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 mDot() const
 {
-    
+
     return Pair<tmp<volScalarField>>
     (
-        tmp<volScalarField>(mDotc_), 
+        tmp<volScalarField>(mDotc_),
         tmp<volScalarField>(mDote_)
     );
 }
@@ -179,20 +179,20 @@ mDotDeltaT() const
         (
             mesh_.lookupObject<basicThermo>(basicThermo::dictName)
         );
-        
+
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
     const dimensionedScalar& TSat = thermo.TSat();
 
     Pair<tmp<volScalarField>> mDotce(mDot());
-    
+
     return Pair<tmp<volScalarField>>
     (
         mDotc_*pos(TSat - T.oldTime())/(TSat - T.oldTime()),
        -mDote_*pos(T.oldTime() - TSat)/(T.oldTime() - TSat)
     );
 
-    
+
 }
 
 
@@ -221,7 +221,7 @@ TSource() const
         );
 
     const dimensionedScalar& TSat = thermo.TSat();
-   
+
     // interface heat resistance
     volScalarField IHRcoeff = interfaceArea_*R_;
 
@@ -234,10 +234,10 @@ TSource() const
 void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 correct()
 {
-   
+
     // Update Interface
     updateInterface();
-    
+
     // Update mDotc_ and mDote_
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
@@ -251,16 +251,16 @@ correct()
     const dimensionedScalar T0(dimTemperature, Zero);
 
     dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
-   
+
     // interface heat resistance
     mDotc_ = interfaceArea_*R_*max(TSat - T, T0)/L;
     mDote_ = interfaceArea_*R_*(T - TSat)/L;
-       
+
     forAll(mDotc_, celli)
     {
         scalar rhobyDt = mixture_.rho1().value()/mesh_.time().deltaTValue();
         scalar maxEvap = mixture_.alpha1()[celli]*rhobyDt; // positive
-        scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative 
+        scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative
         mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap);
         mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap);
     }
@@ -278,7 +278,7 @@ updateInterface()
     );
 
     const dimensionedScalar& TSat = thermo.TSat();
-   
+
     // interface heat resistance
     // Interpolating alpha1 cell centre values to mesh points (vertices)
     scalarField ap
@@ -294,7 +294,7 @@ updateInterface()
         interfaceArea_[celli] = 0;
         if (status == 0) // cell is cut
         {
-            interfaceArea_[celli] = 
+            interfaceArea_[celli] =
                 mag(cutCell.isoFaceArea())/mesh_.V()[celli];
         }
     }
@@ -309,14 +309,14 @@ updateInterface()
             forAll(pp.faceCells(),i)
             {
                 const label pCelli = pp.faceCells()[i];
- 
+
                 if
                 (
-                    (TSat.value() - T[pCelli]) > 0 
+                    (TSat.value() - T[pCelli]) > 0
                   && mixture_.alpha1()[pCelli] < 0.9
                 )
                 {
-                    interfaceArea_[pCelli] = 
+                    interfaceArea_[pCelli] =
                         mag(pp.faceAreas()[i])/mesh_.V()[pCelli];
                 }
             }
@@ -331,23 +331,23 @@ vDot() const
 
     dimensionedScalar D
     (
-        "D", 
-        dimArea, 
+        "D",
+        dimArea,
         spread_/sqr(gAverage(mesh_.nonOrthDeltaCoeffs()))
     );
-    
-   
+
+
     const volScalarField& alpha1 = mixture_.alpha1();
     const volScalarField& alpha2 = mixture_.alpha2();
-    
+
     const dimensionedScalar MDotMin("MdotMin", mDotc_.dimensions(), 1e-3);
-    
+
     Pair<tmp<volScalarField>> mDotSpread
     (
-        tmp<volScalarField>(mDotc_*0.0), 
+        tmp<volScalarField>(mDotc_*0.0),
         tmp<volScalarField>(mDote_*0.0)
     );
-    
+
     if (max(mDotc_) > MDotMin)
     {
         fvc::spreadSource
@@ -373,18 +373,18 @@ vDot() const
             1e-3
         );
     }
-    
+
     dimensionedScalar pCoeff(1.0/mixture_.rho1() - 1.0/mixture_.rho2());
-    
+
     if (mesh_.time().outputTime())
     {
         volScalarField mDotS("mDotSpread", mDotSpread[1].ref());
         mDotS.write();
     }
-    
+
     return Pair<tmp<volScalarField>>
     (
-        pCoeff*mDotSpread[0], 
+        pCoeff*mDotSpread[0],
        -pCoeff*mDotSpread[1]
     );
 }
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
index 836f5a6f72..98b9e33b49 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
@@ -28,9 +28,9 @@ Class
     Foam::phaseChangeTwoPhaseMixtures::interfaceHeatResistance
 
 Description
-    Interface Heat Resistance type of condensation/saturation model using 
+    Interface Heat Resistance type of condensation/saturation model using
     spread source distribution following:
-    
+
     References:
     \verbatim
         Hardt, S., Wondra, F. (2008).
@@ -73,22 +73,22 @@ class interfaceHeatResistance
 
         //- Interface area
         volScalarField interfaceArea_;
-        
+
         //- Mass condensation source
         volScalarField mDotc_;
-        
+
         //- Mass evaporation source
         volScalarField mDote_;
-              
+
         //- Spread for mass source
         scalar spread_;
 
-        
+
     // Private member functions
 
         //- Update interface area
         void updateInterface();
-        
+
 
 public:
 
@@ -124,17 +124,17 @@ public:
         //  coefficient to multiply (Tsat - T) for the condensation rate
         //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
         virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
-        
+
         //- Source for T equarion
         virtual tmp<fvScalarMatrix> TSource() const;
-        
+
         //- Volumetric source for alpha (used by alphaEq)
         virtual Pair<tmp<volScalarField>> vDotAlphal() const;
-        
+
         //- Return the volumetric condensation and vaporisation rates as
         //  coefficients (used by p_rghEq)
         virtual Pair<tmp<volScalarField>> vDot() const;
-        
+
         //- Correct the interfaceHeatResistance phaseChange model
         virtual void correct();
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C
index 01c057ccf5..275cb4a4f1 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C
@@ -68,7 +68,7 @@ temperaturePhaseChangeTwoPhaseMixture
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixture::vDotAlphal() const
-{ 
+{
     volScalarField alphalCoeff
     (
         1.0/mixture_.rho1() - mixture_.alpha1()
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.H
index 8f14c82c96..14371c69cb 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.H
@@ -141,7 +141,7 @@ public:
         //  coefficient to multiply (Tsat - T) for the condensation rate
         //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
         virtual Pair<tmp<volScalarField>> mDotDeltaT() const = 0;
-        
+
         //- Source for T equarion
         virtual tmp<fvScalarMatrix> TSource() const = 0;
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 76949b8977..58feb3d700 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -114,7 +114,7 @@
     }
 
     rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
-    
+
     // Cache alphaPhi
     alphaPhi10 = talphaPhi();
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
index dad4ad165d..4a1d476c22 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
@@ -332,7 +332,7 @@ void Foam::fvc::spreadSource
 )
 {
     const fvMesh& mesh = alpha1.mesh();
-   
+
     volScalarField mDotSmear
     (
         IOobject
@@ -351,43 +351,43 @@ void Foam::fvc::spreadSource
     //- Smearing of source term field
     fvScalarMatrix mSourceEqn
     (
-        fvm::Sp(scalar(1), mDotSmear) 
-      - fvm::laplacian(D, mDotSmear) 
-      == 
+        fvm::Sp(scalar(1), mDotSmear)
+      - fvm::laplacian(D, mDotSmear)
+      ==
         mDotIn
     );
 
     mSourceEqn.solve();
 
-    // Cut cells with cutoff < alpha1 < 1-cutoff and rescale remaining 
+    // Cut cells with cutoff < alpha1 < 1-cutoff and rescale remaining
     // source term field
-    
+
     dimensionedScalar intvDotLiquid("intvDotLiquid", dimMass/dimTime, 0.0);
     dimensionedScalar intvDotVapor ("intvDotVapor", dimMass/dimTime, 0.0);
-    
+
     const scalarField& Vol = mesh.V();
-    
+
     forAll(mesh.C(), celli)
     {
         if (alpha1[celli] < cutoff)
         {
-            intvDotVapor.value() += 
+            intvDotVapor.value() +=
                 alpha2[celli]*mDotSmear[celli]*Vol[celli];
         }
         else if (alpha1[celli] > 1.0 - cutoff)
         {
-            intvDotLiquid.value() += 
+            intvDotLiquid.value() +=
                 alpha1[celli]*mDotSmear[celli]*Vol[celli];
         }
     }
-    
+
     reduce(intvDotVapor.value(), sumOp<scalar>());
     reduce(intvDotLiquid.value(), sumOp<scalar>());
 
     //- Calculate Nl and Nv
     dimensionedScalar Nl ("Nl", dimless, Zero);
-    dimensionedScalar Nv ("Nv", dimless, Zero); 
-    
+    dimensionedScalar Nv ("Nv", dimless, Zero);
+
     const dimensionedScalar intmSource0(fvc::domainIntegrate(mDotIn));
 
     if (intvDotVapor.value() > VSMALL)
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
index b9181a603b..6c2cb05fee 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
@@ -44,11 +44,11 @@ Description
            gradient of alpha is large (where the difference between the values
            in neighbouring cells is larger than alphaDiff) away from that
            starting point of the sweep.
-           
-           
+
+
     spreadSource: spread a source field (mDotIn) for two phase multiphase using
                   a laplacian operator and diffussivity D.
-                  The spread source (mDotOut) is distributed from alpha1 < cutoff 
+                  The spread source (mDotOut) is distributed from alpha1 < cutoff
                   to alpha1 > 1 - cutoff, and it is zero across the interface
 
 SourceFiles
@@ -91,7 +91,7 @@ namespace fvc
         const label nLayers,
         const scalar alphaDiff = 0.2
     );
-    
+
     void spreadSource
     (
         volScalarField& mDotOut,
diff --git a/src/functionObjects/field/interfaceHeight/interfaceHeight.C b/src/functionObjects/field/interfaceHeight/interfaceHeight.C
index 889f89c968..cbb331a146 100644
--- a/src/functionObjects/field/interfaceHeight/interfaceHeight.C
+++ b/src/functionObjects/field/interfaceHeight/interfaceHeight.C
@@ -55,7 +55,7 @@ void Foam::functionObjects::interfaceHeight::writePositions()
     const uniformDimensionedVectorField& g =
         mesh_.time().lookupObject<uniformDimensionedVectorField>("g");
     vector gHat = vector::zero;
-    
+
     if (mag(direction_) > 0.0)
     {
         gHat = direction_/mag(direction_);
@@ -237,7 +237,7 @@ Foam::functionObjects::interfaceHeight::interfaceHeight
 {
     read(dict);
     resetNames({"height", "position"});
-    
+
     writeFileHeader(fileID::heightFile);
     writeFileHeader(fileID::positionFile);
 }
@@ -253,7 +253,7 @@ Foam::functionObjects::interfaceHeight::~interfaceHeight()
 
 bool Foam::functionObjects::interfaceHeight::read(const dictionary& dict)
 {
-    
+
     dict.readIfPresent("alpha", alphaName_);
     dict.readIfPresent("liquid", liquid_);
     dict.lookup("locations") >> locations_;
diff --git a/src/functionObjects/field/interfaceHeight/interfaceHeight.H b/src/functionObjects/field/interfaceHeight/interfaceHeight.H
index b2d403baa6..7d882cca78 100644
--- a/src/functionObjects/field/interfaceHeight/interfaceHeight.H
+++ b/src/functionObjects/field/interfaceHeight/interfaceHeight.H
@@ -96,9 +96,8 @@ class interfaceHeight
 
         //- Interpolation scheme
         word interpolationScheme_;
-        
-        
-        //- Direction of interface motion 
+
+        //- Direction of interface motion
         vector direction_;
 
 
diff --git a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.gas b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.gas
index 4e623feb42..ccc3545d73 100755
--- a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.gas
+++ b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.gas
@@ -24,7 +24,7 @@ boundaryField
     {
         type            zeroGradient;
     }
-    
+
     outlet
     {
         type            inletOutlet;
diff --git a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.liquid b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.liquid
index 9908be1fb5..0cbcd69ae5 100755
--- a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.liquid
+++ b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.liquid
@@ -24,7 +24,7 @@ boundaryField
     {
         type            zeroGradient;
     }
-    
+
     outlet
     {
         type            inletOutlet;
diff --git a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allrun b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allrun
index ecd1f14402..ac16a2c2e9 100755
--- a/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allrun
+++ b/tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allrun
@@ -1,11 +1,12 @@
 #!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
+cd ${0%/*} || exit 1                                # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions            # Tutorial run functions
+#------------------------------------------------------------------------------
 
 runApplication blockMesh
 restore0Dir
+
+# copy 0 folder to 1.36
 cp -r 0 1.36
 cp system/setAlphaFieldDict.liquid system/setAlphaFieldDict
 runApplication setAlphaField
diff --git a/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allrun b/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allrun
index bf83717464..f508a34e2b 100755
--- a/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allrun
+++ b/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allrun
@@ -1,8 +1,7 @@
 #!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
+cd ${0%/*} || exit 1                                # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions            # Tutorial run functions
+#------------------------------------------------------------------------------
 
 runApplication blockMesh
 restore0Dir
diff --git a/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSchemes b/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSchemes
index fba9626cc6..aa26beac5f 100644
--- a/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSchemes
+++ b/tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSchemes
@@ -31,7 +31,7 @@ divSchemes
     div(rhoPhi,U)      Gauss vanLeerV;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
     div(rhoCpPhi,T) Gauss vanLeer;
-    
+
     div((interpolate(cp)*rhoPhi),T) Gauss linear;
 }