diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 11a1809575..e24de6658d 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -74,24 +74,26 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
-
-        // Do any mesh changes
-        mesh.update();
-
-        if (mesh.changing())
         {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
+            // Store divU from the previous mesh for the correctPhi
+            volScalarField divU = fvc::div(phi);
 
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
+            scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        if (mesh.changing() && correctPhi)
-        {
-            //***HGW#include "correctPhi.H"
+            // Do any mesh changes
+            mesh.update();
+
+            if (mesh.changing())
+            {
+                Info<< "Execution time for mesh.update() = "
+                    << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                    << " s" << endl;
+            }
+
+            if (mesh.changing() && correctPhi)
+            {
+                #include "correctPhi.H"
+            }
         }
 
         // Make the fluxes relative to the mesh motion
@@ -102,6 +104,12 @@ int main(int argc, char *argv[])
             #include "meshCourantNo.H"
         }
 
+        if (mesh.changing())
+        {
+            gh = g & mesh.C();
+            ghf = g & mesh.Cf();
+        }
+
         turbulence->correct();
 
         // --- Outer-corrector loop
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
new file mode 100644
index 0000000000..4152105784
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
@@ -0,0 +1,39 @@
+{
+    #include "continuityErrs.H"
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+
+    adjustPhi(phi, U, pcorr);
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
+        );
+
+        pcorrEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
index 1f579d245b..dc04fb4546 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
@@ -150,3 +150,14 @@
     (
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
+
+
+    wordList pcorrTypes(pd.boundaryField().types());
+
+    for (label i=0; i<pd.boundaryField().size(); i++)
+    {
+        if (pd.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }