diff --git a/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H b/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H
index 5ba812f51..918bafa75 100644
--- a/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H
+++ b/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H
@@ -15,7 +15,7 @@
         // Warning.  Possible problem with a zone which is both MRF and porous
         // The solution is to do the loop below everywhere
         // Currently only inserting zones for efficiency.  HJ, 18/Nov/2017
-        register label cellI;
+        label cellI;
 
         forAll (mrfZones, mrfZoneI)
         {
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C
index d7737b890..845e97121 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.C
+++ b/src/ODE/ODESolvers/KRR4/KRR4.C
@@ -101,11 +101,11 @@ void Foam::KRR4::solve
 
     scalar h = hTry;
 
-    for (register label jtry=0; jtry<maxtry; jtry++)
+    for (label jtry=0; jtry<maxtry; jtry++)
     {
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
-            for (register label j=0; j<nEqns; j++)
+            for (label j=0; j<nEqns; j++)
             {
                 a_[i][j] = -dfdy_[i][j];
             }
@@ -115,14 +115,14 @@ void Foam::KRR4::solve
 
         simpleMatrix<scalar>::LUDecompose(a_, pivotIndices_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
         }
 
         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g1_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             y[i] = yTemp_[i] + a21*g1_[i];
         }
@@ -130,14 +130,14 @@ void Foam::KRR4::solve
         x = xTemp + a2X*h;
         ode_.derivatives(x, y, dydx_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
         }
 
         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g2_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
         }
@@ -145,14 +145,14 @@ void Foam::KRR4::solve
         x = xTemp + a3X*h;
         ode_.derivatives(x, y, dydx_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
         }
 
         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g3_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
                 + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
@@ -160,7 +160,7 @@ void Foam::KRR4::solve
 
         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g4_);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
             yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
@@ -176,7 +176,7 @@ void Foam::KRR4::solve
         }
 
         scalar maxErr = 0.0;
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
         }
diff --git a/src/ODE/ODESolvers/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C
index 2def00d47..5afc6fa26 100644
--- a/src/ODE/ODESolvers/RK/RK.C
+++ b/src/ODE/ODESolvers/RK/RK.C
@@ -166,7 +166,7 @@ void Foam::RK::solve
         solve(x, y, dydx, h, yTemp2_, yErr_);
 
         maxErr = 0.0;
-        for (register label i=0; i<ode_.nEqns(); i++)
+        for (label i=0; i<ode_.nEqns(); i++)
         {
             maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
         }
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index 5075d370e..f7f9c2ff5 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -95,14 +95,14 @@ void Foam::SIBS::solve
         scalar eps1 = safe1*eps;
         a_[0] = nSeq_[0] + 1;
 
-        for (register label k=0; k<kMaxX_; k++)
+        for (label k=0; k<kMaxX_; k++)
         {
             a_[k + 1] = a_[k] + nSeq_[k + 1];
         }
 
-        for (register label iq = 1; iq<kMaxX_; iq++)
+        for (label iq = 1; iq<kMaxX_; iq++)
         {
-            for (register label k=0; k<iq; k++)
+            for (label k=0; k<iq; k++)
             {
                 alpha_[k][iq] =
                     pow(eps1, (a_[k + 1] - a_[iq + 1])
@@ -113,7 +113,7 @@ void Foam::SIBS::solve
         epsOld_ = eps;
         a_[0] += nEqns;
 
-        for (register label k=0; k<kMaxX_; k++)
+        for (label k=0; k<kMaxX_; k++)
         {
             a_[k + 1] = a_[k] + nSeq_[k + 1];
         }
@@ -168,7 +168,7 @@ void Foam::SIBS::solve
             if (k != 0)
             {
                 maxErr = SMALL;
-                for (register label i=0; i<nEqns; i++)
+                for (label i=0; i<nEqns; i++)
                 {
                     maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
                 }
@@ -225,7 +225,7 @@ void Foam::SIBS::solve
     scalar wrkmin = GREAT;
     scalar scale = 1.0;
 
-    for (register label kk=0; kk<=km; kk++)
+    for (label kk=0; kk<=km; kk++)
     {
         scalar fact = max(err_[kk], scaleMX);
         scalar work = fact*a_[kk + 1];
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 28494f101..32bd20c35 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -44,9 +44,9 @@ void Foam::SIBS::SIMPR
     const label nEqns = ode_.nEqns();
 
     scalarSquareMatrix a(nEqns);
-    for (register label i=0; i<nEqns; i++)
+    for (label i=0; i<nEqns; i++)
     {
-        for (register label j=0; j<nEqns; j++)
+        for (label j=0; j<nEqns; j++)
         {
             a[i][j] = -h*dfdy[i][j];
         }
@@ -56,7 +56,7 @@ void Foam::SIBS::SIMPR
     labelList pivotIndices(nEqns);
     scalarSquareMatrix::LUDecompose(a, pivotIndices);
 
-    for (register label i=0; i<nEqns; i++)
+    for (label i=0; i<nEqns; i++)
     {
         yEnd[i] = h*(dydx[i] + h*dfdx[i]);
     }
@@ -66,7 +66,7 @@ void Foam::SIBS::SIMPR
     scalarField del(yEnd);
     scalarField ytemp(nEqns);
 
-    for (register label i=0; i<nEqns; i++)
+    for (label i=0; i<nEqns; i++)
     {
         ytemp[i] = y[i] + del[i];
     }
@@ -75,16 +75,16 @@ void Foam::SIBS::SIMPR
 
     ode_.derivatives(x, ytemp, yEnd);
 
-    for (register label nn=2; nn<=nSteps; nn++)
+    for (label nn=2; nn<=nSteps; nn++)
     {
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             yEnd[i] = h*yEnd[i] - del[i];
         }
 
         simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
 
-        for (register label i=0; i<nEqns; i++)
+        for (label i=0; i<nEqns; i++)
         {
             ytemp[i] += (del[i] += 2.0*yEnd[i]);
         }
@@ -93,14 +93,14 @@ void Foam::SIBS::SIMPR
 
         ode_.derivatives(x, ytemp, yEnd);
     }
-    for (register label i=0; i<nEqns; i++)
+    for (label i=0; i<nEqns; i++)
     {
         yEnd[i] = h*yEnd[i] - del[i];
     }
 
     simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
 
-    for (register label i=0; i<nEqns; i++)
+    for (label i=0; i<nEqns; i++)
     {
         yEnd[i] += ytemp[i];
     }
diff --git a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C b/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
index c016596b6..3c5bec63a 100644
--- a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
+++ b/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
@@ -42,14 +42,14 @@ void Foam::SIBS::polyExtrapolate
 
     x[iest] = xest;
 
-    for (register label j=0; j<n; j++)
+    for (label j=0; j<n; j++)
     {
         dy[j] = yz[j] = yest[j];
     }
 
     if (iest == 0)
     {
-        for (register label j=0; j<n; j++)
+        for (label j=0; j<n; j++)
         {
             d[j][0] = yest[j];
         }
@@ -58,13 +58,13 @@ void Foam::SIBS::polyExtrapolate
     {
         scalarField c(yest);
 
-        for (register label k1=0; k1<iest; k1++)
+        for (label k1=0; k1<iest; k1++)
         {
             scalar delta = 1.0/(x[iest - k1 - 1] - xest);
             scalar f1 = xest*delta;
             scalar f2 = x[iest - k1 - 1]*delta;
 
-            for (register label j=0; j<n; j++)
+            for (label j=0; j<n; j++)
             {
                 scalar q = d[j][k1];
                 d[j][k1] = dy[j];
@@ -75,7 +75,7 @@ void Foam::SIBS::polyExtrapolate
             }
         }
 
-        for (register label j=0; j<n; j++)
+        for (label j=0; j<n; j++)
         {
             d[j][iest] = dy[j];
         }
diff --git a/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C b/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C
index 342c00beb..11c6129da 100644
--- a/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C
+++ b/src/coupledMatrix/coupledLduPrecon/GaussSeidelPrecon/coupledGaussSeidelPrecon.C
@@ -69,7 +69,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweep
     const label nRows = x.size();
     label fStart, fEnd;
 
-    for (register label rowI = 0; rowI < nRows; rowI++)
+    for (label rowI = 0; rowI < nRows; rowI++)
     {
         // lRow is equal to rowI
         scalar& curX = x[rowI];
@@ -82,7 +82,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweep
         fEnd = ownStartAddr[rowI + 1];
 
         // Accumulate the owner product side
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             curX -= upper[curCoeff]*x[upperAddr[curCoeff]];
         }
@@ -91,7 +91,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweep
         curX /= diag[rowI];
 
         // Distribute the neighbour side using current x
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             bPrime[upperAddr[curCoeff]] -= lower[curCoeff]*curX;
         }
@@ -116,7 +116,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweep
     const label nRows = x.size();
     label fStart, fEnd;
 
-    for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+    for (label rowI = nRows - 1; rowI >= 0; rowI--)
     {
         // lRow is equal to rowI
         scalar& curX = x[rowI];
@@ -129,7 +129,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweep
         fEnd = ownStartAddr[rowI + 1];
 
         // Accumulate the owner product side
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             curX -= upper[curCoeff]*x[upperAddr[curCoeff]];
         }
@@ -138,7 +138,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweep
         curX /= diag[rowI];
 
         // Distribute the neighbour side using current x
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             bPrime[upperAddr[curCoeff]] -= lower[curCoeff]*curX;
         }
@@ -163,7 +163,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweepTranspose
     const label nRows = x.size();
     label fStart, fEnd;
 
-    for (register label rowI = 0; rowI < nRows; rowI++)
+    for (label rowI = 0; rowI < nRows; rowI++)
     {
         // lRow is equal to rowI
         scalar& curX = x[rowI];
@@ -176,7 +176,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweepTranspose
         fEnd = ownStartAddr[rowI + 1];
 
         // Accumulate the owner product side
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             // Transpose multiplication.  HJ, 19/Jan/2009
             curX -= lower[curCoeff]*x[upperAddr[curCoeff]];
@@ -186,7 +186,7 @@ void Foam::coupledGaussSeidelPrecon::forwardSweepTranspose
         curX /= diag[rowI];
 
         // Distribute the neighbour side using current x
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             // Transpose multiplication.  HJ, 19/Jan/2009
             bPrime[upperAddr[curCoeff]] -= upper[curCoeff]*curX;
@@ -212,7 +212,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweepTranspose
     const label nRows = x.size();
     label fStart, fEnd;
 
-    for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+    for (label rowI = nRows - 1; rowI >= 0; rowI--)
     {
         // lRow is equal to rowI
         scalar& curX = x[rowI];
@@ -225,7 +225,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweepTranspose
         fEnd = ownStartAddr[rowI + 1];
 
         // Accumulate the owner product side
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             // Transpose multiplication.  HJ, 19/Jan/2009
             curX -= lower[curCoeff]*x[upperAddr[curCoeff]];
@@ -235,7 +235,7 @@ void Foam::coupledGaussSeidelPrecon::reverseSweepTranspose
         curX /= diag[rowI];
 
         // Distribute the neighbour side using current x
-        for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+        for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
         {
             // Transpose multiplication.  HJ, 19/Jan/2009
             bPrime[upperAddr[curCoeff]] -= upper[curCoeff]*curX;
diff --git a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index 6578c9ca5..e5658aafe 100644
--- a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -73,7 +73,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     // one extra cell each
     for (j = 0; j < fstProcessorGroup; j++)
     {
-        for (register label k = 0; k < jumpb; k++)
+        for (label k = 0; k < jumpb; k++)
         {
             processorGroup[ind++] = j;
         }
@@ -82,7 +82,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     // and now to the `normal' processor groups
     for (; j<nProcGroup; j++)
     {
-        for (register label k = 0; k < jump; k++)
+        for (label k = 0; k < jump; k++)
         {
             processorGroup[ind++] = j;
         }
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
index 6196b3ea6..055254420 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/enrichedPatch/enrichedPatchPointPoints.C
@@ -57,7 +57,7 @@ void Foam::enrichedPatch::calcPointPoints() const
 
     const faceList& lf = localFaces();
 
-    register bool found = false;
+    bool found = false;
 
     forAll (lf, faceI)
     {
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
index 1faebce4a..cb81181b6 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
@@ -96,8 +96,8 @@ leastSquaresFaGrad<Type>::grad
 
     forAll(own, edgei)
     {
-        register label ownEdgeI = own[edgei];
-        register label neiEdgeI = nei[edgei];
+        label ownEdgeI = own[edgei];
+        label neiEdgeI = nei[edgei];
 
         Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI];
 
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
index 412efbc0f..5a9723d11 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
@@ -207,7 +207,7 @@ edgeInterpolationScheme<Type>::interpolate
 
     Field<Type>& sfi = sf.internalField();
 
-    for (register label fi=0; fi<P.size(); fi++)
+    for (label fi=0; fi<P.size(); fi++)
     {
         // ZT, 22/Apr/2003
         const tensorField& curT = mesh.edgeTransformTensors()[fi];
@@ -328,7 +328,7 @@ edgeInterpolationScheme<Type>::interpolate
 
     Field<Type>& sfi = sf.internalField();
 
-    for (register label eI = 0; eI < P.size(); eI++)
+    for (label eI = 0; eI < P.size(); eI++)
     {
         // ZT, 22/Apr/2003
         const tensorField& curT = mesh.edgeTransformTensors()[eI];
@@ -447,7 +447,7 @@ edgeInterpolationScheme<Type>::euclidianInterpolate
 
     Field<Type>& sfi = sf.internalField();
 
-    for (register label eI = 0; eI < P.size(); eI++)
+    for (label eI = 0; eI < P.size(); eI++)
     {
         sfi[eI] = lambda[eI]*vfi[P[eI]] + (1 - lambda[eI])*vfi[N[eI]];
     }
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index b357cc314..537c5fc5e 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -147,7 +147,7 @@ void Foam::MRFZone::setMRFFaces()
     labelList nIncludedFaces(patches.size(), 0);
     labelList nExcludedFaces(patches.size(), 0);
 
-    register label faceI;
+    label faceI;
 
     forAll (patches, patchi)
     {
@@ -367,7 +367,7 @@ void Foam::MRFZone::calcMeshVelocity() const
         );
 
     // Calculate mesh velocity for all moving faces
-    register label faceI, patchFaceI;
+    label faceI, patchFaceI;
 
     scalarField& meshVelIn = meshVel.internalField();
 
@@ -714,7 +714,7 @@ void Foam::MRFZone::meshPhi
 {
     const surfaceScalarField& meshVel = meshVelocity();
 
-    register label faceI, patchFaceI;
+    label faceI, patchFaceI;
 
     scalarField& phiIn = phi.internalField();
     const scalarField& meshVelIn = meshVel.internalField();
@@ -757,7 +757,7 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
     const vector& origin = origin_.value();
     const vector rotVel = Omega();
 
-    register label patchFaceI;
+    label patchFaceI;
 
     // Included patches
     forAll (includedFaces_, patchI)
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
index 3e5f178a1..be1817d5f 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -41,7 +41,7 @@ void Foam::MRFZone::relativeRhoFlux
     // HJ, 6/Jun/2017
     const surfaceScalarField& meshVel = meshVelocity();
 
-    register label faceI, patchFaceI;
+    label faceI, patchFaceI;
 
     // Internal faces
     scalarField& phiIn = phi.internalField();
@@ -94,7 +94,7 @@ void Foam::MRFZone::absoluteRhoFlux
     // HJ, 6/Jun/2017
     const surfaceScalarField& meshVel = meshVelocity();
 
-    register label faceI, patchFaceI;
+    label faceI, patchFaceI;
 
     // Internal faces
     scalarField& phiIn = phi.internalField();
diff --git a/src/finiteVolume/finiteVolume/divSchemes/leastSquaresDivScheme/vectorLeastSquaresDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/leastSquaresDivScheme/vectorLeastSquaresDivScheme.C
index 138dfcbd8..e6c73042f 100644
--- a/src/finiteVolume/finiteVolume/divSchemes/leastSquaresDivScheme/vectorLeastSquaresDivScheme.C
+++ b/src/finiteVolume/finiteVolume/divSchemes/leastSquaresDivScheme/vectorLeastSquaresDivScheme.C
@@ -98,7 +98,7 @@ tmp<BlockLduSystem<vector, scalar> > leastSquaresDivScheme<vector>::fvmUDiv
     const vectorField& neiLsIn = neiLs.internalField();
 //     const scalarField& wIn = weights.internalField();
 
-    register label owner, neighbour;
+    label owner, neighbour;
 
     forAll (nei, faceI)
     {
@@ -231,8 +231,8 @@ tmp<BlockLduSystem<vector, scalar> > leastSquaresDivScheme<vector>::fvmUDiv
 
     forAll (nei, faceI)
     {
-        register label owner = own[faceI];
-        register label neighbour = nei[faceI];
+        label owner = own[faceI];
+        label neighbour = nei[faceI];
 
         u[faceI] = cellVIn[owner]*ownLsIn[faceI]*fluxIn[faceI];
         l[faceI] = cellVIn[neighbour]*neiLsIn[faceI]*fluxIn[faceI];
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
index 4670f7cd6..e4dc336c1 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
@@ -104,7 +104,7 @@ leastSquaresGrad<Type>::calcGrad
     const vectorField& ownLsIn = ownLs.internalField();
     const vectorField& neiLsIn = neiLs.internalField();
 
-    register label ownFaceI, neiFaceI;
+    label ownFaceI, neiFaceI;
 
     forAll (own, facei)
     {
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/scalarLeastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/scalarLeastSquaresGrad.C
index 144830693..ca202189a 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/scalarLeastSquaresGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/scalarLeastSquaresGrad.C
@@ -94,7 +94,7 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
     const vectorField& ownLsIn = ownLs.internalField();
     const vectorField& neiLsIn = neiLs.internalField();
 
-    register label owner, neighbour;
+    label owner, neighbour;
 
     forAll (nei, faceI)
     {
diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C
index 8cf82a0db..9e2a521fa 100644
--- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C
@@ -941,7 +941,7 @@ void Foam::fvBlockMatrix<Type>::insertAdjointConvection
     // (boundary) contributions - need to check. VV, 7/Apr/2016.
 
     // Loop through faces
-    register label own, nei;
+    label own, nei;
     forAll (neighbour, faceI)
     {
         own = owner[faceI];
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C
index 611900f37..e181bd205 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C
@@ -342,7 +342,7 @@ void Foam::ggiFvPatch::expandCrMatrixToZone(crMatrix& patchP) const
         List<labelField> cols(patchCrAddr.nRows());
         List<scalarField> coeffs(patchCrAddr.nRows());
 
-        for (register label faceI = 0; faceI < patchCrAddr.nRows(); faceI++)
+        for (label faceI = 0; faceI < patchCrAddr.nRows(); faceI++)
         {
             // Unpack row
             const label rowStart = patchRowStart[faceI];
@@ -354,7 +354,7 @@ void Foam::ggiFvPatch::expandCrMatrixToZone(crMatrix& patchP) const
             coeffs[faceI].setSize(rowLength);
             scalarField& curCoeffs = coeffs[faceI];
 
-            for (register label coeffI = 0; coeffI < rowLength; coeffI++)
+            for (label coeffI = 0; coeffI < rowLength; coeffI++)
             {
                 curCols[coeffI] = patchCol[rowStart + coeffI];
                 curCoeffs[coeffI] = patchCoeff[rowStart + coeffI];
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H
index 1b6348475..02e337ec5 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/biLinearFit/biLinearFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H
index 5b895e598..5383eafc0 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubicUpwindFit/cubicUpwindFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H
index efdb9d5d3..1b8a49c68 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H
index b49d7303a..cd02b1a84 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H
index 8ebaddca3..0c2139836 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H
@@ -72,7 +72,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H
index 8d81d7af2..334303198 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearUpwindFit/quadraticLinearUpwindFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H
index daa9d6d91..797cd6248 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticUpwindFit/quadraticUpwindFitPolynomial.H
@@ -70,7 +70,7 @@ public:
             const direction dim
         )
         {
-            register label curIdx = 0;
+            label curIdx = 0;
 
             coeffs[curIdx++] = weight;
             coeffs[curIdx++] = weight*d.x();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index b8b75c549..4b0e01828 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -208,7 +208,7 @@ surfaceInterpolationScheme<Type>::interpolate
 
     Field<Type>& sfi = sf.internalField();
 
-    for (register label fi=0; fi<P.size(); fi++)
+    for (label fi=0; fi<P.size(); fi++)
     {
         sfi[fi] = lambda[fi]*vfi[P[fi]] + y[fi]*vfi[N[fi]];
     }
@@ -283,7 +283,7 @@ surfaceInterpolationScheme<Type>::interpolate
 
     Field<Type>& sfi = sf.internalField();
 
-    for (register label fi=0; fi<P.size(); fi++)
+    for (label fi=0; fi<P.size(); fi++)
     {
         sfi[fi] = lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]];
     }
diff --git a/src/foam/containers/Lists/FixedList/FixedListI.H b/src/foam/containers/Lists/FixedList/FixedListI.H
index 76140a128..31b8c7019 100644
--- a/src/foam/containers/Lists/FixedList/FixedListI.H
+++ b/src/foam/containers/Lists/FixedList/FixedListI.H
@@ -37,7 +37,7 @@ inline Foam::FixedList<T, Size>::FixedList()
 template<class T, unsigned Size>
 inline Foam::FixedList<T, Size>::FixedList(const T v[Size])
 {
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = v[i];
     }
@@ -47,7 +47,7 @@ inline Foam::FixedList<T, Size>::FixedList(const T v[Size])
 template<class T, unsigned Size>
 inline Foam::FixedList<T, Size>::FixedList(const T& t)
 {
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = t;
     }
@@ -59,7 +59,7 @@ inline Foam::FixedList<T, Size>::FixedList(const UList<T>& lst)
 {
     checkSize(lst.size());
 
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = lst[i];
     }
@@ -71,7 +71,7 @@ inline Foam::FixedList<T, Size>::FixedList(const SLList<T>& lst)
 {
     checkSize(lst.size());
 
-    register label i = 0;
+    label i = 0;
     for
     (
         typename SLList<T>::const_iterator iter = lst.begin();
@@ -87,7 +87,7 @@ inline Foam::FixedList<T, Size>::FixedList(const SLList<T>& lst)
 template<class T, unsigned Size>
 inline Foam::FixedList<T, Size>::FixedList(const FixedList<T, Size>& lst)
 {
-    for (register unsigned i = 0; i < Size; i++)
+    for (unsigned i = 0; i < Size; i++)
     {
         v_[i] = lst[i];
     }
@@ -184,7 +184,7 @@ inline void Foam::FixedList<T, Size>::setSize(const label s)
 template<class T, unsigned Size>
 inline void Foam::FixedList<T, Size>::transfer(const FixedList<T, Size>& lst)
 {
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = lst[i];
     }
@@ -234,7 +234,7 @@ inline const T& Foam::FixedList<T, Size>::operator[](const label i) const
 template<class T, unsigned Size>
 inline void Foam::FixedList<T, Size>::operator=(const T lst[Size])
 {
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = lst[i];
     }
@@ -245,7 +245,7 @@ inline void Foam::FixedList<T, Size>::operator=(const UList<T>& lst)
 {
     checkSize(lst.size());
 
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i=0; i<Size; i++)
     {
         v_[i] = lst[i];
     }
@@ -256,7 +256,7 @@ inline void Foam::FixedList<T, Size>::operator=(const SLList<T>& lst)
 {
     checkSize(lst.size());
 
-    register label i = 0;
+    label i = 0;
     for
     (
         typename SLList<T>::const_iterator iter = lst.begin();
@@ -271,7 +271,7 @@ inline void Foam::FixedList<T, Size>::operator=(const SLList<T>& lst)
 template<class T, unsigned Size>
 inline void Foam::FixedList<T, Size>::operator=(const T& t)
 {
-    for (register unsigned i=0; i<Size; i++)
+    for (unsigned i = 0; i < Size; i++)
     {
         v_[i] = t;
     }
@@ -415,7 +415,7 @@ inline unsigned Foam::FixedList<T, Size>::Hash<HashT>::operator()
         // hash incrementally
         unsigned val = seed;
 
-        for (register unsigned i=0; i<Size; i++)
+        for (unsigned i = 0; i<Size; i++)
         {
             val = HashT()(lst[i], val);
         }
diff --git a/src/foam/containers/Lists/FixedList/FixedListIO.C b/src/foam/containers/Lists/FixedList/FixedListIO.C
index 1d855a628..0fc75e0d4 100644
--- a/src/foam/containers/Lists/FixedList/FixedListIO.C
+++ b/src/foam/containers/Lists/FixedList/FixedListIO.C
@@ -85,7 +85,7 @@ Foam::Istream& Foam::operator>>(Foam::Istream& is, FixedList<T, Size>& L)
 
         if (delimiter == token::BEGIN_LIST)
         {
-            for (register unsigned i=0; i<Size; i++)
+            for (unsigned i = 0; i < Size; i++)
             {
                 is >> L[i];
 
@@ -107,7 +107,7 @@ Foam::Istream& Foam::operator>>(Foam::Istream& is, FixedList<T, Size>& L)
                 "reading the single entry"
             );
 
-            for (register unsigned i=0; i<Size; i++)
+            for (unsigned i = 0; i < Size; i++)
             {
                 L[i] = element;
             }
diff --git a/src/foam/containers/Lists/List/List.C b/src/foam/containers/Lists/List/List.C
index ec98d2e44..5c1d5a256 100644
--- a/src/foam/containers/Lists/List/List.C
+++ b/src/foam/containers/Lists/List/List.C
@@ -354,7 +354,7 @@ void Foam::List<T>::setSize(const label newSize)
 
             if (this->size_)
             {
-                register label i = min(this->size_, newSize);
+                label i = min(this->size_, newSize);
 
 #               ifdef USEMEMCPY
                 if (contiguous<T>())
@@ -364,8 +364,8 @@ void Foam::List<T>::setSize(const label newSize)
                 else
 #               endif
                 {
-                    register T* vv = &this->v_[i];
-                    register T* av = &nv[i];
+                    T* vv = &this->v_[i];
+                    T* av = &nv[i];
                     while (i--) *--av = *--vv;
                 }
             }
@@ -390,8 +390,8 @@ void Foam::List<T>::setSize(const label newSize, const T& a)
 
     if (newSize > oldSize)
     {
-        register label i = newSize - oldSize;
-        register T* vv = &this->v_[newSize];
+        label i = newSize - oldSize;
+        T* vv = &this->v_[newSize];
         while (i--) *--vv = a;
     }
 }
diff --git a/src/foam/containers/Lists/List/ListIO.C b/src/foam/containers/Lists/List/ListIO.C
index 9ad637395..44a13c786 100644
--- a/src/foam/containers/Lists/List/ListIO.C
+++ b/src/foam/containers/Lists/List/ListIO.C
@@ -80,7 +80,7 @@ Foam::Istream& Foam::operator>>(Istream& is, List<T>& L)
             {
                 if (delimiter == token::BEGIN_LIST)
                 {
-                    for (register label i=0; i<s; i++)
+                    for (label i=0; i<s; i++)
                     {
                         is >> L[i];
 
@@ -101,7 +101,7 @@ Foam::Istream& Foam::operator>>(Istream& is, List<T>& L)
                         "reading the single entry"
                     );
 
-                    for (register label i=0; i<s; i++)
+                    for (label i=0; i<s; i++)
                     {
                         L[i] = element;
                     }
diff --git a/src/foam/containers/Lists/List/ListLoopM.H b/src/foam/containers/Lists/List/ListLoopM.H
index cd931bf3a..410a5783f 100644
--- a/src/foam/containers/Lists/List/ListLoopM.H
+++ b/src/foam/containers/Lists/List/ListLoopM.H
@@ -37,8 +37,8 @@ Description
 // Element access looping using [] for vector machines
 
 #define List_FOR_ALL(f, i)                      \
-        register const label _n##i = (f).size();\
-        for (register label i=0; i<_n##i; i++)  \
+        const label _n##i = (f).size();         \
+        for (label i = 0; i < _n##i; i++)       \
         {
 
 #define List_END_FOR_ALL  }
@@ -56,7 +56,7 @@ Description
 // Pointer looping for scalar machines
 
 #define List_FOR_ALL(f, i)                      \
-        register label i = (f).size();          \
+        label i = (f).size();                   \
         while (i--)                             \
         {                                       \
 
@@ -65,10 +65,10 @@ Description
 #define List_ELEM(f, fp, i)  (*fp++)
 
 #define List_ACCESS(type, f, fp) \
-    register type* __restrict__ fp = (f).begin()
+    type* __restrict__ fp = (f).begin()
 
 #define List_CONST_ACCESS(type, f, fp) \
-    register const type* __restrict__ fp = (f).begin()
+    const type* __restrict__ fp = (f).begin()
 
 #endif
 
diff --git a/src/foam/containers/Lists/PtrList/PtrList.C b/src/foam/containers/Lists/PtrList/PtrList.C
index 11d35672f..87cbf5030 100644
--- a/src/foam/containers/Lists/PtrList/PtrList.C
+++ b/src/foam/containers/Lists/PtrList/PtrList.C
@@ -156,7 +156,7 @@ void Foam::PtrList<T>::setSize(const label newSize)
     }
     else if (newSize < oldSize)
     {
-        register label i;
+        label i;
         for (i=newSize; i<oldSize; i++)
         {
             if (ptrs_[i])
@@ -171,7 +171,7 @@ void Foam::PtrList<T>::setSize(const label newSize)
     {
         ptrs_.setSize(newSize);
 
-        register label i;
+        label i;
         for (i=oldSize; i<newSize; i++)
         {
             ptrs_[i] = nullptr;
diff --git a/src/foam/containers/Lists/UPtrList/UPtrList.C b/src/foam/containers/Lists/UPtrList/UPtrList.C
index 62cc13acf..37cdfb822 100644
--- a/src/foam/containers/Lists/UPtrList/UPtrList.C
+++ b/src/foam/containers/Lists/UPtrList/UPtrList.C
@@ -76,7 +76,7 @@ void Foam::UPtrList<T>::setSize(const label newSize)
     if (newSize > oldSize)
     {
         // set new elements to nullptr
-        for (register label i = oldSize; i < newSize; i++)
+        for (label i = oldSize; i < newSize; i++)
         {
             ptrs_[i] = nullptr;
         }
diff --git a/src/foam/db/IOstreams/Pstreams/IPstream.C b/src/foam/db/IOstreams/Pstreams/IPstream.C
index 593056d5b..b49c9a506 100644
--- a/src/foam/db/IOstreams/Pstreams/IPstream.C
+++ b/src/foam/db/IOstreams/Pstreams/IPstream.C
@@ -65,9 +65,9 @@ inline void Foam::IPstream::readFromBuffer
         bufPosition_ = align + ((bufPosition_ - 1) & ~(align - 1));
     }
 
-    register const char* bufPtr = &buf_[bufPosition_];
-    register char* dataPtr = reinterpret_cast<char*>(data);
-    register size_t i = count;
+    const char* bufPtr = &buf_[bufPosition_];
+    char* dataPtr = reinterpret_cast<char*>(data);
+    size_t i = count;
     while (i--) *dataPtr++ = *bufPtr++;
     bufPosition_ += count;
     checkEof();
diff --git a/src/foam/db/IOstreams/Pstreams/OPstream.C b/src/foam/db/IOstreams/Pstreams/OPstream.C
index a44cf34a2..3c0d4d39e 100644
--- a/src/foam/db/IOstreams/Pstreams/OPstream.C
+++ b/src/foam/db/IOstreams/Pstreams/OPstream.C
@@ -75,9 +75,9 @@ inline void Foam::OPstream::writeToBuffer
         enlargeBuffer(bufPosition_ - oldPos + count);
     }
 
-    register char* bufPtr = &buf_[bufPosition_];
-    register const char* dataPtr = reinterpret_cast<const char*>(data);
-    register size_t i = count;
+    char* bufPtr = &buf_[bufPosition_];
+    const char* dataPtr = reinterpret_cast<const char*>(data);
+    size_t i = count;
     while (i--) *bufPtr++ = *dataPtr++;
 
     bufPosition_ += count;
diff --git a/src/foam/db/IOstreams/Sstreams/OSstream.C b/src/foam/db/IOstreams/Sstreams/OSstream.C
index b76417108..352b7d0c2 100644
--- a/src/foam/db/IOstreams/Sstreams/OSstream.C
+++ b/src/foam/db/IOstreams/Sstreams/OSstream.C
@@ -80,10 +80,10 @@ Foam::Ostream& Foam::OSstream::write(const string& str)
 {
     os_ << token::BEGIN_STRING;
 
-    register int backslash = 0;
+    int backslash = 0;
     for (string::const_iterator iter = str.begin(); iter != str.end(); ++iter)
     {
-        register char c = *iter;
+        char c = *iter;
 
         if (c == '\\')
         {
@@ -131,7 +131,7 @@ Foam::Ostream& Foam::OSstream::writeQuoted
     {
         os_ << token::BEGIN_STRING;
 
-        register int backslash = 0;
+        int backslash = 0;
         for
         (
             string::const_iterator iter = str.begin();
@@ -139,7 +139,7 @@ Foam::Ostream& Foam::OSstream::writeQuoted
             ++iter
         )
         {
-            register char c = *iter;
+            char c = *iter;
 
             if (c == '\\')
             {
@@ -236,7 +236,7 @@ Foam::Ostream& Foam::OSstream::write(const char* buf, std::streamsize count)
 
 void Foam::OSstream::indent()
 {
-    for (register unsigned short i = 0; i < indentLevel_*indentSize_; i++)
+    for (unsigned short i = 0; i < indentLevel_*indentSize_; i++)
     {
         os_ << ' ';
     }
diff --git a/src/foam/db/IOstreams/Sstreams/readHexLabel.C b/src/foam/db/IOstreams/Sstreams/readHexLabel.C
index 0ec0c7a66..685d14416 100644
--- a/src/foam/db/IOstreams/Sstreams/readHexLabel.C
+++ b/src/foam/db/IOstreams/Sstreams/readHexLabel.C
@@ -43,7 +43,7 @@ Foam::label Foam::readHexLabel(ISstream& is)
     while (is.get(c) && isspace(c))
     {}
 
-    register label result = 0;
+    label result = 0;
     do
     {
         if (isspace(c) || c == 0) break;
diff --git a/src/foam/graph/curve/curve.C b/src/foam/graph/curve/curve.C
index 8d6aff4ea..1c0d4a67c 100644
--- a/src/foam/graph/curve/curve.C
+++ b/src/foam/graph/curve/curve.C
@@ -46,7 +46,7 @@ curve::curve(const curve& Curve, const label nFacets)
 {
     // Calculate curve length
     scalar curveLength=0;
-    register label i;
+    label i;
     for (i=0; i<Curve.size()-1; i++)
     {
         curveLength += distance(Curve[i], Curve[i+1]);
@@ -129,7 +129,7 @@ curve grad(const curve& Curve)
 {
     curve gradCurve(Curve);
 
-    register label i;
+    label i;
     for (i=1; i<Curve.size()-1; i++)
     {
         scalar deltaIm1 = Curve[i].x() - Curve[i-1].x();
diff --git a/src/foam/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/foam/matrices/LUscalarMatrix/LUscalarMatrix.C
index be29cbdfc..c0ce29e6e 100644
--- a/src/foam/matrices/LUscalarMatrix/LUscalarMatrix.C
+++ b/src/foam/matrices/LUscalarMatrix/LUscalarMatrix.C
@@ -137,15 +137,15 @@ void Foam::LUscalarMatrix::convert
     const scalar* __restrict__ upperPtr = ldum.upper().begin();
     const scalar* __restrict__ lowerPtr = ldum.lower().begin();
 
-    register const label nCells = ldum.diag().size();
-    register const label nFaces = ldum.upper().size();
+    const label nCells = ldum.diag().size();
+    const label nFaces = ldum.upper().size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         operator[](cell)[cell] = diagPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         label uCell = uPtr[face];
         label lCell = lPtr[face];
@@ -164,9 +164,9 @@ void Foam::LUscalarMatrix::convert
             const scalar* __restrict__ upperLowerPtr =
                 interfaceCoeffs[inti].begin();
 
-            register label inFaces = interface.faceCells().size()/2;
+            label inFaces = interface.faceCells().size()/2;
 
-            for (register label face=0; face<inFaces; face++)
+            for (label face=0; face<inFaces; face++)
             {
                 label uCell = ulPtr[face];
                 label lCell = ulPtr[face + inFaces];
@@ -206,16 +206,16 @@ void Foam::LUscalarMatrix::convert
         const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
         const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
 
-        register const label nCells = lduMatrixi.size();
-        register const label nFaces = lduMatrixi.upper_.size();
+        const label nCells = lduMatrixi.size();
+        const label nFaces = lduMatrixi.upper_.size();
 
-        for (register label cell=0; cell<nCells; cell++)
+        for (label cell=0; cell<nCells; cell++)
         {
             label globalCell = cell + offset;
             operator[](globalCell)[globalCell] = diagPtr[cell];
         }
 
-        for (register label face=0; face<nFaces; face++)
+        for (label face=0; face<nFaces; face++)
         {
             label uCell = uPtr[face] + offset;
             label lCell = lPtr[face] + offset;
@@ -238,9 +238,9 @@ void Foam::LUscalarMatrix::convert
                 const scalar* __restrict__ upperLowerPtr =
                     interface.coeffs_.begin();
 
-                register label inFaces = interface.faceCells_.size()/2;
+                label inFaces = interface.faceCells_.size()/2;
 
-                for (register label face=0; face<inFaces; face++)
+                for (label face=0; face<inFaces; face++)
                 {
                     label uCell = ulPtr[face] + offset;
                     label lCell = ulPtr[face + inFaces] + offset;
@@ -283,10 +283,10 @@ void Foam::LUscalarMatrix::convert
                 const scalar* __restrict__ upperPtr = interface.coeffs_.begin();
                 const scalar* __restrict__ lowerPtr = neiInterface.coeffs_.begin();
 
-                register label inFaces = interface.faceCells_.size();
+                label inFaces = interface.faceCells_.size();
                 label neiOffset = procOffsets_[interface.neighbProcNo_];
 
-                for (register label face=0; face<inFaces; face++)
+                for (label face=0; face<inFaces; face++)
                 {
                     label uCell = uPtr[face] + offset;
                     label lCell = lPtr[face] + neiOffset;
diff --git a/src/foam/matrices/Matrix/Matrix.C b/src/foam/matrices/Matrix/Matrix.C
index 7749ab974..053951a44 100644
--- a/src/foam/matrices/Matrix/Matrix.C
+++ b/src/foam/matrices/Matrix/Matrix.C
@@ -41,7 +41,7 @@ void Foam::Matrix<Form, Type>::allocate()
         v_ = new Type*[n_];
         v_[0] = new Type[n_*m_];
 
-        for (register label i=1; i<n_; i++)
+        for (label i=1; i<n_; i++)
         {
             v_[i] = v_[i-1] + m_;
         }
@@ -109,7 +109,7 @@ Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
 
         label nm = n_*m_;
 
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             v[i] = a;
         }
@@ -131,7 +131,7 @@ Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
         const Type* av = a.v_[0];
 
         label nm = n_*m_;
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             v[i] = av[i];
         }
@@ -175,9 +175,9 @@ Form Foam::Matrix<Form, Type>::T() const
     const Matrix<Form, Type>& A = *this;
     Form At(m(), n());
 
-    for (register label i=0; i<n(); i++)
+    for (label i=0; i<n(); i++)
     {
-        for (register label j=0; j<m(); j++)
+        for (label j=0; j<m(); j++)
         {
             At[j][i] = A[i][j];
         }
@@ -197,7 +197,7 @@ void Foam::Matrix<Form, Type>::operator=(const Type& t)
         Type* v = v_[0];
 
         label nm = n_*m_;
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             v[i] = t;
         }
@@ -232,7 +232,7 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
         const Type* av = a.v_[0];
 
         label nm = n_*m_;
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             v[i] = av[i];
         }
@@ -252,7 +252,7 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
         label curMaxI = 0;
         const Type* v = a[0];
 
-        for (register label i=1; i<nm; i++)
+        for (label i=1; i<nm; i++)
         {
             if (v[i] > v[curMaxI])
             {
@@ -284,7 +284,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
         label curMinI = 0;
         const Type* v = a[0];
 
-        for (register label i=1; i<nm; i++)
+        for (label i=1; i<nm; i++)
         {
             if (v[i] < v[curMinI])
             {
@@ -319,7 +319,7 @@ Form Foam::operator-(const Matrix<Form, Type>& a)
         const Type* av = a[0];
 
         label nm = a.n()*a.m();
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             nav[i] = -av[i];
         }
@@ -361,7 +361,7 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     const Type* bv = b[0];
 
     label nm = a.n()*a.m();
-    for (register label i=0; i<nm; i++)
+    for (label i=0; i<nm; i++)
     {
         abv[i] = av[i] + bv[i];
     }
@@ -402,7 +402,7 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     const Type* bv = b[0];
 
     label nm = a.n()*a.m();
-    for (register label i=0; i<nm; i++)
+    for (label i=0; i<nm; i++)
     {
         abv[i] = av[i] - bv[i];
     }
@@ -422,7 +422,7 @@ Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
         const Type* av = a[0];
 
         label nm = a.n()*a.m();
-        for (register label i=0; i<nm; i++)
+        for (label i=0; i<nm; i++)
         {
             sav[i] = s*av[i];
         }
diff --git a/src/foam/matrices/Matrix/MatrixIO.C b/src/foam/matrices/Matrix/MatrixIO.C
index b8ac56d61..b45bd86af 100644
--- a/src/foam/matrices/Matrix/MatrixIO.C
+++ b/src/foam/matrices/Matrix/MatrixIO.C
@@ -77,11 +77,11 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
                     label k = 0;
 
                     // loop over rows
-                    for (register label i=0; i<M.n(); i++)
+                    for (label i=0; i<M.n(); i++)
                     {
                         listDelimiter = is.readBeginList("MatrixRow");
 
-                        for (register label j=0; j<M.m(); j++)
+                        for (label j=0; j<M.m(); j++)
                         {
                             is >> v[k++];
 
@@ -106,7 +106,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
                         "reading the single entry"
                     );
 
-                    for (register label i=0; i<nm; i++)
+                    for (label i=0; i<nm; i++)
                     {
                         v[i] = element;
                     }
@@ -165,7 +165,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
             {
                 uniform = true;
 
-                for (register label i=0; i< nm; i++)
+                for (label i=0; i< nm; i++)
                 {
                     if (v[i] != v[0])
                     {
@@ -196,12 +196,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
                 label k = 0;
 
                 // loop over rows
-                for (register label i=0; i< M.n(); i++)
+                for (label i=0; i< M.n(); i++)
                 {
                     os  << nl << token::BEGIN_LIST;
 
                     // Write row
-                    for (register label j=0; j< M.m(); j++)
+                    for (label j=0; j< M.m(); j++)
                     {
                         if (j > 0) os << token::SPACE;
                         os << v[k++];
@@ -221,12 +221,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
                 label k = 0;
 
                 // loop over rows
-                for (register label i=0; i< M.n(); i++)
+                for (label i=0; i< M.n(); i++)
                 {
                     os  << nl << token::BEGIN_LIST;
 
                     // Write row
-                    for (register label j=0; j< M.m(); j++)
+                    for (label j=0; j< M.m(); j++)
                     {
                         os << nl << v[k++];
                     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C
index adbaa2747..5947420b6 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C
@@ -186,12 +186,12 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
     scalarField epsilonStrongCoeff(nRows, 0);
 
     // Select the strongest coefficient in each row
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(normDiag[i]);
 
         // Do lower triangle coefficient for the row first
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij =
                 mag(min(signDiag*normLower[losortAddr[jp]], 0));
@@ -203,7 +203,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         }
 
         // Do upper triangle coefficient for the row
-        for (register label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
+        for (label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
         {
             const scalar magAij = mag(min(signDiag*normUpper[ip], 0));
 
@@ -222,12 +222,12 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
     // positive and negative coefficients.  HJ, 28/Feb/2017
     labelList strongCoeffCounter(nRows, 0);
 
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(normDiag[i]);
 
         // Do lower triangle coefficient for the row first
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij =
                 mag(min(signDiag*normLower[losortAddr[jp]], 0));
@@ -239,7 +239,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         }
 
         // Do upper triangle coefficient for the row
-        for (register label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
+        for (label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
         {
             const scalar magAij = mag(min(signDiag*normUpper[ip], 0));
 
@@ -266,12 +266,12 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
     // Reset strongCoeffCounter for re-use
     strongCoeffCounter = 0;
 
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(normDiag[i]);
 
         // Do lower triangle coefficient for the row first
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij =
                 mag(min(signDiag*normLower[losortAddr[jp]], 0));
@@ -511,7 +511,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
     // Start row assembly
     pRow[0] = 0;
 
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         label rowCount = pRow[i];
 
@@ -520,7 +520,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         signDii = sign(normDiag[i]);
 
         // Do lower triangle coefficient for the row first
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             // Adjust sign of off-diag coeff
             Dij = sign(normDiag[i])*normLower[losortAddr[jp]];
@@ -532,7 +532,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         }
 
         // Do upper triangle coefficient for the row
-        for (register label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
+        for (label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
         {
             // Adjust sign of off-diag coeff
             Dij = signDii*normUpper[ip];
@@ -551,7 +551,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
 
             for
             (
-                register label sip = strongRow[i];
+                label sip = strongRow[i];
                 sip < strongRow[i + 1];
                 sip++
             )
@@ -566,7 +566,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
 
             for
             (
-                register label sip = strongRow[i];
+                label sip = strongRow[i];
                 sip < strongRow[i + 1];
                 sip++
             )
@@ -614,7 +614,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         const labelList& prolongationRow = Pptr_->crAddr().rowStart();
         const scalarField& prolongationCoeff = Pptr_->coeffs();
 
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             for
             (
@@ -1505,9 +1505,9 @@ void Foam::BlockMatrixSelection<Type>::restrictResidual
     // Multiply the residual with restriction weights to obtain the initial
     // coarse residual
 
-    for (register label i = 0; i < nCoarseEqns_; i++)
+    for (label i = 0; i < nCoarseEqns_; i++)
     {
-        for (register label k = rowR[i]; k < rowR[i + 1]; k++)
+        for (label k = rowR[i]; k < rowR[i + 1]; k++)
         {
             // Multiply each coeff in row of restriction with the corresponding
             // residual coefficient (col index of R is the same as row index
@@ -1546,9 +1546,9 @@ void Foam::BlockMatrixSelection<Type>::prolongateCorrection
     // Multiply the coarse level solution with prolongation and obtain fine
     // level solution
 
-    for (register label i = 0; i < sizeP; i++)
+    for (label i = 0; i < sizeP; i++)
     {
-        for (register label k = rowP[i]; k < rowP[i + 1]; k++)
+        for (label k = rowP[i]; k < rowP[i + 1]; k++)
         {
             // Multiply each coeff in row of prolongation with the corresponding
             // coarse residual coefficient (col index of prolongation must be
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
index f8781964a..95d0441c6 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
@@ -83,7 +83,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -92,7 +92,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -101,7 +101,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 // Use transpose upper coefficient
                 Ax[u[coeffI]] +=
@@ -117,7 +117,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -126,7 +126,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -135,7 +135,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
         {
             const squareTypeField& activeLower = Lower.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -148,7 +148,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
     {
         const scalarTypeField& activeUpper = Upper.asScalar();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -157,7 +157,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
     {
         const linearTypeField& activeUpper = Upper.asLinear();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -166,7 +166,7 @@ void Foam::BlockLduMatrix<Type>::AmulCore
     {
         const squareTypeField& activeUpper = Upper.asSquare();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -222,7 +222,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
     {
         const scalarTypeField& activeUpper = Upper.asScalar();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
         }
@@ -231,7 +231,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
     {
         const linearTypeField& activeUpper = Upper.asLinear();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
         }
@@ -240,7 +240,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
     {
         const squareTypeField& activeUpper = Upper.asSquare();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             // Bug fix: Missing transpose. VV, 31/Aug/2015.
             Tx[u[coeffI]] += mult(activeUpper[coeffI].T(), x[l[coeffI]]);
@@ -255,7 +255,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -264,7 +264,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -273,7 +273,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 // Use transpose upper coefficient
                 Tx[l[coeffI]] +=
@@ -289,7 +289,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
             }
@@ -298,7 +298,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
             }
@@ -307,7 +307,7 @@ void Foam::BlockLduMatrix<Type>::TmulCore
         {
             const squareTypeField& activeLower = Lower.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 // Bug fix: Missing transpose. VV, 31/Aug/2015.
                 Tx[l[coeffI]] += mult(activeLower[coeffI].T(), x[u[coeffI]]);
@@ -362,7 +362,7 @@ void Foam::BlockLduMatrix<Type>::segregateB
             linearType lt;
             squareType st;
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 contractLinear(lt, activeLower[coeffI]);
                 expandLinear(st, lt);
@@ -384,7 +384,7 @@ void Foam::BlockLduMatrix<Type>::segregateB
             linearType lt;
             squareType st;
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 contractLinear(lt, activeUpper[coeffI]);
                 expandLinear(st, lt);
@@ -396,7 +396,7 @@ void Foam::BlockLduMatrix<Type>::segregateB
             // is also needed
             if (symmetric())
             {
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     // Use transpose upper coefficient
                     contractLinear(lt, activeUpper[coeffI]);
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecouple.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecouple.C
index 7b8ecdff2..09e41c22a 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecouple.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecouple.C
@@ -57,7 +57,7 @@ void Foam::BlockLduMatrix<Type>::decoupledSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeUpper[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -72,7 +72,7 @@ void Foam::BlockLduMatrix<Type>::decoupledSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeUpper[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -100,7 +100,7 @@ void Foam::BlockLduMatrix<Type>::decoupledSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeLower[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -117,7 +117,7 @@ void Foam::BlockLduMatrix<Type>::decoupledSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeLower[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -160,7 +160,7 @@ void Foam::BlockLduMatrix<Type>::decoupledNegSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -175,7 +175,7 @@ void Foam::BlockLduMatrix<Type>::decoupledNegSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -203,7 +203,7 @@ void Foam::BlockLduMatrix<Type>::decoupledNegSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeLower[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -220,7 +220,7 @@ void Foam::BlockLduMatrix<Type>::decoupledNegSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeLower[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -263,7 +263,7 @@ void Foam::BlockLduMatrix<Type>::decoupledCheck() const
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiagCopy = DiagCopy.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
@@ -285,7 +285,7 @@ void Foam::BlockLduMatrix<Type>::decoupledCheck() const
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
@@ -317,7 +317,7 @@ void Foam::BlockLduMatrix<Type>::decoupledCheck() const
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiagCopy = DiagCopy.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiagCopy[l[coeffI]] += activeLower[coeffI];
                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
@@ -341,7 +341,7 @@ void Foam::BlockLduMatrix<Type>::decoupledCheck() const
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiagCopy[l[coeffI]] += activeLower[coeffI];
                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
@@ -415,7 +415,7 @@ void Foam::BlockLduMatrix<Type>::decoupledRelax
                 pTraits<typename TypeCoeffField::linearType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeUpper[coeffI]);
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -448,7 +448,7 @@ void Foam::BlockLduMatrix<Type>::decoupledRelax
                 pTraits<typename TypeCoeffField::scalarType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += mag(activeUpper[coeffI]);
                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -491,7 +491,7 @@ void Foam::BlockLduMatrix<Type>::decoupledRelax
                 pTraits<typename TypeCoeffField::linearType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeLower[coeffI]);
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -520,7 +520,7 @@ void Foam::BlockLduMatrix<Type>::decoupledRelax
                 pTraits<typename TypeCoeffField::scalarType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += mag(activeLower[coeffI]);
                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -559,7 +559,7 @@ void Foam::BlockLduMatrix<Type>::decoupledMultEqOp(const scalarField& sf)
         {
             scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeUpper[coeffI] *= sf[l[coeffI]];
             }
@@ -568,7 +568,7 @@ void Foam::BlockLduMatrix<Type>::decoupledMultEqOp(const scalarField& sf)
         {
             linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeUpper[coeffI] *= sf[l[coeffI]];
             }
@@ -585,7 +585,7 @@ void Foam::BlockLduMatrix<Type>::decoupledMultEqOp(const scalarField& sf)
         {
             scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 activeLower[coeffI] *= sf[u[coeffI]];
             }
@@ -594,7 +594,7 @@ void Foam::BlockLduMatrix<Type>::decoupledMultEqOp(const scalarField& sf)
         {
             linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 activeLower[coeffI] *= sf[u[coeffI]];
             }
@@ -636,7 +636,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -645,7 +645,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -659,7 +659,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
         {
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -668,7 +668,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
         {
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -682,7 +682,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
     {
         const scalarTypeField& activeUpper = Upper.asScalar();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -691,7 +691,7 @@ void Foam::BlockLduMatrix<Type>::decoupledAmulCore
     {
         const linearTypeField& activeUpper = Upper.asLinear();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -729,7 +729,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
     {
         const scalarTypeField& activeUpper = Upper.asScalar();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
         }
@@ -738,7 +738,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
     {
         const linearTypeField& activeUpper = Upper.asLinear();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
         }
@@ -752,7 +752,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -761,7 +761,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -775,7 +775,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
         {
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
             }
@@ -784,7 +784,7 @@ void Foam::BlockLduMatrix<Type>::decoupledTmulCore
         {
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
             }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecoupledHOps.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecoupledHOps.C
index c3c41e476..3cfec26aa 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecoupledHOps.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecoupledHOps.C
@@ -60,7 +60,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -69,7 +69,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
             }
@@ -83,7 +83,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
         {
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -92,7 +92,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
         {
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
             }
@@ -106,7 +106,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
     {
         const scalarTypeField& activeUpper = Upper.asScalar();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -115,7 +115,7 @@ Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
     {
         const linearTypeField& activeUpper = Upper.asLinear();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
         }
@@ -152,7 +152,7 @@ Foam::BlockLduMatrix<Type>::decoupledFaceH(const Field<Type>& x) const
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 // This can be optimised with a subtraction.
                 // Currently not done for clarity.  HJ, 31/Oct/2007
@@ -165,7 +165,7 @@ Foam::BlockLduMatrix<Type>::decoupledFaceH(const Field<Type>& x) const
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 // This can be optimised with a subtraction.
                 // Currently not done for clarity.  HJ, 31/Oct/2007
@@ -184,7 +184,7 @@ Foam::BlockLduMatrix<Type>::decoupledFaceH(const Field<Type>& x) const
             const scalarTypeField& activeUpper = Upper.asScalar();
             const scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[coeffI] =
                     mult(activeUpper[coeffI], x[u[coeffI]])
@@ -196,7 +196,7 @@ Foam::BlockLduMatrix<Type>::decoupledFaceH(const Field<Type>& x) const
             const linearTypeField& activeUpper = Upper.asLinear();
             const linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[coeffI] =
                     mult(activeUpper[coeffI], x[u[coeffI]])
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixHOps.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixHOps.C
index c5378e134..217eeb8a6 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixHOps.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixHOps.C
@@ -62,7 +62,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const scalarTypeField& activeUpper = Upper.asScalar();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
                 }
@@ -71,7 +71,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const linearTypeField& activeUpper = Upper.asLinear();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
                 }
@@ -80,7 +80,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const squareTypeField& activeUpper = Upper.asSquare();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     // Use transpose upper coefficient
                     result[u[coeffI]] -=
@@ -96,7 +96,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const scalarTypeField& activeLower = Lower.asScalar();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
                 }
@@ -105,7 +105,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const linearTypeField& activeLower = Lower.asLinear();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
                 }
@@ -114,7 +114,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
             {
                 const squareTypeField& activeLower = Lower.asSquare();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
                 }
@@ -128,7 +128,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
         {
             const scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -137,7 +137,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
         {
             const linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -146,7 +146,7 @@ Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
         {
             const squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
             }
@@ -187,7 +187,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
             {
                 const scalarTypeField& activeUpper = Upper.asScalar();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     // This can be optimised with a subtraction.
                     // Currently not done for clarity.  HJ, 31/Oct/2007
@@ -200,7 +200,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
             {
                 const linearTypeField& activeUpper = Upper.asLinear();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     // This can be optimised with a subtraction.
                     // Currently not done for clarity.  HJ, 31/Oct/2007
@@ -213,7 +213,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
             {
                 const squareTypeField& activeUpper = Upper.asSquare();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     // Use transpose upper coefficient
                     result[coeffI] =
@@ -231,7 +231,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
                 const scalarTypeField& activeUpper = Upper.asScalar();
                 const scalarTypeField& activeLower = Lower.asScalar();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[coeffI] =
                         mult(activeUpper[coeffI], x[u[coeffI]])
@@ -243,7 +243,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
                 const linearTypeField& activeUpper = Upper.asLinear();
                 const linearTypeField& activeLower = Lower.asLinear();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[coeffI] =
                         mult(activeUpper[coeffI], x[u[coeffI]])
@@ -255,7 +255,7 @@ Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
                 const squareTypeField& activeUpper = Upper.asSquare();
                 const squareTypeField& activeLower = Lower.asSquare();
 
-                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                for (label coeffI = 0; coeffI < u.size(); coeffI++)
                 {
                     result[coeffI] =
                         mult(activeUpper[coeffI], x[u[coeffI]])
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
index f9a193ec8..b7a7ccfa1 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
@@ -56,7 +56,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const squareTypeField& activeUpper = Upper.asSquare();
             squareTypeField& activeDiag = Diag.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeUpper[coeffI].T();
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -71,7 +71,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeUpper[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -86,7 +86,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeUpper[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -114,7 +114,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const squareTypeField& activeUpper = Upper.asSquare();
             squareTypeField& activeDiag = Diag.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeLower[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -131,7 +131,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeLower[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -148,7 +148,7 @@ void Foam::BlockLduMatrix<Type>::sumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] += activeLower[coeffI];
                 activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -192,7 +192,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const squareTypeField& activeUpper = Upper.asSquare();
             squareTypeField& activeDiag = Diag.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeUpper[coeffI].T();
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -207,7 +207,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -222,7 +222,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -250,7 +250,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const squareTypeField& activeUpper = Upper.asSquare();
             squareTypeField& activeDiag = Diag.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeLower[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -267,7 +267,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeDiag = Diag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeLower[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -284,7 +284,7 @@ void Foam::BlockLduMatrix<Type>::negSumDiag()
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeDiag = Diag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeDiag[l[coeffI]] -= activeLower[coeffI];
                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -362,7 +362,7 @@ void Foam::BlockLduMatrix<Type>::check() const
 
             const squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] +=
                     sumMagToDiag(activeUpper[coeffI].T());
@@ -395,7 +395,7 @@ void Foam::BlockLduMatrix<Type>::check() const
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
                 activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -422,7 +422,7 @@ void Foam::BlockLduMatrix<Type>::check() const
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
                 activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -481,7 +481,7 @@ void Foam::BlockLduMatrix<Type>::check() const
             const squareTypeField& activeLower = Lower.asSquare();
             const squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] +=
                     sumMagToDiag(activeLower[coeffI]);
@@ -516,7 +516,7 @@ void Foam::BlockLduMatrix<Type>::check() const
             const linearTypeField& activeUpper = Upper.asLinear();
             linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
                 activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -545,7 +545,7 @@ void Foam::BlockLduMatrix<Type>::check() const
             const scalarTypeField& activeUpper = Upper.asScalar();
             scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
                 activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -623,7 +623,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::squareType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeUpper[coeffI].T());
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -658,7 +658,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::linearType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeUpper[coeffI]);
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -691,7 +691,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::scalarType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += mag(activeUpper[coeffI]);
                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -741,7 +741,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::squareType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeLower[coeffI]);
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -778,7 +778,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::linearType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += cmptMag(activeLower[coeffI]);
                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
@@ -813,7 +813,7 @@ void Foam::BlockLduMatrix<Type>::relax
                 pTraits<typename TypeCoeffField::scalarType>::zero
             );
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 sumOff[u[coeffI]] += mag(activeLower[coeffI]);
                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -1047,7 +1047,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             scalarTypeField& activeUpper = Upper.asScalar();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeUpper[coeffI] *= sf[l[coeffI]];
             }
@@ -1056,7 +1056,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             linearTypeField& activeUpper = Upper.asLinear();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeUpper[coeffI] *= sf[l[coeffI]];
             }
@@ -1065,7 +1065,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             squareTypeField& activeUpper = Upper.asSquare();
 
-            for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+            for (label coeffI = 0; coeffI < l.size(); coeffI++)
             {
                 activeUpper[coeffI] *= sf[l[coeffI]];
             }
@@ -1082,7 +1082,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             scalarTypeField& activeLower = Lower.asScalar();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 activeLower[coeffI] *= sf[u[coeffI]];
             }
@@ -1091,7 +1091,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             linearTypeField& activeLower = Lower.asLinear();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 activeLower[coeffI] *= sf[u[coeffI]];
             }
@@ -1100,7 +1100,7 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
         {
             squareTypeField& activeLower = Lower.asSquare();
 
-            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            for (label coeffI = 0; coeffI < u.size(); coeffI++)
             {
                 activeLower[coeffI] *= sf[u[coeffI]];
             }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C
index 885a0a7c6..5f32b951f 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C
@@ -43,7 +43,7 @@ void BlockLduMatrix<scalar>::sumDiag()
     {
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiag[l[coeffI]] += activeUpper[coeffI];
             activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -54,7 +54,7 @@ void BlockLduMatrix<scalar>::sumDiag()
         const scalarField& activeLower = lower();
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiag[l[coeffI]] += activeLower[coeffI];
             activeDiag[u[coeffI]] += activeUpper[coeffI];
@@ -81,7 +81,7 @@ void BlockLduMatrix<scalar>::negSumDiag()
     {
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiag[l[coeffI]] -= activeUpper[coeffI];
             activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -92,7 +92,7 @@ void BlockLduMatrix<scalar>::negSumDiag()
         const scalarField& activeLower = lower();
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiag[l[coeffI]] -= activeLower[coeffI];
             activeDiag[u[coeffI]] -= activeUpper[coeffI];
@@ -120,7 +120,7 @@ void BlockLduMatrix<scalar>::check() const
     {
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiagCopy[l[coeffI]] -= activeUpper[coeffI];
             activeDiagCopy[u[coeffI]] -= activeUpper[coeffI];
@@ -138,7 +138,7 @@ void BlockLduMatrix<scalar>::check() const
         const scalarField& activeLower = lower();
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeDiagCopy[l[coeffI]] -= activeLower[coeffI];
             activeDiagCopy[u[coeffI]] -= activeUpper[coeffI];
@@ -180,7 +180,7 @@ void BlockLduMatrix<scalar>::relax
     {
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             sumOff[u[coeffI]] += mag(activeUpper[coeffI]);
             sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -191,7 +191,7 @@ void BlockLduMatrix<scalar>::relax
         const scalarField& activeLower = lower();
         const scalarField& activeUpper = upper();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             sumOff[u[coeffI]] += mag(activeLower[coeffI]);
             sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
@@ -220,7 +220,7 @@ void BlockLduMatrix<scalar>::operator*=(const scalarField& sf)
 
         const unallocLabelList& l = lduAddr().lowerAddr();
 
-        for (register label coeffI = 0; coeffI < l.size(); coeffI++)
+        for (label coeffI = 0; coeffI < l.size(); coeffI++)
         {
             activeUpper[coeffI] *= sf[l[coeffI]];
         }
@@ -232,7 +232,7 @@ void BlockLduMatrix<scalar>::operator*=(const scalarField& sf)
 
         const unallocLabelList& u = lduAddr().upperAddr();
 
-        for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+        for (label coeffI = 0; coeffI < u.size(); coeffI++)
         {
             activeLower[coeffI] *= sf[u[coeffI]];
         }
@@ -258,8 +258,8 @@ void BlockLduMatrix<scalar>::AmulCore
     {
         const scalar* const __restrict__ diagPtr = diag().begin();
 
-        register const label nCells = diag().size();
-        for (register label cell=0; cell<nCells; cell++)
+        const label nCells = diag().size();
+        for (label cell=0; cell<nCells; cell++)
         {
             // AmulCore must be additive to account for initialisation step
             // in ldu interfaces.  HJ, 6/Nov/2007
@@ -273,7 +273,7 @@ void BlockLduMatrix<scalar>::AmulCore
         {
             const scalar* const __restrict__ Upper = upper().begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 AX[U[coeffI]] += Upper[coeffI]*X[L[coeffI]];
                 AX[L[coeffI]] += Upper[coeffI]*X[U[coeffI]];
@@ -283,7 +283,7 @@ void BlockLduMatrix<scalar>::AmulCore
         {
             const scalar* const __restrict__ Lower = lower().begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 AX[U[coeffI]] += Lower[coeffI]*X[L[coeffI]];
                 AX[L[coeffI]] += Lower[coeffI]*X[U[coeffI]];
@@ -295,7 +295,7 @@ void BlockLduMatrix<scalar>::AmulCore
         const scalar* const __restrict__ Lower = lower().begin();
         const scalar* const __restrict__  Upper = upper().begin();
 
-        for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+        for (label coeffI = 0; coeffI < upper().size(); coeffI++)
         {
             AX[U[coeffI]] += Lower[coeffI]*X[L[coeffI]];
             AX[L[coeffI]] += Upper[coeffI]*X[U[coeffI]];
@@ -322,8 +322,8 @@ void BlockLduMatrix<scalar>::TmulCore
     {
         const scalar* const __restrict__ diagPtr = diag().begin();
 
-        register const label nCells = diag().size();
-        for (register label cell=0; cell<nCells; cell++)
+        const label nCells = diag().size();
+        for (label cell=0; cell<nCells; cell++)
         {
             // AmulCore must be additive to account for initialisation step
             // in ldu interfaces.  HJ, 6/Nov/2007
@@ -337,7 +337,7 @@ void BlockLduMatrix<scalar>::TmulCore
         {
             const scalar* const __restrict__ Upper = upper().begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 TX[U[coeffI]] += Upper[coeffI]*X[L[coeffI]];
                 TX[L[coeffI]] += Upper[coeffI]*X[U[coeffI]];
@@ -347,7 +347,7 @@ void BlockLduMatrix<scalar>::TmulCore
         {
             const scalar* const __restrict__ Lower = lower().begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 TX[U[coeffI]] += Lower[coeffI]*X[L[coeffI]];
                 TX[L[coeffI]] += Lower[coeffI]*X[U[coeffI]];
@@ -359,7 +359,7 @@ void BlockLduMatrix<scalar>::TmulCore
         const scalar* const __restrict__ Lower = lower().begin();
         const scalar* const __restrict__ Upper = upper().begin();
 
-        for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+        for (label coeffI = 0; coeffI < upper().size(); coeffI++)
         {
             TX[U[coeffI]] += Upper[coeffI]*X[L[coeffI]];
             TX[L[coeffI]] += Lower[coeffI]*X[U[coeffI]];
@@ -409,12 +409,12 @@ tmp<scalarField> BlockLduMatrix<scalar>::H(const scalarField& x) const
 
                 scalar* __restrict__ R = result.begin();
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]];
                 }
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[L[coeffI]] -= Upper[coeffI]*X[L[coeffI]];
                 }
@@ -425,12 +425,12 @@ tmp<scalarField> BlockLduMatrix<scalar>::H(const scalarField& x) const
 
                 scalar* __restrict__ R = result.begin();
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[U[coeffI]] -= Lower[coeffI]*X[U[coeffI]];
                 }
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]];
                 }
@@ -443,12 +443,12 @@ tmp<scalarField> BlockLduMatrix<scalar>::H(const scalarField& x) const
 
             scalar* __restrict__ R = result.begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]];
             }
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
                 R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]];
             }
@@ -480,7 +480,7 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
 
                 scalar* __restrict__ R = result.begin();
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[coeffI] = Upper[coeffI]*(X[U[coeffI]] - X[L[coeffI]]);
                 }
@@ -496,7 +496,7 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
 
                 scalar* __restrict__ R = result.begin();
 
-                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                for (label coeffI = 0; coeffI < upper().size(); coeffI++)
                 {
                     R[coeffI] = Lower[coeffI]*(X[U[coeffI]] - X[L[coeffI]]);
                 }
@@ -514,9 +514,11 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
 
             scalar* __restrict__ R = result.begin();
 
-            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            for (label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
-                R[coeffI] = Upper[coeffI]*X[U[coeffI]] - Lower[coeffI]*X[L[coeffI]];
+                R[coeffI] =
+                    Upper[coeffI]*X[U[coeffI]]
+                  - Lower[coeffI]*X[L[coeffI]];
             }
         }
     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C
index 89f53d002..1a021e5a6 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C
@@ -826,7 +826,7 @@ void Foam::BlockCholeskyPrecon<Type>::ILUmultiply
     const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr();
     const unallocLabelList& losortAddr = this->matrix_.lduAddr().losortAddr();
 
-    register label losortCoeff;
+    label losortCoeff;
 
     forAll (lower, coeffI)
     {
@@ -875,7 +875,7 @@ void Foam::BlockCholeskyPrecon<Type>::ILUmultiplyTranspose
     const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr();
     const unallocLabelList& losortAddr = this->matrix_.lduAddr().losortAddr();
 
-    register label losortCoeff;
+    label losortCoeff;
 
     //HJ Not sure if the coefficient itself needs to be transposed.
     // HJ, 30/Oct/2007
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
index 391cb1535..1d9fd649f 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
@@ -148,10 +148,10 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
         true             // switch to lhs of system
     );
 
-    register label fStart, fEnd, curCoeff;
+    label fStart, fEnd, curCoeff;
 
     // Forward sweep
-    for (register label rowI = 0; rowI < nRows; rowI++)
+    for (label rowI = 0; rowI < nRows; rowI++)
     {
         Type& curX = x[rowI];
 
@@ -181,7 +181,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
     }
 
     // Reverse sweep
-    for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+    for (label rowI = nRows - 1; rowI >= 0; rowI--)
     {
         Type& curX = x[rowI];
 
@@ -249,10 +249,10 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
         true             // switch to lhs of system
     );
 
-    register label fStart, fEnd, curCoeff;
+    label fStart, fEnd, curCoeff;
 
     // Forward sweep
-    for (register label rowI = 0; rowI < nRows; rowI++)
+    for (label rowI = 0; rowI < nRows; rowI++)
     {
         Type& curX = x[rowI];
 
@@ -280,7 +280,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
     }
 
     // Reverse sweep
-    for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+    for (label rowI = nRows - 1; rowI >= 0; rowI--)
     {
         Type& curX = x[rowI];
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C
index e1cc04e8b..e1895c9b4 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C
@@ -88,12 +88,12 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
 
         // Define start and end face ("virtual" face when extended addressing
         // is used) of this row/column.
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -103,7 +103,7 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
             zDiag = diagPtr[rowI];
 
             // Initialize temporary working row and column fields
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper)
@@ -118,13 +118,13 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
             // Lower/upper coeff loop (i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
@@ -137,14 +137,14 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
                 );
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortCoeff + 1 = start)
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -168,10 +168,10 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
             diagRowI = mult.inverse(zDiag);
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -192,23 +192,23 @@ void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -369,8 +369,8 @@ void Foam::BlockILUC0Precon<Type>::LUSubstitute
     // Initialize x field
     x = b;
 
-    register label losortCoeffI;
-    register label rowI;
+    label losortCoeffI;
+    label rowI;
 
     // Forward substitution loop
     forAll (lower, coeffI)
@@ -451,8 +451,8 @@ void Foam::BlockILUC0Precon<Type>::LUSubstituteT
         );
     }
 
-    register label losortCoeffI;
-    register label rowI;
+    label losortCoeffI;
+    label rowI;
 
     // Forward substitution loop
     forAll (upper, coeffI)
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C
index cfbfdce59..e921d1c67 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C
@@ -90,12 +90,12 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
 
         // Define start and end face ("virtual" face when extended addressing is
         // used) of this row/column.
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -105,7 +105,7 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
             zDiag = diagPtr[rowI];
 
             // Initialize temporary working row field
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper), w addressed by neighbour of face (row index for
@@ -121,13 +121,13 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
             // Lower/upper coeff loop (i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
@@ -136,14 +136,14 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
                 zDiag -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortCoeff + 1 = start)
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -158,10 +158,10 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
             diagRowI = 1.0/zDiag;
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -178,23 +178,23 @@ void BlockILUC0Precon<scalar>::calcActiveTypeFactorization
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
index 5177a122b..d38012486 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
@@ -91,12 +91,12 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
 
         // Define start and end face ("virtual" face when extended addressing
         // is used) of this row/column.
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -106,7 +106,7 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
             zDiag = diagPtr[rowI];
 
             // Initialize temporary working row and column fields
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper)
@@ -121,13 +121,13 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
             // Lower/upper coeff loop (i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
@@ -140,14 +140,14 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
                 );
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortCoeff + 1 = start)
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -171,10 +171,10 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
             diagRowI = mult.inverse(zDiag);
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -195,23 +195,23 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -375,8 +375,8 @@ void Foam::BlockILUCpPrecon<Type>::LUSubstitute
     // Initialize x field
     x = b;
 
-    register label losortCoeffI;
-    register label rowI;
+    label losortCoeffI;
+    label rowI;
 
     // Forward substitution loop
     forAll (lower, coeffI)
@@ -457,8 +457,8 @@ void Foam::BlockILUCpPrecon<Type>::LUSubstituteT
         );
     }
 
-    register label losortCoeffI;
-    register label rowI;
+    label losortCoeffI;
+    label rowI;
 
     // Forward substitution loop
     forAll (upper, coeffI)
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C
index f8dd938b2..ee812b4a2 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C
@@ -90,12 +90,12 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
 
         // Define start and end face ("virtual" face when extended addressing is
         // used) of this row/column.
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -105,7 +105,7 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
             zDiag = diagPtr[rowI];
 
             // Initialize temporary working row field
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper), w addressed by neighbour of face (row index for
@@ -121,13 +121,13 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
             // Lower/upper coeff loop (i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
@@ -136,14 +136,14 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
                 zDiag -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortCoeff + 1 = start)
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -158,10 +158,10 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
             diagRowI = 1.0/zDiag;
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -178,23 +178,23 @@ void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
index 8557e0498..36a62ae5f 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
+++ b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
@@ -237,7 +237,7 @@ void Foam::extendedLduAddressing::calcFaceMap() const
         // Loop through extended neighbouring faces (upper matrix coeffs)
         for
         (
-            register label extFaceI = startLabel;
+            label extFaceI = startLabel;
             extFaceI < endLabel;
             ++extFaceI
         )
diff --git a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
index a01767461..60670c230 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
+++ b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
@@ -86,8 +86,8 @@ void Foam::lduMatrix::AmulCore
     {
         const scalar* const __restrict__ diagPtr = diag().begin();
 
-        register const label nCells = diag().size();
-        for (register label cell=0; cell<nCells; cell++)
+        const label nCells = diag().size();
+        for (label cell=0; cell<nCells; cell++)
         {
             // AmulCore must be additive to account for initialisation step
             // in ldu interfaces.  HJ, 6/Nov/2007
@@ -105,9 +105,9 @@ void Foam::lduMatrix::AmulCore
         const scalar* const __restrict__ upperPtr = upper().begin();
         const scalar* const __restrict__ lowerPtr = lower().begin();
 
-        register const label nFaces = upper().size();
+        const label nFaces = upper().size();
 
-        for (register label face=0; face<nFaces; face++)
+        for (label face=0; face<nFaces; face++)
         {
             AxPtr[uPtr[face]] += lowerPtr[face]*xPtr[lPtr[face]];
             AxPtr[lPtr[face]] += upperPtr[face]*xPtr[uPtr[face]];
@@ -171,8 +171,8 @@ void Foam::lduMatrix::TmulCore
     {
         const scalar* const __restrict__ diagPtr = diag().begin();
 
-        register const label nCells = diag().size();
-        for (register label cell=0; cell<nCells; cell++)
+        const label nCells = diag().size();
+        for (label cell=0; cell<nCells; cell++)
         {
             // TmulCore must be additive to account for initialisation step
             // in ldu interfaces.  HJ, 6/Nov/2007
@@ -190,8 +190,8 @@ void Foam::lduMatrix::TmulCore
         const scalar* const __restrict__ lowerPtr = lower().begin();
         const scalar* const __restrict__ upperPtr = upper().begin();
 
-        register const label nFaces = upper().size();
-        for (register label face=0; face<nFaces; face++)
+        const label nFaces = upper().size();
+        for (label face=0; face<nFaces; face++)
         {
             TxPtr[uPtr[face]] += upperPtr[face]*xPtr[lPtr[face]];
             TxPtr[lPtr[face]] += lowerPtr[face]*xPtr[uPtr[face]];
@@ -217,15 +217,15 @@ void Foam::lduMatrix::sumA
     const scalar* __restrict__ lowerPtr = lower().begin();
     const scalar* __restrict__ upperPtr = upper().begin();
 
-    register const label nCells = diag().size();
-    register const label nFaces = upper().size();
+    const label nCells = diag().size();
+    const label nFaces = upper().size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         sumAPtr[cell] = diagPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         sumAPtr[uPtr[face]] += lowerPtr[face];
         sumAPtr[lPtr[face]] += upperPtr[face];
diff --git a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
index 66dcf17a8..1099020d3 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
+++ b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
@@ -39,7 +39,7 @@ void Foam::lduMatrix::sumDiag()
     const unallocLabelList& l = lduAddr().lowerAddr();
     const unallocLabelList& u = lduAddr().upperAddr();
 
-    for (register label odcI = 0; odcI < l.size(); odcI++)
+    for (label odcI = 0; odcI < l.size(); odcI++)
     {
         Diag[l[odcI]] += Lower[odcI];
         Diag[u[odcI]] += Upper[odcI];
@@ -56,7 +56,7 @@ void Foam::lduMatrix::negSumDiag()
     const unallocLabelList& l = lduAddr().lowerAddr();
     const unallocLabelList& u = lduAddr().upperAddr();
 
-    for (register label odcI = 0; odcI < l.size(); odcI++)
+    for (label odcI = 0; odcI < l.size(); odcI++)
     {
         Diag[l[odcI]] -= Lower[odcI];
         Diag[u[odcI]] -= Upper[odcI];
@@ -75,7 +75,7 @@ void Foam::lduMatrix::sumMagOffDiag
     const unallocLabelList& l = lduAddr().lowerAddr();
     const unallocLabelList& u = lduAddr().upperAddr();
 
-    for (register label odcI = 0; odcI < l.size(); odcI++)
+    for (label odcI = 0; odcI < l.size(); odcI++)
     {
         sumOff[u[odcI]] += mag(Lower[odcI]);
         sumOff[l[odcI]] += mag(Upper[odcI]);
@@ -302,7 +302,7 @@ void Foam::lduMatrix::operator*=(const scalarField& sf)
 
         const unallocLabelList& l = lduAddr().lowerAddr();
 
-        for (register label odcI = 0; odcI < upper.size(); odcI++)
+        for (label odcI = 0; odcI < upper.size(); odcI++)
         {
             upper[odcI] *= sf[l[odcI]];
         }
@@ -314,7 +314,7 @@ void Foam::lduMatrix::operator*=(const scalarField& sf)
 
         const unallocLabelList& u = lduAddr().upperAddr();
 
-        for (register label odcI = 0; odcI < lower.size(); odcI++)
+        for (label odcI = 0; odcI < lower.size(); odcI++)
         {
             lower[odcI] *= sf[u[odcI]];
         }
@@ -360,9 +360,9 @@ Foam::tmp<Foam::scalarField> Foam::lduMatrix::H1() const
         const scalar* __restrict__ lowerPtr = lower().begin();
         const scalar* __restrict__ upperPtr = upper().begin();
 
-        register const label nOdcIs = upper().size();
+        const label nOdcIs = upper().size();
 
-        for (register label odcI = 0; odcI < nOdcIs; odcI++)
+        for (label odcI = 0; odcI < nOdcIs; odcI++)
         {
             H1Ptr[uPtr[odcI]] -= lowerPtr[odcI];
             H1Ptr[lPtr[odcI]] -= upperPtr[odcI];
diff --git a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
index d363dbcf7..e034f039d 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
+++ b/src/foam/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
@@ -52,9 +52,9 @@ Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const
         const scalar* __restrict__ lowerPtr = lower().begin();
         const scalar* __restrict__ upperPtr = upper().begin();
 
-        register const label nFaces = upper().size();
+        const label nFaces = upper().size();
 
-        for (register label face=0; face<nFaces; face++)
+        for (label face=0; face<nFaces; face++)
         {
             HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
             HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
@@ -92,7 +92,7 @@ Foam::lduMatrix::faceH(const Field<Type>& psi) const
         const unallocLabelList& l = lduAddr().lowerAddr();
         const unallocLabelList& u = lduAddr().upperAddr();
 
-        for (register label face=0; face<l.size(); face++)
+        for (label face=0; face<l.size(); face++)
         {
             faceHpsi[face] = Upper[face]*psi[u[face]]
                 - Lower[face]*psi[l[face]];
diff --git a/src/foam/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C b/src/foam/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C
index 8afa89b94..350cbb2d5 100644
--- a/src/foam/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C
+++ b/src/foam/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C
@@ -52,17 +52,17 @@ void Foam::DICPreconditioner::calcReciprocalD
     const scalar* const __restrict__ upperPtr = matrix.upper().begin();
 
     // Calculate the DIC diagonal
-    register const label nFaces = matrix.upper().size();
-    for (register label face=0; face<nFaces; face++)
+    const label nFaces = matrix.upper().size();
+    for (label face=0; face<nFaces; face++)
     {
         rDPtr[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rDPtr[lPtr[face]];
     }
 
 
     // Calculate the reciprocal of the preconditioned diagonal
-    register const label nCells = rD.size();
+    const label nCells = rD.size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         rDPtr[cell] = 1.0/rDPtr[cell];
     }
@@ -112,21 +112,21 @@ void Foam::DICPreconditioner::precondition
         matrix_.lduAddr().lowerAddr().begin();
     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 
-    register label nCells = wA.size();
-    register label nFaces = matrix_.upper().size();
-    register label nFacesM1 = nFaces - 1;
+    label nCells = wA.size();
+    label nFaces = matrix_.upper().size();
+    label nFacesM1 = nFaces - 1;
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
     }
 
-    for (register label face=nFacesM1; face>=0; face--)
+    for (label face=nFacesM1; face>=0; face--)
     {
         wAPtr[lPtr[face]] -= rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
     }
diff --git a/src/foam/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C b/src/foam/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C
index 2ac1277ac..80228c084 100644
--- a/src/foam/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C
+++ b/src/foam/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C
@@ -83,17 +83,17 @@ void Foam::DILUPreconditioner::calcReciprocalD
     const scalar* const __restrict__ upperPtr = matrix.upper().begin();
     const scalar* const __restrict__ lowerPtr = matrix.lower().begin();
 
-    register label nFaces = matrix.upper().size();
-    for (register label face=0; face<nFaces; face++)
+    label nFaces = matrix.upper().size();
+    for (label face=0; face<nFaces; face++)
     {
         rDPtr[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rDPtr[lPtr[face]];
     }
 
 
     // Calculate the reciprocal of the preconditioned diagonal
-    register label nCells = rD.size();
+    label nCells = rD.size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         rDPtr[cell] = 1.0/rDPtr[cell];
     }
@@ -121,26 +121,26 @@ void Foam::DILUPreconditioner::precondition
     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
     const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 
-    register label nCells = wA.size();
-    register label nFaces = matrix_.upper().size();
-    register label nFacesM1 = nFaces - 1;
+    label nCells = wA.size();
+    label nFaces = matrix_.upper().size();
+    label nFacesM1 = nFaces - 1;
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
     }
 
 
-    register label sface;
+    label sface;
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         sface = losortPtr[face];
         wAPtr[uPtr[sface]] -=
             rDPtr[uPtr[sface]]*lowerPtr[sface]*wAPtr[lPtr[sface]];
     }
 
-    for (register label face=nFacesM1; face>=0; face--)
+    for (label face=nFacesM1; face>=0; face--)
     {
         wAPtr[lPtr[face]] -=
             rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
@@ -169,25 +169,25 @@ void Foam::DILUPreconditioner::preconditionT
     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
     const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
 
-    register label nCells = wT.size();
-    register label nFaces = matrix_.upper().size();
-    register label nFacesM1 = nFaces - 1;
+    label nCells = wT.size();
+    label nFaces = matrix_.upper().size();
+    label nFacesM1 = nFaces - 1;
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wTPtr[cell] = rDPtr[cell]*rTPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         wTPtr[uPtr[face]] -=
             rDPtr[uPtr[face]]*upperPtr[face]*wTPtr[lPtr[face]];
     }
 
 
-    register label sface;
+    label sface;
 
-    for (register label face=nFacesM1; face>=0; face--)
+    for (label face=nFacesM1; face>=0; face--)
     {
         sface = losortPtr[face];
         wTPtr[lPtr[sface]] -=
diff --git a/src/foam/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.C b/src/foam/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.C
index 30aea1653..3e83232e0 100644
--- a/src/foam/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.C
+++ b/src/foam/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.C
@@ -69,21 +69,21 @@ Foam::FDICPreconditioner::FDICPreconditioner
         matrix_.lduAddr().lowerAddr().begin();
     const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
 
-    register label nCells = rD_.size();
-    register label nFaces = matrix_.upper().size();
+    label nCells = rD_.size();
+    label nFaces = matrix_.upper().size();
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         rDPtr[uPtr[face]] -= sqr(upperPtr[face])/rDPtr[lPtr[face]];
     }
 
     // Generate reciprocal FDIC
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         rDPtr[cell] = 1.0/rDPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         rDuUpperPtr[face] = rDPtr[uPtr[face]]*upperPtr[face];
         rDlUpperPtr[face] = rDPtr[lPtr[face]]*upperPtr[face];
@@ -112,21 +112,21 @@ void Foam::FDICPreconditioner::precondition
     const scalar* const __restrict__ rDuUpperPtr = rDuUpper_.begin();
     const scalar* const __restrict__ rDlUpperPtr = rDlUpper_.begin();
 
-    register label nCells = wA.size();
-    register label nFaces = matrix_.upper().size();
-    register label nFacesM1 = nFaces - 1;
+    label nCells = wA.size();
+    label nFaces = matrix_.upper().size();
+    label nFacesM1 = nFaces - 1;
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
     }
 
-    for (register label face=0; face<nFaces; face++)
+    for (label face=0; face<nFaces; face++)
     {
         wAPtr[uPtr[face]] -= rDuUpperPtr[face]*wAPtr[lPtr[face]];
     }
 
-    for (register label face=nFacesM1; face>=0; face--)
+    for (label face=nFacesM1; face>=0; face--)
     {
         wAPtr[lPtr[face]] -= rDlUpperPtr[face]*wAPtr[uPtr[face]];
     }
diff --git a/src/foam/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.C b/src/foam/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.C
index 8b2914df5..0ae60d5aa 100644
--- a/src/foam/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.C
+++ b/src/foam/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.C
@@ -64,10 +64,10 @@ Foam::diagonalPreconditioner::diagonalPreconditioner
     scalar* __restrict__ rDPtr = rD.begin();
     const scalar* __restrict__ DPtr = matrix_.diag().begin();
 
-    register label nCells = rD.size();
+    label nCells = rD.size();
 
     // Generate reciprocal diagonal
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         rDPtr[cell] = 1.0/DPtr[cell];
     }
@@ -87,9 +87,9 @@ void Foam::diagonalPreconditioner::precondition
     const scalar* __restrict__ rAPtr = rA.begin();
     const scalar* __restrict__ rDPtr = rD.begin();
 
-    register label nCells = wA.size();
+    label nCells = wA.size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
     }
diff --git a/src/foam/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.C b/src/foam/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.C
index cfae0eaf9..4d5876a43 100644
--- a/src/foam/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.C
+++ b/src/foam/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.C
@@ -74,9 +74,9 @@ void Foam::noPreconditioner::precondition
     scalar* __restrict__ wAPtr = wA.begin();
     const scalar* __restrict__ rAPtr = rA.begin();
 
-    register label nCells = wA.size();
+    label nCells = wA.size();
 
-    for (register label cell=0; cell<nCells; cell++)
+    for (label cell=0; cell<nCells; cell++)
     {
         wAPtr[cell] = rAPtr[cell];
     }
diff --git a/src/foam/matrices/lduMatrix/smoothers/DIC/DICSmoother.C b/src/foam/matrices/lduMatrix/smoothers/DIC/DICSmoother.C
index 8361c3dd5..c1c7490f2 100644
--- a/src/foam/matrices/lduMatrix/smoothers/DIC/DICSmoother.C
+++ b/src/foam/matrices/lduMatrix/smoothers/DIC/DICSmoother.C
@@ -96,17 +96,17 @@ void Foam::DICSmoother::smooth
 
         rA *= rD_;
 
-        register label nFaces = matrix_.upper().size();
-        for (register label face=0; face<nFaces; face++)
+        label nFaces = matrix_.upper().size();
+        for (label face=0; face<nFaces; face++)
         {
-            register label u = uPtr[face];
+            label u = uPtr[face];
             rAPtr[u] -= rDPtr[u]*upperPtr[face]*rAPtr[lPtr[face]];
         }
 
-        register label nFacesM1 = nFaces - 1;
-        for (register label face=nFacesM1; face>=0; face--)
+        label nFacesM1 = nFaces - 1;
+        for (label face=nFacesM1; face>=0; face--)
         {
-            register label l = lPtr[face];
+            label l = lPtr[face];
             rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
         }
 
diff --git a/src/foam/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C b/src/foam/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C
index 9132733fc..725b05002 100644
--- a/src/foam/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C
+++ b/src/foam/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C
@@ -98,17 +98,17 @@ void Foam::DILUSmoother::smooth
 
         rA *= rD_;
 
-        register label nFaces = matrix_.upper().size();
-        for (register label face=0; face<nFaces; face++)
+        label nFaces = matrix_.upper().size();
+        for (label face=0; face<nFaces; face++)
         {
-            register label u = uPtr[face];
+            label u = uPtr[face];
             rAPtr[u] -= rDPtr[u]*lowerPtr[face]*rAPtr[lPtr[face]];
         }
 
-        register label nFacesM1 = nFaces - 1;
-        for (register label face=nFacesM1; face>=0; face--)
+        label nFacesM1 = nFaces - 1;
+        for (label face=nFacesM1; face>=0; face--)
         {
-            register label l = lPtr[face];
+            label l = lPtr[face];
             rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
         }
 
diff --git a/src/foam/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C b/src/foam/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
index 11c9f6a2d..55696e73e 100644
--- a/src/foam/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
+++ b/src/foam/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
@@ -72,23 +72,23 @@ void Foam::GaussSeidelSmoother::smooth
     const label nSweeps
 )
 {
-    register scalar* __restrict__ xPtr = x.begin();
+    scalar* __restrict__ xPtr = x.begin();
 
-    register const label nCells = x.size();
+    const label nCells = x.size();
 
     scalarField bPrime(nCells);
-    register scalar* __restrict__ bPrimePtr = bPrime.begin();
+    scalar* __restrict__ bPrimePtr = bPrime.begin();
 
-    register const scalar* const __restrict__ diagPtr = matrix.diag().begin();
-    register const scalar* const __restrict__ upperPtr =
+    const scalar* const __restrict__ diagPtr = matrix.diag().begin();
+    const scalar* const __restrict__ upperPtr =
         matrix.upper().begin();
-    register const scalar* const __restrict__ lowerPtr =
+    const scalar* const __restrict__ lowerPtr =
         matrix.lower().begin();
 
-    register const label* const __restrict__ uPtr =
+    const label* const __restrict__ uPtr =
         matrix.lduAddr().upperAddr().begin();
 
-    register const label* const __restrict__ ownStartPtr =
+    const label* const __restrict__ ownStartPtr =
         matrix.lduAddr().ownerStartAddr().begin();
 
 
@@ -133,11 +133,11 @@ void Foam::GaussSeidelSmoother::smooth
             true         // switch to lhs
         );
 
-        register scalar curX;
-        register label fStart;
-        register label fEnd = ownStartPtr[0];
+        scalar curX;
+        label fStart;
+        label fEnd = ownStartPtr[0];
 
-        for (register label cellI = 0; cellI < nCells; cellI++)
+        for (label cellI = 0; cellI < nCells; cellI++)
         {
             // Start and end of this row
             fStart = fEnd;
@@ -147,7 +147,7 @@ void Foam::GaussSeidelSmoother::smooth
             curX = bPrimePtr[cellI];
 
             // Accumulate the owner product side
-            for (register label curFace = fStart; curFace < fEnd; curFace++)
+            for (label curFace = fStart; curFace < fEnd; curFace++)
             {
                 curX -= upperPtr[curFace]*xPtr[uPtr[curFace]];
             }
@@ -156,7 +156,7 @@ void Foam::GaussSeidelSmoother::smooth
             curX /= diagPtr[cellI];
 
             // Distribute the neighbour side using current x
-            for (register label curFace = fStart; curFace < fEnd; curFace++)
+            for (label curFace = fStart; curFace < fEnd; curFace++)
             {
                 bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curX;
             }
diff --git a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C
index e49555890..687dbcb82 100644
--- a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C
+++ b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C
@@ -341,7 +341,7 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
 
                 for
                 (
-                    register label pfI = patchRowStart[patchFace];
+                    label pfI = patchRowStart[patchFace];
                     pfI < patchRowStart[patchFace + 1];
                     pfI++
                 )
@@ -1484,7 +1484,7 @@ void Foam::ggiSAMGInterface::expandCrMatrixToZone(crMatrix& patchP) const
         List<labelField> cols(patchCrAddr.nRows());
         List<scalarField> coeffs(patchCrAddr.nRows());
 
-        for (register label faceI = 0; faceI < patchCrAddr.nRows(); faceI++)
+        for (label faceI = 0; faceI < patchCrAddr.nRows(); faceI++)
         {
             // Unpack row
             const label rowStart = patchRowStart[faceI];
@@ -1496,7 +1496,7 @@ void Foam::ggiSAMGInterface::expandCrMatrixToZone(crMatrix& patchP) const
             coeffs[faceI].setSize(rowLength);
             scalarField& curCoeffs = coeffs[faceI];
 
-            for (register label coeffI = 0; coeffI < rowLength; coeffI++)
+            for (label coeffI = 0; coeffI < rowLength; coeffI++)
             {
                 curCols[coeffI] = patchCol[rowStart + coeffI];
                 curCoeffs[coeffI] = patchCoeff[rowStart + coeffI];
diff --git a/src/foam/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/foam/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
index a1cf6ae5a..7825a9996 100644
--- a/src/foam/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
+++ b/src/foam/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
@@ -72,7 +72,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
     // --- Setup class containing solver performance data
     lduSolverPerformance solverPerf(typeName, fieldName());
 
-    register label nCells = x.size();
+    label nCells = x.size();
 
     scalar* __restrict__ xPtr = x.begin();
 
@@ -147,7 +147,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
 
             if (solverPerf.nIterations() == 0)
             {
-                for (register label cell=0; cell<nCells; cell++)
+                for (label cell=0; cell<nCells; cell++)
                 {
                     pAPtr[cell] = wAPtr[cell];
                     pTPtr[cell] = wTPtr[cell];
@@ -157,7 +157,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
             {
                 scalar beta = wArT/wArTold;
 
-                for (register label cell=0; cell<nCells; cell++)
+                for (label cell=0; cell<nCells; cell++)
                 {
                     pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
                     pTPtr[cell] = wTPtr[cell] + beta*pTPtr[cell];
@@ -180,7 +180,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
 
             scalar alpha = wArT/wApT;
 
-            for (register label cell=0; cell<nCells; cell++)
+            for (label cell=0; cell<nCells; cell++)
             {
                 xPtr[cell] += alpha*pAPtr[cell];
                 rAPtr[cell] -= alpha*wAPtr[cell];
diff --git a/src/foam/matrices/lduMatrix/solvers/PCG/PCG.C b/src/foam/matrices/lduMatrix/solvers/PCG/PCG.C
index 059babef2..46161fe35 100644
--- a/src/foam/matrices/lduMatrix/solvers/PCG/PCG.C
+++ b/src/foam/matrices/lduMatrix/solvers/PCG/PCG.C
@@ -72,7 +72,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
     // --- Setup class containing solver performance data
     lduSolverPerformance solverPerf(typeName, fieldName());
 
-    register label nCells = x.size();
+    label nCells = x.size();
 
     scalar* __restrict__ xPtr = x.begin();
 
@@ -137,7 +137,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
 
             if (solverPerf.nIterations() == 0)
             {
-                for (register label cell=0; cell<nCells; cell++)
+                for (label cell=0; cell<nCells; cell++)
                 {
                     pAPtr[cell] = wAPtr[cell];
                 }
@@ -146,7 +146,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
             {
                 scalar beta = wArA/wArAold;
 
-                for (register label cell=0; cell<nCells; cell++)
+                for (label cell=0; cell<nCells; cell++)
                 {
                     pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
                 }
@@ -167,7 +167,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
 
             scalar alpha = wArA/wApA;
 
-            for (register label cell=0; cell<nCells; cell++)
+            for (label cell=0; cell<nCells; cell++)
             {
                 xPtr[cell] += alpha*pAPtr[cell];
                 rAPtr[cell] -= alpha*wAPtr[cell];
diff --git a/src/foam/matrices/scalarMatrices/scalarMatrices.C b/src/foam/matrices/scalarMatrices/scalarMatrices.C
index 9a28e99a7..fcbfbe8f8 100644
--- a/src/foam/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/foam/matrices/scalarMatrices/scalarMatrices.C
@@ -50,11 +50,11 @@ void Foam::multiply
 
     ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
 
-    for(register label i = 0; i < A.n(); i++)
+    for(label i = 0; i < A.n(); i++)
     {
-        for(register label j = 0; j < B.m(); j++)
+        for(label j = 0; j < B.m(); j++)
         {
-            for(register label l = 0; l < B.n(); l++)
+            for(label l = 0; l < B.n(); l++)
             {
                 ans[i][j] += A[i][l]*B[l][j];
             }
@@ -101,14 +101,14 @@ void Foam::multiply
 
     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 
-    for(register label i = 0; i < A.n(); i++)
+    for(label i = 0; i < A.n(); i++)
     {
-        for(register label g = 0; g < C.m(); g++)
+        for(label g = 0; g < C.m(); g++)
         {
-            for(register label l = 0; l < C.n(); l++)
+            for(label l = 0; l < C.n(); l++)
             {
                 scalar ab = 0;
-                for(register label j = 0; j < A.m(); j++)
+                for(label j = 0; j < A.m(); j++)
                 {
                     ab += A[i][j]*B[j][l];
                 }
@@ -157,11 +157,11 @@ void Foam::multiply
 
     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 
-    for(register label i = 0; i < A.n(); i++)
+    for(label i = 0; i < A.n(); i++)
     {
-        for(register label g = 0; g < C.m(); g++)
+        for(label g = 0; g < C.m(); g++)
         {
-            for(register label l = 0; l < C.n(); l++)
+            for(label l = 0; l < C.n(); l++)
             {
                 ans[i][g] += C[l][g] * A[i][l]*B[l];
             }
diff --git a/src/foam/matrices/scalarMatrices/scalarSquareMatrix.C b/src/foam/matrices/scalarMatrices/scalarSquareMatrix.C
index f2ab6f806..249bddea4 100644
--- a/src/foam/matrices/scalarMatrices/scalarSquareMatrix.C
+++ b/src/foam/matrices/scalarMatrices/scalarSquareMatrix.C
@@ -76,13 +76,13 @@ void Foam::scalarSquareMatrix::LUDecompose
     label n = matrix.n();
     scalar vv[n];
 
-    for (register label i=0; i<n; i++)
+    for (label i=0; i<n; i++)
     {
         scalar largestCoeff = 0.0;
         scalar temp;
         const scalar* __restrict__ matrixi = matrix[i];
 
-        for (register label j=0; j<n; j++)
+        for (label j=0; j<n; j++)
         {
             if ((temp = mag(matrixi[j])) > largestCoeff)
             {
@@ -102,16 +102,16 @@ void Foam::scalarSquareMatrix::LUDecompose
         vv[i] = 1.0/largestCoeff;
     }
 
-    for (register label j=0; j<n; j++)
+    for (label j=0; j<n; j++)
     {
         scalar* __restrict__ matrixj = matrix[j];
 
-        for (register label i=0; i<j; i++)
+        for (label i=0; i<j; i++)
         {
             scalar* __restrict__ matrixi = matrix[i];
 
             scalar sum = matrixi[j];
-            for (register label k=0; k<i; k++)
+            for (label k=0; k<i; k++)
             {
                 sum -= matrixi[k]*matrix[k][j];
             }
@@ -121,12 +121,12 @@ void Foam::scalarSquareMatrix::LUDecompose
         label iMax = 0;
 
         scalar largestCoeff = 0.0;
-        for (register label i=j; i<n; i++)
+        for (label i=j; i<n; i++)
         {
             scalar* __restrict__ matrixi = matrix[i];
             scalar sum = matrixi[j];
 
-            for (register label k=0; k<j; k++)
+            for (label k=0; k<j; k++)
             {
                 sum -= matrixi[k]*matrix[k][j];
             }
@@ -147,7 +147,7 @@ void Foam::scalarSquareMatrix::LUDecompose
         {
             scalar* __restrict__ matrixiMax = matrix[iMax];
 
-            for (register label k=0; k<n; k++)
+            for (label k = 0; k < n; k++)
             {
                 Swap(matrixj[k], matrixiMax[k]);
             }
@@ -164,7 +164,7 @@ void Foam::scalarSquareMatrix::LUDecompose
         {
             scalar rDiag = 1.0/matrixj[j];
 
-            for (register label i=j+1; i<n; i++)
+            for (label i = j + 1; i < n; i++)
             {
                 matrix[i][j] *= rDiag;
             }
diff --git a/src/foam/matrices/scalarMatrices/scalarSquareMatrixTemplates.C b/src/foam/matrices/scalarMatrices/scalarSquareMatrixTemplates.C
index cc1315199..524beae87 100644
--- a/src/foam/matrices/scalarMatrices/scalarSquareMatrixTemplates.C
+++ b/src/foam/matrices/scalarMatrices/scalarSquareMatrixTemplates.C
@@ -38,13 +38,13 @@ void Foam::scalarSquareMatrix::solve
     label n = tmpMatrix.n();
 
     // Elimination
-    for (register label i=0; i<n; i++)
+    for (label i=0; i<n; i++)
     {
         label iMax = i;
         scalar largestCoeff = mag(tmpMatrix[iMax][i]);
 
         // Swap entries around to find a good pivot
-        for (register label j=i+1; j<n; j++)
+        for (label j=i+1; j<n; j++)
         {
             if (mag(tmpMatrix[j][i]) > largestCoeff)
             {
@@ -57,7 +57,7 @@ void Foam::scalarSquareMatrix::solve
         {
             //Info<< "Pivoted on " << i << " " << iMax << endl;
 
-            for (register label k=i; k<n; k++)
+            for (label k=i; k<n; k++)
             {
                 Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
             }
@@ -73,11 +73,11 @@ void Foam::scalarSquareMatrix::solve
         }
 
         // Reduce to upper triangular form
-        for (register label j=i+1; j<n; j++)
+        for (label j=i+1; j<n; j++)
         {
             sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
 
-            for (register label k=n-1; k>=i; k--)
+            for (label k=n-1; k>=i; k--)
             {
                 tmpMatrix[j][k] -=
                     tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i];
@@ -86,11 +86,11 @@ void Foam::scalarSquareMatrix::solve
     }
 
     // Back-substitution
-    for (register label j=n-1; j>=0; j--)
+    for (label j=n-1; j>=0; j--)
     {
         Type ntempvec = pTraits<Type>::zero;
 
-        for (register label k=j+1; k<n; k++)
+        for (label k=j+1; k<n; k++)
         {
             ntempvec += tmpMatrix[j][k]*sourceSol[k];
         }
@@ -125,7 +125,7 @@ void Foam::scalarSquareMatrix::LUBacksubstitute
 
     label ii = 0;
 
-    for (register label i=0; i<n; i++)
+    for (label i=0; i<n; i++)
     {
         label ip = pivotIndices[i];
         Type sum = sourceSol[ip];
@@ -147,12 +147,12 @@ void Foam::scalarSquareMatrix::LUBacksubstitute
         sourceSol[i] = sum;
     }
 
-    for (register label i=n-1; i>=0; i--)
+    for (label i=n-1; i>=0; i--)
     {
         Type sum = sourceSol[i];
         const scalar* __restrict__ luMatrixi = luMatrix[i];
 
-        for (register label j=i+1; j<n; j++)
+        for (label j=i+1; j<n; j++)
         {
             sum -= luMatrixi[j]*sourceSol[j];
         }
diff --git a/src/foam/meshes/meshShapes/cell/cell.C b/src/foam/meshes/meshShapes/cell/cell.C
index d02dbf71c..1de2e0018 100644
--- a/src/foam/meshes/meshShapes/cell/cell.C
+++ b/src/foam/meshes/meshShapes/cell/cell.C
@@ -81,7 +81,7 @@ Foam::labelList Foam::cell::labels(const unallocFaceList& f) const
 
             bool found = false;
 
-            for (register label checkI = 0; checkI < maxVert; checkI++)
+            for (label checkI = 0; checkI < maxVert; checkI++)
             {
                 if (curPoint == p[checkI])
                 {
diff --git a/src/foam/meshes/meshShapes/face/face.C b/src/foam/meshes/meshShapes/face/face.C
index 78119073f..1d5dcebd9 100644
--- a/src/foam/meshes/meshShapes/face/face.C
+++ b/src/foam/meshes/meshShapes/face/face.C
@@ -503,7 +503,7 @@ Foam::point Foam::face::centre(const pointField& meshPoints) const
     label nPoints = size();
 
     point centrePoint = point::zero;
-    for (register label pI = 0; pI < nPoints; pI++)
+    for (label pI = 0; pI < nPoints; pI++)
     {
         centrePoint += meshPoints[operator[](pI)];
     }
@@ -512,7 +512,7 @@ Foam::point Foam::face::centre(const pointField& meshPoints) const
     scalar sumA = 0;
     vector sumAc = vector::zero;
 
-    for (register label pI = 0; pI < nPoints; pI++)
+    for (label pI = 0; pI < nPoints; pI++)
     {
         const point& nextPoint = meshPoints[operator[]((pI + 1) % nPoints)];
 
@@ -567,7 +567,7 @@ Foam::vector Foam::face::normal(const pointField& p) const
 
     label nPoints = size();
 
-    register label pI;
+    label pI;
 
     point centrePoint = vector::zero;
     for (pI = 0; pI < nPoints; pI++)
@@ -690,7 +690,7 @@ Foam::scalar Foam::face::sweptVol
     point nextOldPoint = centreOldPoint;
     point nextNewPoint = centreNewPoint;
 
-    register label pI;
+    label pI;
 
     for (pI = 0; pI < nPoints; pI++)
     {
diff --git a/src/foam/meshes/polyMesh/polyMeshFromShapeMesh.C b/src/foam/meshes/polyMesh/polyMeshFromShapeMesh.C
index b75b8150e..24ffd2c57 100644
--- a/src/foam/meshes/polyMesh/polyMeshFromShapeMesh.C
+++ b/src/foam/meshes/polyMesh/polyMeshFromShapeMesh.C
@@ -78,7 +78,7 @@ Foam::labelList Foam::polyMesh::facePatchFaceCells
     const label patchID
 ) const
 {
-    register bool found;
+    bool found;
 
     labelList FaceCells(patchFaces.size());
 
diff --git a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshEdges.C b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshEdges.C
index dc14ad668..53a889d9e 100644
--- a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshEdges.C
+++ b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshEdges.C
@@ -60,7 +60,7 @@ meshEdges
     // create the storage
     labelList meshEdges(PatchEdges.size());
 
-    register bool found = false;
+    bool found = false;
 
     // get reference to the points on the patch
     const labelList& pp = meshPoints();
diff --git a/src/foam/primitives/BlockCoeff/blockCoeffs.H b/src/foam/primitives/BlockCoeff/blockCoeffs.H
index 776023c49..adf85beb0 100644
--- a/src/foam/primitives/BlockCoeff/blockCoeffs.H
+++ b/src/foam/primitives/BlockCoeff/blockCoeffs.H
@@ -61,7 +61,7 @@ scalar sumMult
 
     typename BlockCoeff<Type>::multiply mult;
 
-    for (register label i = 0; i < f1.size(); i++)
+    for (label i = 0; i < f1.size(); i++)
     {
         result += mult(f1[i], f2[i]);
     }
diff --git a/src/foam/primitives/VectorN/DiagTensorNI.H b/src/foam/primitives/VectorN/DiagTensorNI.H
index 6a1147e7a..13b8007e3 100644
--- a/src/foam/primitives/VectorN/DiagTensorNI.H
+++ b/src/foam/primitives/VectorN/DiagTensorNI.H
@@ -502,8 +502,8 @@ inline DiagTensorN<Cmpt, length> diag(const DiagTensorN<Cmpt, length>& dt)
 // template <class Cmpt, int length>
 // inline Cmpt sum(const DiagTensorN<Cmpt, length>& dt)
 // {
-//     Cmpt result=Cmpt::zero;
-//     for(register label i=0; i<DiagTensorN<Cmpt, length>::nComponents; i++)
+//     Cmpt result = Cmpt::zero;
+//     for(label i = 0; i < DiagTensorN<Cmpt, length>::nComponents; i++)
 //     {
 //         result += dt[i];
 //     }
diff --git a/src/foam/primitives/VectorN/SphericalTensorN.H b/src/foam/primitives/VectorN/SphericalTensorN.H
index 05e0e42d2..b1cedf4fb 100644
--- a/src/foam/primitives/VectorN/SphericalTensorN.H
+++ b/src/foam/primitives/VectorN/SphericalTensorN.H
@@ -101,10 +101,16 @@ public:
         // Access
 
             //- Return component
-            inline const Cmpt& ii() const;
+            inline const Cmpt& ii() const
+            {
+                return this->v_[II];
+            }
 
             //- Return access to component
-            inline Cmpt& ii();
+            inline Cmpt& ii()
+            {
+                return this->v_[II];
+            }
 
 
         //- Diagonal
diff --git a/src/foam/primitives/VectorN/VectorNI.H b/src/foam/primitives/VectorN/VectorNI.H
index efbc453c2..be300e65f 100644
--- a/src/foam/primitives/VectorN/VectorNI.H
+++ b/src/foam/primitives/VectorN/VectorNI.H
@@ -210,7 +210,7 @@ scaleRow(const VectorN<Cmpt, length>& v, const Cmpt& c)
 // inline Cmpt sum(const VectorN<Cmpt, length>& v)
 // {
 //     Cmpt result = pTraits<Cmpt>::zero;
-//     for(register label i=0; i<VectorN<Cmpt, length>::nComponents; i++)
+//     for(label i = 0; i < VectorN<Cmpt, length>::nComponents; i++)
 //     {
 //         result += v[i];
 //     }
diff --git a/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H b/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
index b0db92d85..f2be4e4dd 100644
--- a/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
+++ b/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
@@ -43,7 +43,7 @@ inline void contractScalar(Cmpt& result, const VectorN<Cmpt, length>& t)
 {
     result = pTraits<Cmpt>::zero;
 
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result += t[i];
     }
@@ -69,8 +69,8 @@ inline void contractScalar(Cmpt& result, const TensorN<Cmpt, length>& t)
 {
     result = pTraits<Cmpt>::zero;
 
-    register int j = 0;
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    int j = 0;
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result += t[j];
         j += TensorN<Cmpt, length>::rowLength + 1;
@@ -97,7 +97,7 @@ inline void contractScalar(Cmpt& result, const DiagTensorN<Cmpt, length>& t)
 {
     result = pTraits<Cmpt>::zero;
 
-    for (register int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
     {
         result += t[i];
     }
@@ -146,8 +146,8 @@ inline void contractLinear
     const TensorN<Cmpt, length>& t
 )
 {
-    register int j = 0;
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    int j = 0;
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result[i] = t[j];
         j += TensorN<Cmpt, length>::rowLength + 1;
@@ -173,7 +173,7 @@ inline void contractLinear
     const DiagTensorN<Cmpt, length>& t
 )
 {
-    for (register int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
     {
         result[i] = t[i];
     }
@@ -200,7 +200,7 @@ inline void contractLinear
 {
     for
     (
-        register int i = 0;
+        int i = 0;
         i < SphericalTensorN<Cmpt, length>::rowLength;
         i++
     )
@@ -227,7 +227,7 @@ inline VectorN<Cmpt, length> contractLinear
 template <class Cmpt, int length>
 inline void expandScalar(VectorN<Cmpt, length>& result, const Cmpt& v)
 {
-    for (register int i = 0; i < VectorN<Cmpt, length>::nComponents; i++)
+    for (int i = 0; i < VectorN<Cmpt, length>::nComponents; i++)
     {
         result[i] = v;
     }
@@ -240,8 +240,8 @@ inline void expandScalar(TensorN<Cmpt, length>& result, const Cmpt& v)
 {
     result = TensorN<Cmpt, length>::zero;
 
-    register int j = 0;
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    int j = 0;
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result[j] = v;
         j += TensorN<Cmpt, length>::rowLength + 1;
@@ -253,7 +253,7 @@ inline void expandScalar(TensorN<Cmpt, length>& result, const Cmpt& v)
 template <class Cmpt, int length>
 inline void expandScalar(DiagTensorN<Cmpt, length>& result, const Cmpt& v)
 {
-    for (register int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
     {
         result[i] = v;
     }
@@ -278,8 +278,8 @@ inline void expandLinear
 {
     result = TensorN<Cmpt, length>::zero;
 
-    register int j = 0;
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    int j = 0;
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result[j] = v[i];
         j += TensorN<Cmpt, length>::rowLength + 1;
@@ -294,7 +294,7 @@ inline void expandLinear
     const VectorN<Cmpt, length>& v
 )
 {
-    for (register int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
     {
         result[i] = v[i];
     }
@@ -311,7 +311,7 @@ inline void expandLinear
 {
     result[0] = pTraits<Cmpt>::zero;
 
-    for (register int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < DiagTensorN<Cmpt, length>::rowLength; i++)
     {
         result[0] += v[i];
     }
@@ -330,9 +330,9 @@ inline void sumToDiag
 {
     result = VectorN<Cmpt, length>::zero;
 
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
-        for (register int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
+        for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
         {
             result[i] += t(i, j);
         }
@@ -360,9 +360,9 @@ inline void sumMagToDiag
 {
     result = VectorN<Cmpt, length>::zero;
 
-    for (register int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
-        for (register int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
+        for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
         {
             result[i] += Foam::mag(t(i, j));
         }
diff --git a/src/foam/primitives/strings/fileName/fileName.C b/src/foam/primitives/strings/fileName/fileName.C
index 89a4e1e88..c34a8ef02 100644
--- a/src/foam/primitives/strings/fileName/fileName.C
+++ b/src/foam/primitives/strings/fileName/fileName.C
@@ -105,7 +105,7 @@ bool Foam::fileName::clean()
         src < maxLen;
     )
     {
-        register char c = operator[](src++);
+        char c = operator[](src++);
 
         if (prev == '/')
         {
diff --git a/src/foam/primitives/strings/string/foamString.C b/src/foam/primitives/strings/string/foamString.C
index 208921aca..1783e6e4e 100644
--- a/src/foam/primitives/strings/string/foamString.C
+++ b/src/foam/primitives/strings/string/foamString.C
@@ -38,7 +38,7 @@ const Foam::string Foam::string::null;
 
 Foam::string::size_type Foam::string::count(const char c) const
 {
-    register size_type cCount = 0;
+    size_type cCount = 0;
 
     for (const_iterator iter = begin(); iter != end(); ++iter)
     {
@@ -105,10 +105,10 @@ bool Foam::string::removeRepeated(const char character)
 
     if (character && find(character) != npos)
     {
-        register string::size_type nChar=0;
+        string::size_type nChar=0;
         iterator iter2 = begin();
 
-        register char prev = 0;
+        char prev = 0;
 
         for
         (
@@ -117,7 +117,7 @@ bool Foam::string::removeRepeated(const char character)
             iter1++
         )
         {
-            register char c = *iter1;
+            char c = *iter1;
 
             if (prev == c && c == character)
             {
diff --git a/src/foam/primitives/strings/string/foamStringI.H b/src/foam/primitives/strings/string/foamStringI.H
index 51f701283..7fb6e4018 100644
--- a/src/foam/primitives/strings/string/foamStringI.H
+++ b/src/foam/primitives/strings/string/foamStringI.H
@@ -79,7 +79,7 @@ inline bool Foam::string::stripInvalid(string& str)
 {
     if (!valid<String>(str))
     {
-        register size_type nValid = 0;
+        size_type nValid = 0;
         iterator iter2 = str.begin();
 
         for
@@ -89,7 +89,7 @@ inline bool Foam::string::stripInvalid(string& str)
             iter1++
         )
         {
-            register char c = *iter1;
+            char c = *iter1;
 
             if (String::valid(c))
             {
diff --git a/src/lduSolvers/amg/amgPolicy/selectionAmgPolicy.C b/src/lduSolvers/amg/amgPolicy/selectionAmgPolicy.C
index 5b27aa5bd..55dab38eb 100644
--- a/src/lduSolvers/amg/amgPolicy/selectionAmgPolicy.C
+++ b/src/lduSolvers/amg/amgPolicy/selectionAmgPolicy.C
@@ -175,12 +175,12 @@ void Foam::selectionAmgPolicy::calcCoarsening()
     // Select the strongest coefficient in each row
     // (There is no need to operate in rows in the first 2 loops, but the
     // code is left row-wise for clarity)
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(diag[i]);
 
         // Do lower triangle coefficient for the row
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij = mag(min(signDiag*lower[losortAddr[jp]], 0));
 
@@ -191,7 +191,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
         }
 
         // Do upper triangle coefficient for the row
-        for (register label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
+        for (label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
         {
             const scalar magAij = mag(min(signDiag*upper[ip], 0));
 
@@ -209,12 +209,12 @@ void Foam::selectionAmgPolicy::calcCoarsening()
     // (epsilonStrongCoeff_)
     labelList strongCoeffCounter(nRows, 0);
 
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(diag[i]);
 
         // Do lower triangle coefficient for the row
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij = mag(min(signDiag*lower[losortAddr[jp]], 0));
 
@@ -251,12 +251,12 @@ void Foam::selectionAmgPolicy::calcCoarsening()
     strongCoeffCounter = 0;
 
     // Fill in the empty arrays of the strong matrix
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         const scalar signDiag = sign(diag[i]);
 
         // Do lower triangle coefficient for the row
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             const scalar magAij = mag(min(signDiag*lower[losortAddr[jp]], 0));
 
@@ -518,7 +518,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
     // Start row assembly
     pRow[0] = 0;
 
-    for (register label i = 0; i < nRows; i++)
+    for (label i = 0; i < nRows; i++)
     {
         label rowCount = pRow[i];
 
@@ -529,7 +529,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
         signDii = sign(diag[i]);
 
         // Do lower triangle coefficient for the row
-        for (register label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
+        for (label jp = losortStart[i]; jp < losortStart[i + 1]; jp++)
         {
             // Adjust sign of off-diag coeff
             Dij = signDii*lower[losortAddr[jp]];
@@ -541,7 +541,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
         }
 
         // Do upper triangle coefficient for the row
-        for (register label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
+        for (label ip = ownerStart[i]; ip < ownerStart[i + 1]; ip++)
         {
             // Adjust sign of off-diag coeff
             Dij = signDii*upper[ip];
@@ -561,7 +561,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
 
             for
             (
-                register label sip = strongRow[i];
+                label sip = strongRow[i];
                 sip < strongRow[i + 1];
                 sip++
             )
@@ -577,7 +577,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
             // Calculate interpolation weights.
             for
             (
-                register label sip = strongRow[i];
+                label sip = strongRow[i];
                 sip < strongRow[i + 1];
                 sip++
             )
@@ -644,7 +644,7 @@ void Foam::selectionAmgPolicy::calcCoarsening()
         const labelList& prolongationRow = Pptr_->crAddr().rowStart();
         const scalarField& prolongationCoeff = Pptr_->coeffs();
 
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             if (prolongationRow[rowI] == prolongationRow[rowI + 1])
             {
diff --git a/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C b/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C
index cd34c6ae9..77dfd036f 100644
--- a/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C
+++ b/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C
@@ -165,7 +165,7 @@ void Foam::CholeskyPrecon::precondition
 
         // Note: multiplication over-write x: no need to initialise
         // HJ, and VV, 19/Jun/2017
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             xPtr[rowI] = bPtr[rowI]*preconDiagPtr[rowI];
         }
@@ -190,7 +190,7 @@ void Foam::CholeskyPrecon::precondition
         const label nCoeffs = matrix_.upper().size();
 
         // Forward sweep
-        for (register label coeffI = 0; coeffI < nCoeffs; coeffI++)
+        for (label coeffI = 0; coeffI < nCoeffs; coeffI++)
         {
             xPtr[uPtr[coeffI]] -=
                 preconDiagPtr[uPtr[coeffI]]*
@@ -198,7 +198,7 @@ void Foam::CholeskyPrecon::precondition
         }
 
         // Reverse sweep
-        for (register label coeffI = nCoeffs - 1; coeffI >= 0; coeffI--)
+        for (label coeffI = nCoeffs - 1; coeffI >= 0; coeffI--)
         {
             xPtr[lPtr[coeffI]] -=
                 preconDiagPtr[lPtr[coeffI]]*
diff --git a/src/lduSolvers/lduPrecon/ILU0/ILU0.C b/src/lduSolvers/lduPrecon/ILU0/ILU0.C
index 06db8143c..31676ce51 100644
--- a/src/lduSolvers/lduPrecon/ILU0/ILU0.C
+++ b/src/lduSolvers/lduPrecon/ILU0/ILU0.C
@@ -202,7 +202,7 @@ void Foam::ILU0::precondition
 
         // Note: multiplication over-write x: no need to initialise
         // HJ, and VV, 19/Jun/2017
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             xPtr[rowI] = bPtr[rowI]*preconDiagPtr[rowI];
         }
@@ -308,7 +308,7 @@ void Foam::ILU0::preconditionT
 
         // Note: multiplication over-write x: no need to initialise
         // HJ, and VV, 19/Jun/2017
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             xPtr[rowI] = bPtr[rowI]*preconDiagPtr[rowI];
         }
diff --git a/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C b/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C
index ca830b52a..fb8ecc1d7 100644
--- a/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C
+++ b/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C
@@ -94,12 +94,12 @@ void Foam::ILUC0::calcFactorization()
 
         // Define start and end face of this row/column, and number of non zero
         // off diagonal entries
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -109,7 +109,7 @@ void Foam::ILUC0::calcFactorization()
             zDiag_ = diagPtr[rowI];
 
             // Initialize temporary working row field
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper), w addressed by neighbour of face (row index for
@@ -125,13 +125,13 @@ void Foam::ILUC0::calcFactorization()
             // Lower coeff loop (first i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortIndex = lsrPtr[faceLsrI];
+                const label losortIndex = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortIndex];
@@ -140,14 +140,14 @@ void Foam::ILUC0::calcFactorization()
                 zDiag_ -= lowerPtr[losortIndex]*upperPtr[losortIndex];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortIndex + 1 = start)
-                    register label faceI = losortIndex + 1;
+                    label faceI = losortIndex + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -162,10 +162,10 @@ void Foam::ILUC0::calcFactorization()
             diagRowI = 1.0/zDiag_;
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -182,23 +182,23 @@ void Foam::ILUC0::calcFactorization()
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortIndex = lsrPtr[faceLsrI];
+                const label losortIndex = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortIndex];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortIndex + 1;
+                    label faceI = losortIndex + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -303,8 +303,8 @@ void Foam::ILUC0::precondition
         // Initialize x field
         x = b;
 
-        register label losortIndexI;
-        register label rowI;
+        label losortIndexI;
+        label rowI;
 
         // Forward substitution loop
         forAll (preconLower_, coeffI)
@@ -379,8 +379,8 @@ void Foam::ILUC0::preconditionT
             x[i] = b[i]*preconDiag_[i];
         }
 
-        register label losortIndexI;
-        register label rowI;
+        label losortIndexI;
+        label rowI;
 
         // Forward substitution loop
         forAll (preconUpper_, coeffI)
diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
index 628f83c26..4ce1c5602 100644
--- a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
+++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
@@ -98,12 +98,12 @@ void Foam::ILUCp::calcFactorization()
 
         // Define start and end face ("virtual" face when extended addressing
         // is used) of this row/column.
-        register label fStart, fEnd, fLsrStart, fLsrEnd;
+        label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
 
         // Row by row loop (k - loop).
-        for (register label rowI = 0; rowI < nRows; ++rowI)
+        for (label rowI = 0; rowI < nRows; ++rowI)
         {
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
@@ -113,7 +113,7 @@ void Foam::ILUCp::calcFactorization()
             zDiag_ = diagPtr[rowI];
 
             // Initialize temporary working row field
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Note: z addressed by neighbour of face (column index for
                 // upper), w addressed by neighbour of face (row index for
@@ -129,13 +129,13 @@ void Foam::ILUCp::calcFactorization()
             // Lower/upper coeff loop (i - loop)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
@@ -144,14 +144,14 @@ void Foam::ILUCp::calcFactorization()
                 zDiag_ -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 // Upper coeff loop (additional loop to avoid checking the
                 // existence of certain upper coeffs)
                 for
                 (
                     // Diagonal is already updated (losortCoeff + 1 = start)
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -166,10 +166,10 @@ void Foam::ILUCp::calcFactorization()
             diagRowI = 1.0/zDiag_;
 
             // Index for updating L and U
-            register label zwIndex;
+            label zwIndex;
 
             // Update upper and lower coeffs
-            for (register label faceI = fStart; faceI < fEnd; ++faceI)
+            for (label faceI = fStart; faceI < fEnd; ++faceI)
             {
                 // Get index for current face
                 zwIndex = uPtr[faceI];
@@ -186,23 +186,23 @@ void Foam::ILUCp::calcFactorization()
             // this step (for this row and column)
             for
             (
-                register label faceLsrI = fLsrStart;
+                label faceLsrI = fLsrStart;
                 faceLsrI < fLsrEnd;
                 ++faceLsrI
             )
             {
                 // Get losort coefficient for this face
-                const register label losortCoeff = lsrPtr[faceLsrI];
+                const label losortCoeff = lsrPtr[faceLsrI];
 
                 // Get corresponding row index for upper (i label)
                 const label i = lPtr[losortCoeff];
 
                 // Get end of row for cell i
-                const register label fEndRowi = ownStartPtr[i + 1];
+                const label fEndRowi = ownStartPtr[i + 1];
 
                 for
                 (
-                    register label faceI = losortCoeff + 1;
+                    label faceI = losortCoeff + 1;
                     faceI < fEndRowi;
                     ++faceI
                 )
@@ -291,8 +291,8 @@ void Foam::ILUCp::precondition
         // Initialize x field
         x = b;
 
-        register label losortCoeffI;
-        register label rowI;
+        label losortCoeffI;
+        label rowI;
 
         // Forward substitution loop
         forAll (lower, coeffI)
@@ -370,8 +370,8 @@ void Foam::ILUCp::preconditionT
             x[i] = b[i]*preconDiag_[i];
         }
 
-        register label losortCoeffI;
-        register label rowI;
+        label losortCoeffI;
+        label rowI;
 
         // Forward substitution loop
         forAll (upper, coeffI)
diff --git a/src/lduSolvers/lduPrecon/symGaussSeidelPrecon/symGaussSeidelPrecon.C b/src/lduSolvers/lduPrecon/symGaussSeidelPrecon/symGaussSeidelPrecon.C
index 961f216a5..c68d79059 100644
--- a/src/lduSolvers/lduPrecon/symGaussSeidelPrecon/symGaussSeidelPrecon.C
+++ b/src/lduSolvers/lduPrecon/symGaussSeidelPrecon/symGaussSeidelPrecon.C
@@ -112,7 +112,7 @@ void Foam::symGaussSeidelPrecon::precondition
 
         const label nRows = x.size();
 
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             xPtr[rowI] = bPtr[rowI]/diagPtr[rowI];
         }
@@ -164,7 +164,7 @@ void Foam::symGaussSeidelPrecon::precondition
         }
 
         // Forward sweep
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             // lRow is equal to rowI
             scalar& curX = xPtr[rowI];
@@ -177,7 +177,7 @@ void Foam::symGaussSeidelPrecon::precondition
             fEnd = ownStartPtr[rowI + 1];
 
             // Accumulate the owner product side
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 curX -= upperPtr[curCoeff]*xPtr[uPtr[curCoeff]];
             }
@@ -186,14 +186,14 @@ void Foam::symGaussSeidelPrecon::precondition
             curX /= diagPtr[rowI];
 
             // Distribute the neighbour side using current x
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 bPrimePtr[uPtr[curCoeff]] -= lowerPtr[curCoeff]*curX;
             }
         }
 
         // Reverse sweep
-        for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+        for (label rowI = nRows - 1; rowI >= 0; rowI--)
         {
             // lRow is equal to rowI
             scalar& curX = xPtr[rowI];
@@ -206,7 +206,7 @@ void Foam::symGaussSeidelPrecon::precondition
             fEnd = ownStartPtr[rowI + 1];
 
             // Accumulate the owner product side
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 curX -= upperPtr[curCoeff]*xPtr[uPtr[curCoeff]];
             }
@@ -279,7 +279,7 @@ void Foam::symGaussSeidelPrecon::preconditionT
         }
 
         // Forward sweep
-        for (register label rowI = 0; rowI < nRows; rowI++)
+        for (label rowI = 0; rowI < nRows; rowI++)
         {
             // lRow is equal to rowI
             scalar& curX = xPtr[rowI];
@@ -292,7 +292,7 @@ void Foam::symGaussSeidelPrecon::preconditionT
             fEnd = ownStartPtr[rowI + 1];
 
             // Accumulate the owner product side
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 // Transpose multiplication.  HJ, 10/Jul/2007
                 curX -= lowerPtr[curCoeff]*xPtr[uPtr[curCoeff]];
@@ -302,7 +302,7 @@ void Foam::symGaussSeidelPrecon::preconditionT
             curX /= diagPtr[rowI];
 
             // Distribute the neighbour side using current x
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 // Transpose multiplication.  HJ, 10/Jul/2007
                 bPrimePtr[uPtr[curCoeff]] -= upperPtr[curCoeff]*curX;
@@ -310,7 +310,7 @@ void Foam::symGaussSeidelPrecon::preconditionT
         }
 
         // Reverse sweep
-        for (register label rowI = nRows - 1; rowI >= 0; rowI--)
+        for (label rowI = nRows - 1; rowI >= 0; rowI--)
         {
             // lRow is equal to rowI
             scalar& curX = xPtr[rowI];
@@ -323,7 +323,7 @@ void Foam::symGaussSeidelPrecon::preconditionT
             fEnd = ownStartPtr[rowI + 1];
 
             // Accumulate the owner product side
-            for (register label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
+            for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)
             {
                 // Transpose multiplication.  HJ, 10/Jul/2007
                 curX -= lowerPtr[curCoeff]*xPtr[uPtr[curCoeff]];
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
index b7cbe707e..a9080982f 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
@@ -256,7 +256,7 @@ void meshToMesh::cellAddresses
     // point than the current cell
 
     // set curCell label to zero (start)
-    register label curCell = 0;
+    label curCell = 0;
 
     // set reference to cell to cell addressing
     const vectorField& centresFrom = fromMesh.cellCentres();
diff --git a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
index 574fd0646..f62948f94 100644
--- a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
+++ b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
@@ -692,7 +692,7 @@ void dirichletNeumann::calcSlavePointPoints(const label slavePatchID)
 
   const faceList& lf = mesh.boundaryMesh()[slavePatchID].localFaces();
 
-  register bool found = false;
+  bool found = false;
 
   forAll (lf, faceI)
     {
diff --git a/src/solidModels/contactModels/normalContactModels/iterativePenalty/iterativePenalty.C b/src/solidModels/contactModels/normalContactModels/iterativePenalty/iterativePenalty.C
index 31b20b7cb..05811448d 100644
--- a/src/solidModels/contactModels/normalContactModels/iterativePenalty/iterativePenalty.C
+++ b/src/solidModels/contactModels/normalContactModels/iterativePenalty/iterativePenalty.C
@@ -455,7 +455,7 @@ iterativePenalty::iterativePenalty
 
   const faceList& lf = mesh.boundaryMesh()[slavePatchID].localFaces();
 
-  register bool found = false;
+  bool found = false;
 
   forAll (lf, faceI)
     {
diff --git a/src/solidModels/contactModels/normalContactModels/standardPenalty/standardPenalty.C b/src/solidModels/contactModels/normalContactModels/standardPenalty/standardPenalty.C
index ff6e0e56e..7f6ce61f3 100644
--- a/src/solidModels/contactModels/normalContactModels/standardPenalty/standardPenalty.C
+++ b/src/solidModels/contactModels/normalContactModels/standardPenalty/standardPenalty.C
@@ -362,7 +362,7 @@ standardPenalty::standardPenalty
 
   const faceList& lf = mesh.boundaryMesh()[slavePatchID].localFaces();
 
-  register bool found = false;
+  bool found = false;
 
   forAll (lf, faceI)
     {
diff --git a/src/solidModels/finiteVolume/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceGrad.C b/src/solidModels/finiteVolume/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceGrad.C
index 66dc4988e..48d0b2ba2 100644
--- a/src/solidModels/finiteVolume/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceGrad.C
+++ b/src/solidModels/finiteVolume/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceGrad.C
@@ -146,8 +146,8 @@ leastSquaresSolidInterfaceGrad<Type>::calcGrad
 
     forAll(own, facei)
     {
-        register label ownFaceI = own[facei];
-        register label neiFaceI = nei[facei];
+        label ownFaceI = own[facei];
+        label neiFaceI = nei[facei];
 
         if (interfaceFacesMap[facei] > -SMALL)
         {
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
index e5578a2b2..b80f1bdcf 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
@@ -71,7 +71,7 @@ Foam::janafThermo<equationOfState>::janafThermo(Istream& is)
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -81,7 +81,7 @@ Foam::janafThermo<equationOfState>::janafThermo(Istream& is)
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -112,7 +112,7 @@ Foam::Ostream& Foam::operator<<
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -124,7 +124,7 @@ Foam::Ostream& Foam::operator<<
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index 6b18e8f89..416fe16dc 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -44,7 +44,7 @@ inline Foam::janafThermo<equationOfState>::janafThermo
     Thigh_(Thigh),
     Tcommon_(Tcommon)
 {
-    for (register label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
+    for (label coefLabel = 0; coefLabel < nCoeffs_; coefLabel++)
     {
         highCpCoeffs_[coefLabel] = highCpCoeffs[coefLabel];
         lowCpCoeffs_[coefLabel] = lowCpCoeffs[coefLabel];
@@ -108,7 +108,7 @@ inline Foam::janafThermo<equationOfState>::janafThermo
     Thigh_(jt.Thigh_),
     Tcommon_(jt.Tcommon_)
 {
-    for (register label coefLabel = 0; coefLabel < nCoeffs_; coefLabel++)
+    for (label coefLabel = 0; coefLabel < nCoeffs_; coefLabel++)
     {
         highCpCoeffs_[coefLabel] = jt.highCpCoeffs_[coefLabel];
         lowCpCoeffs_[coefLabel] = jt.lowCpCoeffs_[coefLabel];
@@ -209,7 +209,7 @@ inline void Foam::janafThermo<equationOfState>::operator+=
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -244,7 +244,7 @@ inline void Foam::janafThermo<equationOfState>::operator-=
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -280,7 +280,7 @@ inline Foam::janafThermo<equationOfState> Foam::operator+
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )
@@ -324,7 +324,7 @@ inline Foam::janafThermo<equationOfState> Foam::operator-
 
     for
     (
-        register label coefLabel=0;
+        label coefLabel = 0;
         coefLabel<janafThermo<equationOfState>::nCoeffs_;
         coefLabel++
     )