diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
index e86b9e9bb280bbf913d2e4daf228f5c57020ba59..84ecd5d124f5cd16436d4caeae4a80be28484a14 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
@@ -35,12 +35,13 @@ template
     template<class> class SlicedPatchField,
     class GeoMesh
 >
-Foam::tmp<Foam::FieldField<PatchField, Type> > 
+Foam::tmp<Foam::FieldField<PatchField, Type> >
 Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
 slicedBoundaryField
 (
     const Mesh& mesh,
-    const Field<Type>& completeField
+    const Field<Type>& completeField,
+    const bool preserveCouples
 )
 {
     tmp<FieldField<PatchField, Type> > tbf
@@ -52,7 +53,7 @@ slicedBoundaryField
 
     forAll (mesh.boundary(), patchi)
     {
-        if (mesh.boundary()[patchi].coupled())
+        if (preserveCouples && mesh.boundary()[patchi].coupled())
         {
             // For coupled patched construct the correct patch field type
             bf.set
@@ -143,7 +144,8 @@ SlicedGeometricField
     const IOobject& io,
     const Mesh& mesh,
     const dimensionSet& ds,
-    const Field<Type>& completeField
+    const Field<Type>& completeField,
+    const bool preserveCouples
 )
 :
     GeometricField<Type, PatchField, GeoMesh>
@@ -152,7 +154,7 @@ SlicedGeometricField
         mesh,
         ds,
         Field<Type>(),
-        slicedBoundaryField(mesh, completeField)
+        slicedBoundaryField(mesh, completeField, preserveCouples)
     )
 {
     // Set the internalField to the slice of the complete field
@@ -179,7 +181,8 @@ SlicedGeometricField
     const Mesh& mesh,
     const dimensionSet& ds,
     const Field<Type>& completeIField,
-    const Field<Type>& completeBField
+    const Field<Type>& completeBField,
+    const bool preserveCouples
 )
 :
     GeometricField<Type, PatchField, GeoMesh>
@@ -188,7 +191,7 @@ SlicedGeometricField
         mesh,
         ds,
         Field<Type>(),
-        slicedBoundaryField(mesh, completeBField)
+        slicedBoundaryField(mesh, completeBField, preserveCouples)
     )
 {
     // Set the internalField to the slice of the complete field
diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
index 37eea9a3b58a37ecbb2ff11218b2f23d0db5539e..18f9abe05ccb4590721432cdd1b571d1d53ed31f 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
@@ -28,11 +28,11 @@ Class
 Description
     Specialization of GeometricField which holds slices of given complete
     fields in a form that they act as a GeometricField.
- 
+
     The destructor is wrapped to avoid deallocation of the storage of the
     complete fields when this is destroyed.
 
-    SlicedGeometricField can only be instantiated with a valid form of 
+    SlicedGeometricField can only be instantiated with a valid form of
     SlicedPatchField to handle the slicing and storage deallocation of the
     boundary field.
 
@@ -83,7 +83,8 @@ private:
         tmp<FieldField<PatchField, Type> >  slicedBoundaryField
         (
             const Mesh& mesh,
-            const Field<Type>& completeField
+            const Field<Type>& completeField,
+            const bool preserveCouples
         );
 
         //- Disallow default bitwise copy construct
@@ -111,7 +112,8 @@ public:
             const IOobject&,
             const Mesh&,
             const dimensionSet&,
-            const Field<Type>& completeField
+            const Field<Type>& completeField,
+            const bool preserveCouples=true
         );
 
         //- Construct from components and separate fields to slice for the
@@ -122,7 +124,8 @@ public:
             const Mesh&,
             const dimensionSet&,
             const Field<Type>& completeIField,
-            const Field<Type>& completeBField
+            const Field<Type>& completeBField,
+            const bool preserveCouples=true
         );
 
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index e9e012a09d07bab27bd0111a9a14e806365e59e3..daded730b9ba9ba3a76fb7667383275c15711dfb 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -55,6 +55,7 @@ void Foam::MULES::explicitSolve
     Info<< "MULES: Solving for " << psi.name() << endl;
 
     const fvMesh& mesh = psi.mesh();
+    psi.correctBoundaryConditions();
 
     surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
 
@@ -76,7 +77,8 @@ void Foam::MULES::explicitSolve
         ),
         mesh,
         dimless,
-        allLambda
+        allLambda,
+        false   // Use slices for the couples
     );
 
     limiter
@@ -183,7 +185,8 @@ void Foam::MULES::implicitSolve
             ),
             mesh,
             dimless,
-            allCoLambda
+            allCoLambda,
+            false   // Use slices for the couples
         );
 
         CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
@@ -205,7 +208,8 @@ void Foam::MULES::implicitSolve
         ),
         mesh,
         dimless,
-        allLambda
+        allLambda,
+        false   // Use slices for the couples
     );
 
     linear<scalar> CDs(mesh);
@@ -347,7 +351,8 @@ void Foam::MULES::limiter
         ),
         mesh,
         dimless,
-        allLambda
+        allLambda,
+        false   // Use slices for the couples
     );
 
     scalarField& lambdaIf = lambda;
@@ -406,6 +411,9 @@ void Foam::MULES::limiter
             {
                 label pfCelli = pFaceCells[pFacei];
 
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+
                 psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
                 psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
             }