diff --git a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
index 2f7b843d2549d59bd5a64b4106ff28acdc03e9bb..dd4d382befaed467343166d527cb864dc9b6d72e 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
@@ -1,6 +1,5 @@
 {
-    volVectorField HbyA("HbyA", Uc);
-    HbyA = rAUc*UcEqn.H();
+    volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
 
     surfaceScalarField phiHbyA
     (
@@ -8,7 +7,6 @@
         (
            fvc::flux(HbyA)
          + alphacf*rAUcf*fvc::ddtCorr(Uc, Ucf)
-         + phicForces
         )
     );
 
@@ -19,6 +17,8 @@
         fvc::makeAbsolute(phiHbyA, Uc);
     }
 
+    phiHbyA += phicForces;
+
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, Uc, phiHbyA, rAUcf);
 
diff --git a/applications/solvers/lagrangian/DPMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/pEqn.H
index 9e465511221baf25438c93e5b2828c3aba915fd0..ea0bd1c101d9e27a89eea1db2f37bfb66e27d92f 100644
--- a/applications/solvers/lagrangian/DPMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/pEqn.H
@@ -1,6 +1,5 @@
 {
-    volVectorField HbyA("HbyA", Uc);
-    HbyA = rAUc*UcEqn.H();
+    volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
 
     surfaceScalarField phiHbyA
     (
@@ -8,10 +7,16 @@
         (
            fvc::flux(HbyA)
          + alphacf*rAUcf*fvc::ddtCorr(Uc, phic)
-         + phicForces
         )
     );
 
+    if (p.needReference())
+    {
+        adjustPhi(phiHbyA, Uc, p);
+    }
+
+    phiHbyA += phicForces;
+
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, Uc, phiHbyA, rAUcf);