diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index bf79e250ef..b3046c7fbb 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -26,8 +26,9 @@ Application
     interDyMFoam
 
 Description
-    Solver for 2 incompressible fluids, which tracks the interface and includes
-    the option of mesh motion and mesh topology changes including adaptive
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach,
+    with optional mesh motion and mesh topology changes including adaptive
     re-meshing.
 
 \*---------------------------------------------------------------------------*/
@@ -45,18 +46,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createDynamicFvMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,10 +64,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readControls.H"
-#       include "CourantNo.H"
+        #include "readControls.H"
+        #include "CourantNo.H"
 
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -96,7 +96,7 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && correctPhi)
         {
-#           include "correctPhi.H"
+            #include "correctPhi.H"
         }
 
         // Make the fluxes relative to the mesh motion
@@ -107,22 +107,22 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && checkMeshCourantNo)
         {
-#           include "meshCourantNo.H"
+            #include "meshCourantNo.H"
         }
 
         twoPhaseProperties.correct();
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
 
         p = pd + rho*gh;
 
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 239779ce2e..887bcbdead 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -45,7 +45,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -83,14 +83,29 @@
     );
 
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    surfaceScalarField ghf("gh", g & mesh.Cf());
-
-
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 89b7ecc5cd..f693111680 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -45,18 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,10 +63,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -75,17 +74,19 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
+
+        p = pd + rho*gh;
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index de9baf4d7e..b57c02a8cb 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -48,18 +48,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,10 +66,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -78,10 +77,11 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties->correct();
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
         turbulence->correct();
 
+        // --- Outer-corrector loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
             #include "UEqn.H"
diff --git a/applications/solvers/multiphase/lesInterFoam/createFields.H b/applications/solvers/multiphase/lesInterFoam/createFields.H
index 666e655871..c24f8bda35 100644
--- a/applications/solvers/multiphase/lesInterFoam/createFields.H
+++ b/applications/solvers/multiphase/lesInterFoam/createFields.H
@@ -44,7 +44,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -87,9 +87,24 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
 
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
 
diff --git a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
index c006fc6dd6..acb11e1b29 100644
--- a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
+++ b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
@@ -26,8 +26,11 @@ Application
     lesInterFoam
 
 Description
-    Solver for 2 incompressible fluids capturing the interface.  Turbulence is
-    modelled using a runtime selectable incompressible LES model.
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.  Turbulence is modelled using a run-time
+    selectable incompressible LES model.
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,28 +41,21 @@ Description
 #include "twoPhaseMixture.H"
 #include "incompressible/LESmodel/LESmodel.H"
 
-#include "IFstream.H"
-#include "OFstream.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-
-#   include "createFields.H"
-//#   include "createAverages.H"
-
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,31 +63,29 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
         turbulence->correct();
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr < nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
-//#       include "calculateAverages.H"
+        #include "continuityErrs.H"
 
         runTime.write();
-//#       include "writeNaveragingSteps.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/multiphase/rasInterFoam/createFields.H b/applications/solvers/multiphase/rasInterFoam/createFields.H
index 1aef6db579..988e8b3f74 100644
--- a/applications/solvers/multiphase/rasInterFoam/createFields.H
+++ b/applications/solvers/multiphase/rasInterFoam/createFields.H
@@ -44,7 +44,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -88,9 +88,24 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
 
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
 
diff --git a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C b/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
index 9cf47c5f32..1c78f70bd8 100644
--- a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
+++ b/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
@@ -23,11 +23,14 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    lesInterFoam
+    rasInterFoam
 
 Description
-    Solver for 2 incompressible fluids capturing the interface.  Turbulence is
-    modelled using a runtime selectable incompressible RAS model.
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.  Turbulence is modelled using a run-time
+    selectable incompressible RAS model.
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,18 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,26 +63,28 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr < nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
+
+        p = pd + rho*gh;
 
         turbulence->correct();
 
diff --git a/etc/controlDict b/etc/controlDict
index cb5cee6ddf..a52db70f3c 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -434,7 +434,7 @@ DebugSwitches
     featureEdgeMesh     0;
     fieldToCell         0;
     file                0;
-    fileName            0;
+    fileName            2;
     filteredLinear      0;
     filteredLinear2     0;
     filteredLinear2V    0;
@@ -840,7 +840,7 @@ DebugSwitches
     waveTransmissive    0;
     wedge               0;
     weighted            0;
-    word                0;
+    word                2;
     writer              0;
     xmgr                0;
     zeroGradient        0;