diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 0390307f2..9633fff7a 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -139,13 +139,15 @@ void Foam::MRFZone::setMRFFaces()
     labelList nIncludedFaces(patches.size(), 0);
     labelList nExcludedFaces(patches.size(), 0);
 
+    register label faceI;
+
     forAll (patches, patchi)
     {
         const polyPatch& pp = patches[patchi];
 
         forAll (pp, patchFaceI)
         {
-            label faceI = pp.start() + patchFaceI;
+            faceI = pp.start() + patchFaceI;
 
             if (faceType[faceI] == 1)
             {
@@ -174,7 +176,7 @@ void Foam::MRFZone::setMRFFaces()
 
         forAll (pp, patchFaceI)
         {
-            label faceI = pp.start() + patchFaceI;
+            faceI = pp.start() + patchFaceI;
 
             if (faceType[faceI] == 1)
             {
@@ -187,7 +189,6 @@ void Foam::MRFZone::setMRFFaces()
         }
     }
 
-
     if (debug)
     {
         faceSet internalFaces
@@ -489,42 +490,51 @@ void Foam::MRFZone::absoluteFlux
 
 void Foam::MRFZone::meshPhi
 (
-    surfaceVectorField& phi
+    surfaceScalarField& phi
 ) const
 {
     const surfaceVectorField& Cf = mesh_.Cf();
+    const surfaceVectorField& Sf = mesh_.Sf();
 
     const vector& origin = origin_.value();
     const vector rotVel = Omega();
 
     // Internal faces
+    const vectorField& CfIn = Cf.internalField();
+    const vectorField& SfIn = Sf.internalField();
+
+    register label faceI, patchFaceI;
+
     forAll (internalFaces_, i)
     {
-        label facei = internalFaces_[i];
-        phi[facei] = (rotVel ^ (Cf[facei] - origin));
+        faceI = internalFaces_[i];
+        phi[faceI] = SfIn[faceI] & (rotVel ^ (CfIn[faceI] - origin));
     }
 
     // Included patches
-    forAll (includedFaces_, patchi)
-    {
-        forAll (includedFaces_[patchi], i)
-        {
-            label patchFaceI = includedFaces_[patchi][i];
 
-            phi.boundaryField()[patchi][patchFaceI] =
-                (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin));
+    forAll (includedFaces_, patchI)
+    {
+        forAll (includedFaces_[patchI], i)
+        {
+            patchFaceI = includedFaces_[patchI][i];
+
+            phi.boundaryField()[patchI][patchFaceI] =
+                Sf.boundaryField()[patchI][patchFaceI]
+              & (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin));
         }
     }
 
     // Excluded patches
-    forAll (excludedFaces_, patchi)
+    forAll (excludedFaces_, patchI)
     {
-        forAll (excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchI], i)
         {
-            label patchFaceI = excludedFaces_[patchi][i];
+            patchFaceI = excludedFaces_[patchI][i];
 
-            phi.boundaryField()[patchi][patchFaceI] =
-                (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin));
+            phi.boundaryField()[patchI][patchFaceI] =
+                Sf.boundaryField()[patchI][patchFaceI]
+              & (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin));
         }
     }
 }
@@ -534,21 +544,23 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
     const vector& origin = origin_.value();
     const vector rotVel = Omega();
 
+    register label patchFaceI;
+
     // Included patches
-    forAll (includedFaces_, patchi)
+    forAll (includedFaces_, patchI)
     {
-        const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
+        const vectorField& patchC = mesh_.Cf().boundaryField()[patchI];
 
-        vectorField pfld(U.boundaryField()[patchi]);
+        vectorField pfld(U.boundaryField()[patchI]);
 
-        forAll (includedFaces_[patchi], i)
+        forAll (includedFaces_[patchI], i)
         {
-            label patchFaceI = includedFaces_[patchi][i];
+            patchFaceI = includedFaces_[patchI][i];
 
             pfld[patchFaceI] = (rotVel ^ (patchC[patchFaceI] - origin));
         }
 
-        U.boundaryField()[patchi] == pfld;
+        U.boundaryField()[patchI] == pfld;
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index 954bc85b1..0121eed2a 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
@@ -221,7 +221,7 @@ public:
         ) const;
 
         //- Compute the pseudo mesh phi of the MRF region
-        void meshPhi(surfaceVectorField& phi) const;
+        void meshPhi(surfaceScalarField& phi) const;
 
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
index 0ae07364a..a6eabcbfa 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -44,34 +44,42 @@ void Foam::MRFZone::relativeRhoFlux
     const vector rotVel = Omega();
 
     // Internal faces
+    const vectorField& CfIn = Cf.internalField();
+    const vectorField& SfIn = Sf.internalField();
+
+    register label faceI, patchFaceI;
+
     forAll (internalFaces_, i)
     {
-        label facei = internalFaces_[i];
-        phi[facei] -= rho[facei]*(rotVel ^ (Cf[facei] - origin)) & Sf[facei];
+        faceI = internalFaces_[i];
+
+        phi[faceI] -=
+            rho[faceI]*(rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI];
     }
 
-    // Included patches
-    forAll (includedFaces_, patchi)
+    // Included patches: reset the flux to exactly zero to avoid
+    // round-off issues
+    forAll (includedFaces_, patchI)
     {
-        forAll (includedFaces_[patchi], i)
+        forAll (includedFaces_[patchI], i)
         {
-            label patchFacei = includedFaces_[patchi][i];
+            patchFaceI = includedFaces_[patchI][i];
 
-            phi.boundaryField()[patchi][patchFacei] = 0.0;
+            phi.boundaryField()[patchI][patchFaceI] = 0.0;
         }
     }
 
     // Excluded patches
-    forAll (excludedFaces_, patchi)
+    forAll (excludedFaces_, patchI)
     {
-        forAll (excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchI], i)
         {
-            label patchFacei = excludedFaces_[patchi][i];
+            patchFaceI = excludedFaces_[patchI][i];
 
-            phi.boundaryField()[patchi][patchFacei] -=
-                rho.boundaryField()[patchi][patchFacei]
-               *(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
+            phi.boundaryField()[patchI][patchFaceI] -=
+                rho.boundaryField()[patchI][patchFaceI]
+               *(rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
+              & Sf.boundaryField()[patchI][patchFaceI];
         }
     }
 }
@@ -91,35 +99,41 @@ void Foam::MRFZone::absoluteRhoFlux
     const vector rotVel = Omega();
 
     // Internal faces
+    const vectorField& CfIn = Cf.internalField();
+    const vectorField& SfIn = Sf.internalField();
+
+    register label faceI, patchFaceI;
+
     forAll (internalFaces_, i)
     {
-        label facei = internalFaces_[i];
-        phi[facei] += (rotVel ^ (Cf[facei] - origin)) & Sf[facei];
+        faceI = internalFaces_[i];
+        
+        phi[faceI] += (rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI];
     }
 
     // Included patches
-    forAll (includedFaces_, patchi)
+    forAll (includedFaces_, patchI)
     {
-        forAll (includedFaces_[patchi], i)
+        forAll (includedFaces_[patchI], i)
         {
-            label patchFacei = includedFaces_[patchi][i];
+            patchFaceI = includedFaces_[patchI][i];
 
-            phi.boundaryField()[patchi][patchFacei] +=
-                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
+            phi.boundaryField()[patchI][patchFaceI] +=
+                (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
+              & Sf.boundaryField()[patchI][patchFaceI];
         }
     }
 
     // Excluded patches
-    forAll (excludedFaces_, patchi)
+    forAll (excludedFaces_, patchI)
     {
-        forAll (excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchI], i)
         {
-            label patchFacei = excludedFaces_[patchi][i];
+            patchFaceI = excludedFaces_[patchI][i];
 
-            phi.boundaryField()[patchi][patchFacei] +=
-                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
+            phi.boundaryField()[patchI][patchFaceI] +=
+                (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
+              & Sf.boundaryField()[patchI][patchFaceI];
         }
     }
 }
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
index e5f29c690..f9cb77151 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
@@ -86,11 +86,11 @@ Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::fluxCorrection() const
 }
 
 
-Foam::tmp<Foam::surfaceVectorField> Foam::MRFZones::meshPhi() const
+Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::meshPhi() const
 {
-    tmp<surfaceVectorField> tMRFZonesPhiCorr
+    tmp<surfaceScalarField> tMRFZonesPhiCorr
     (
-        new surfaceVectorField
+        new surfaceScalarField
         (
             IOobject
             (
@@ -101,10 +101,10 @@ Foam::tmp<Foam::surfaceVectorField> Foam::MRFZones::meshPhi() const
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedVector("zero", dimVelocity, vector::zero)
+            dimensionedScalar("zero", dimVolume/dimTime, 0)
         )
     );
-    surfaceVectorField& MRFZonesPhiCorr = tMRFZonesPhiCorr();
+    surfaceScalarField& MRFZonesPhiCorr = tMRFZonesPhiCorr();
 
     forAll (*this, i)
     {
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
index ff683525d..082fc61f2 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
@@ -81,7 +81,7 @@ public:
         tmp<surfaceScalarField> fluxCorrection() const;
 
         //- Return pseudo mesh flux
-        tmp<surfaceVectorField> meshPhi() const;
+        tmp<surfaceScalarField> meshPhi() const;
 
         // Incompressible MRF
 
@@ -134,7 +134,7 @@ public:
         void absoluteVelocity(volVectorField& Urel) const;
 
 
-        //- Correct the boundary velocity for the roation of the MRF region
+        //- Correct the boundary velocity for the rotation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 
         //- Compute source term for volScalarFields