diff --git a/applications/solvers/coupled/pUCoupledFoam/boundPU.H b/applications/solvers/coupled/pUCoupledFoam/boundPU.H
new file mode 100644
index 000000000..eca655e8c
--- /dev/null
+++ b/applications/solvers/coupled/pUCoupledFoam/boundPU.H
@@ -0,0 +1,32 @@
+{
+    // Bound the pressure
+    dimensionedScalar p1 = min(p);
+    dimensionedScalar p2 = max(p);
+
+    if (p1 < pMin || p2 > pMax)
+    {
+        Info<< "p: " << p1.value() << " " << p2.value()
+            << ".  Bounding." << endl;
+
+        p.max(pMin);
+        p.min(pMax);
+        p.correctBoundaryConditions();
+    }
+
+    // Bound the velocity
+    volScalarField magU = mag(U);
+    dimensionedScalar U1 = max(magU);
+
+    if (U1 > UMax)
+    {
+        Info<< "U: " << U1.value() << ".  Bounding." << endl;
+
+        volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU)
+            + neg(magU - UMax);
+        Ulimiter.max(scalar(0));
+        Ulimiter.min(scalar(1));
+
+        U *= Ulimiter;
+        U.correctBoundaryConditions();
+    }
+}
diff --git a/applications/solvers/coupled/pUCoupledFoam/convergenceCheck.H b/applications/solvers/coupled/pUCoupledFoam/convergenceCheck.H
new file mode 100644
index 000000000..9c7734131
--- /dev/null
+++ b/applications/solvers/coupled/pUCoupledFoam/convergenceCheck.H
@@ -0,0 +1,8 @@
+// Check convergence
+if (maxResidual < convergenceCriterion)
+{
+    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
+    runTime.writeAndEnd();
+    Info<< "latestTime = " << runTime.timeName() << endl;
+}
+
diff --git a/applications/solvers/coupled/pUCoupledFoam/initConvergenceCheck.H b/applications/solvers/coupled/pUCoupledFoam/initConvergenceCheck.H
new file mode 100644
index 000000000..482bc421d
--- /dev/null
+++ b/applications/solvers/coupled/pUCoupledFoam/initConvergenceCheck.H
@@ -0,0 +1,4 @@
+// initialize values for convergence checks
+
+    scalar maxResidual = 0;
+    scalar convergenceCriterion = 0;
diff --git a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C
index 088d41759..355b11ae4 100644
--- a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C
+++ b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C
@@ -52,11 +52,14 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readBlockSolverControls.H"
+#   include "initConvergenceCheck.H"
 
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
+#       include "readBlockSolverControls.H"
+#       include "readFieldBounds.H"
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         p.storePrevIter();
@@ -74,7 +77,7 @@ int main(int argc, char *argv[])
 #       include "couplingTerms.H"
 
         // Solve the block matrix
-        UpEqn.solve();
+        maxResidual = cmptMax(UpEqn.solve().initialResidual());
 
         // Retrieve solution
         UpEqn.retrieveSolution(0, U.internalField());
@@ -87,6 +90,8 @@ int main(int argc, char *argv[])
 
 #       include "continuityErrs.H"
 
+#       include "boundPU.H"
+
         p.relax();
 
         turbulence->correct();
@@ -95,6 +100,8 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
+
+#       include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H b/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H
index 4ee98e519..93c2bda47 100644
--- a/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H
+++ b/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H
@@ -1,3 +1,15 @@
-label pRefCell = 0;
-scalar pRefValue = 0;
-setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue);
+    label pRefCell = 0;
+    scalar pRefValue = 0;
+    setRefCell
+    (
+        p,
+        mesh.solutionDict().subDict("blockSolver"),
+        pRefCell,
+        pRefValue
+    );
+
+    mesh.solutionDict().subDict("blockSolver").readIfPresent
+    (
+        "convergence",
+        convergenceCriterion
+    );
diff --git a/applications/solvers/coupled/pUCoupledFoam/readFieldBounds.H b/applications/solvers/coupled/pUCoupledFoam/readFieldBounds.H
new file mode 100644
index 000000000..4a6c6f787
--- /dev/null
+++ b/applications/solvers/coupled/pUCoupledFoam/readFieldBounds.H
@@ -0,0 +1,14 @@
+    // Read field bounds
+    dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
+
+    // Pressure bounds
+    dimensionedScalar pMin("pMin", p.dimensions(), 0);
+    dimensionedScalar pMax("pMax", p.dimensions(), 0);
+
+    fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
+
+    // Velocity bound
+    dimensionedScalar UMax("UMax", U.dimensions(), 0);
+
+    fieldBounds.lookup(U.name()) >> UMax.value();
+    dimensionedScalar smallU("smallU", dimVelocity, 1e-10);