diff --git a/run/coFlow/constant/polyMesh/blockMeshDict b/run/coFlow/system/blockMeshDict
similarity index 100%
rename from run/coFlow/constant/polyMesh/blockMeshDict
rename to run/coFlow/system/blockMeshDict
diff --git a/run/counterFlow/constant/polyMesh/blockMeshDict b/run/counterFlow/system/blockMeshDict
similarity index 100%
rename from run/counterFlow/constant/polyMesh/blockMeshDict
rename to run/counterFlow/system/blockMeshDict
diff --git a/run/crossFlow/constant/polyMesh/blockMeshDict b/run/crossFlow/system/blockMeshDict
similarity index 100%
rename from run/crossFlow/constant/polyMesh/blockMeshDict
rename to run/crossFlow/system/blockMeshDict
diff --git a/run/quickTest/constant/polyMesh/blockMeshDict b/run/quickTest/system/blockMeshDict
similarity index 100%
rename from run/quickTest/constant/polyMesh/blockMeshDict
rename to run/quickTest/system/blockMeshDict
diff --git a/run/quickTestStack/constant/polyMesh/blockMeshDict b/run/quickTestStack/system/blockMeshDict
similarity index 100%
rename from run/quickTestStack/constant/polyMesh/blockMeshDict
rename to run/quickTestStack/system/blockMeshDict
diff --git a/src/appSrc/Make/options b/src/appSrc/Make/options
index eb78c5f..3eae2b3 100755
--- a/src/appSrc/Make/options
+++ b/src/appSrc/Make/options
@@ -2,7 +2,7 @@ EXE_INC = -g \
     -I../libSrc/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ 
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
 
 EXE_LIBS = \
     -L$(FOAM_USER_LIBBIN) \
diff --git a/src/appSrc/activationOverpotential.H b/src/appSrc/activationOverpotential.H
index e43375a..454be4d 100644
--- a/src/appSrc/activationOverpotential.H
+++ b/src/appSrc/activationOverpotential.H
@@ -130,17 +130,17 @@
 
 //  Save data for output
 
-    scalarField& etaAIn = etaA.internalField();
+    scalarField& etaAIn = etaA.primitiveFieldRef();
     smearPatchToMesh smearetaA(electrolyteMesh, electrolyteAnodeName);
     etaAIn = smearetaA.smear(etaAnode);
     etaA.correctBoundaryConditions();
-    etaA.boundaryField()[electrolyteAnodeID] == etaAnode;
+    etaA.boundaryFieldRef()[electrolyteAnodeID] == etaAnode;
 
-    scalarField& etaCIn = etaC.internalField();
+    scalarField& etaCIn = etaC.primitiveFieldRef();
     smearPatchToMesh smearetaC(electrolyteMesh, electrolyteAnodeName);
     etaCIn = smearetaC.smear(etaCathode);
     etaC.correctBoundaryConditions();
-    etaC.boundaryField()[electrolyteCathodeID] == etaCathode;
+    etaC.boundaryFieldRef()[electrolyteCathodeID] == etaCathode;
 
 
 
diff --git a/src/appSrc/diffusivityAir.H b/src/appSrc/diffusivityAir.H
index 0e83413..6ab4ec6 100644
--- a/src/appSrc/diffusivityAir.H
+++ b/src/appSrc/diffusivityAir.H
@@ -20,7 +20,7 @@ forAll(airSpecies, a)
     airDiff = 0;         //diffusivity field passed to/from model
 
     //diffSpAir[a] is diffusivity field for specie a, to be used in YEqn
-    scalarField& diffSpAirIn = diffSpAir[a].internalField();
+    scalarField& diffSpAirIn = diffSpAir[a].primitiveFieldRef();
     diffSpAirIn = 0;
 
     if(airSpecies[a].name() != airInertSpecie)
diff --git a/src/appSrc/diffusivityFuel.H b/src/appSrc/diffusivityFuel.H
index 20712f1..f4e170e 100644
--- a/src/appSrc/diffusivityFuel.H
+++ b/src/appSrc/diffusivityFuel.H
@@ -20,7 +20,7 @@ forAll(fuelSpecies, a)
     fuelDiff = 0;         //diffusivity field passed to/from model
 
     //diffSpFuel[a] is diffusivity field for specie a, to be used in YEqn
-    scalarField& diffSpFuelIn = diffSpFuel[a].internalField();
+    scalarField& diffSpFuelIn = diffSpFuel[a].primitiveFieldRef();
     diffSpFuelIn = 0;
 
     if(fuelSpecies[a].name() != fuelInertSpecie)
diff --git a/src/appSrc/electrochemicalHeating.H b/src/appSrc/electrochemicalHeating.H
index 07cc64f..6a14390 100755
--- a/src/appSrc/electrochemicalHeating.H
+++ b/src/appSrc/electrochemicalHeating.H
@@ -3,7 +3,7 @@
     Info << nl << "Calculating electrochemical heating" << endl;
 
     // iEA applied in cell volumes (not patch)
-    scalarField iEA = idensity.internalField();
+    scalarField iEA = idensity.primitiveFieldRef();
 
 
     //scalar Tr = 298.15;  //reference Temperature [K] for enthalpy calculations
@@ -120,7 +120,7 @@
 //        - iEA*V.value()/hE.value();
 //
 
-    scalarField& electrochemicalHeatingIn = electrochemicalHeating.internalField();
+    scalarField& electrochemicalHeatingIn = electrochemicalHeating.primitiveFieldRef();
     electrochemicalHeatingIn =
     (
         -hFormSum*volMolRate - hSource - iEA*voltage/hE.value()
@@ -128,7 +128,7 @@
     electrochemicalHeating.correctBoundaryConditions();
 
     Info << "min,mean,max(electrochemicalHeating): "
-         << gMin(electrochemicalHeating.internalField()) << "  "
-         << gAverage(electrochemicalHeating.internalField()) << "  "
-         << gMax(electrochemicalHeating.internalField()) << endl;
+         << gMin(electrochemicalHeating.primitiveFieldRef()) << "  "
+         << gAverage(electrochemicalHeating.primitiveFieldRef()) << "  "
+         << gMax(electrochemicalHeating.primitiveFieldRef()) << endl;
 
diff --git a/src/appSrc/idensity.H b/src/appSrc/idensity.H
index 7bb1e5b..82b4c8b 100644
--- a/src/appSrc/idensity.H
+++ b/src/appSrc/idensity.H
@@ -1,13 +1,13 @@
 // Current density
 
-    scalarField& idensityIn = idensity.internalField();
+    scalarField& idensityIn = idensity.primitiveFieldRef();
     smearPatchToMesh smearIdensity(electrolyteMesh, electrolyteAnodeName);
     idensityIn = smearIdensity.smear(i);
     idensity.correctBoundaryConditions();
 
 // Nernst Potential
 
-    scalarField& NernstPotIn = NernstPot.internalField();
+    scalarField& NernstPotIn = NernstPot.primitiveFieldRef();
     smearPatchToMesh smearNernstPot(electrolyteMesh, electrolyteAnodeName);
     NernstPotIn = smearNernstPot.smear(Nernst);
     NernstPot.correctBoundaryConditions();
diff --git a/src/appSrc/kAir.H b/src/appSrc/kAir.H
index d89dbb3..a6b6922 100644
--- a/src/appSrc/kAir.H
+++ b/src/appSrc/kAir.H
@@ -1,7 +1,7 @@
 // obtain k(air mixture) from k(species components)
 
 {
-    scalarField& kAIn = kAir.internalField();
+    scalarField& kAIn = kAir.primitiveFieldRef();
     kAIn = 0;
 
     // initially using a simple linear combination weighted by mass fraction
diff --git a/src/appSrc/kFuel.H b/src/appSrc/kFuel.H
index 0ded78f..a102afb 100644
--- a/src/appSrc/kFuel.H
+++ b/src/appSrc/kFuel.H
@@ -1,7 +1,7 @@
 // obtain k(fuel mixture) from k(species components)
 
 {
-    scalarField& kFIn = kFuel.internalField();
+    scalarField& kFIn = kFuel.primitiveFieldRef();
     kFIn = 0;
 
     // initially using a simple linear combination weighted by mass fraction
diff --git a/src/appSrc/mapAirToCell.H b/src/appSrc/mapAirToCell.H
index a0ec54b..67c0baa 100755
--- a/src/appSrc/mapAirToCell.H
+++ b/src/appSrc/mapAirToCell.H
@@ -5,18 +5,18 @@
         // mass based heat capacity cp
         // cp [J/kg-K] = Cp/M where Cp = [J/mol-K], M = [kg/mol]
 
-        scalarField& cp = cpAir.internalField();
+        scalarField& cp = cpAir.primitiveFieldRef();
         cp = 0;
 
         forAll(airSpecies, i)
         {
             scalar rMWi = 1.0/(airSpecies[i].MW()*1e-3);    // 1/[kg/mol]
-            scalarField cpi = molarCpAir[i].polyVal(Tair.internalField())*rMWi;
-            cp += Yair[i].internalField()*cpi;
+            scalarField cpi = molarCpAir[i].polyVal(Tair.primitiveFieldRef())*rMWi;
+            cp += Yair[i].primitiveFieldRef()*cpi;
         }
         cpAir.correctBoundaryConditions();
 
-        scalarField rhoCpF = cp*rhoAir.internalField();
+        scalarField rhoCpF = cp*rhoAir.primitiveFieldRef();
 
         Info << nl;
         Info<< "min mean max rhoCpF(air) = " << gMin(rhoCpF) << " "
@@ -27,7 +27,7 @@
 
 
         //scalarField kIn(airMesh.nCells(), kA.value());
-        scalarField& kIn = kAir.internalField();
+        scalarField& kIn = kAir.primitiveFieldRef();
 
         forAll(cathodeZones, iz)
         {
@@ -69,12 +69,12 @@
         //
 
         scalarField rhoCpPhiF =
-            //CpAir.value()*phiAir.internalField();
-            linearInterpolate(cpAir)*phiAir.internalField();
+            //CpAir.value()*phiAir.primitiveFieldRef();
+            linearInterpolate(cpAir)*phiAir.primitiveFieldRef();
 
         phiCellIn.rmap
         (
-            phiAir.internalField()*airInternalFaceMask,
+            phiAir.primitiveFieldRef()*airInternalFaceMask,
             airInternalFaceMap
         );
 
@@ -131,14 +131,14 @@
                 curFpm -= mesh.boundary()
                     [airPatchesMap[patchI]].patch().start();
 
-                phiCell.boundaryField()[airPatchesMap[patchI]].
+                phiCell.boundaryFieldRef()[airPatchesMap[patchI]].
                     scalarField::rmap
                     (
                         phiAir.boundaryField()[patchI]*curMask,
                         curFpm
                     );
 
-                rhoCpPhiCell.boundaryField()[airPatchesMap[patchI]].
+                rhoCpPhiCell.boundaryFieldRef()[airPatchesMap[patchI]].
                     scalarField::rmap
                     (
                         (
diff --git a/src/appSrc/mapFuelToCell.H b/src/appSrc/mapFuelToCell.H
index 74be52d..1427384 100755
--- a/src/appSrc/mapFuelToCell.H
+++ b/src/appSrc/mapFuelToCell.H
@@ -5,7 +5,7 @@
         // mass based heat capacity cp
         // cp [J/kg-K] = Cp/M where Cp = [J/mol-K], M = [kg/mol]
 
-        scalarField& cp = cpFuel.internalField();
+        scalarField& cp = cpFuel.primitiveFieldRef();
         cp = 0;
 
         forAll(fuelSpecies, i)
@@ -17,13 +17,13 @@
             scalar rMWi = 1.0/(fuelSpecies[i].MW()*1e-3);    // 1/[kg/mol]
             scalarField cpi =
             (
-                 molarCpFuel[i].polyVal(Tfuel.internalField())*rMWi
+                 molarCpFuel[i].polyVal(Tfuel.primitiveFieldRef())*rMWi
             );
-            cp += Yfuel[i].internalField()*cpi;
+            cp += Yfuel[i].primitiveFieldRef()*cpi;
         }
         cpFuel.correctBoundaryConditions();
 
-        scalarField rhoCpF = cp*rhoFuel.internalField();
+        scalarField rhoCpF = cp*rhoFuel.primitiveFieldRef();
 
         Info<< "min mean max rhoCpF(fuel) = " << gMin(rhoCpF) << " "
             << gAverage(rhoCpF) << " " << gMax(rhoCpF) << nl;
@@ -34,7 +34,7 @@
 
 
         //scalarField kIn(fuelMesh.nCells(), kF.value());
-        scalarField& kIn = kFuel.internalField();
+        scalarField& kIn = kFuel.primitiveFieldRef();
 
         forAll(anodeZones, iz)
         {
@@ -76,12 +76,12 @@
         //
 
         scalarField rhoCpPhiF =
-            //CpFuel.value()*phiFuel.internalField();
-            linearInterpolate(cpFuel)*phiFuel.internalField();
+            //CpFuel.value()*phiFuel.primitiveFieldRef();
+            linearInterpolate(cpFuel)*phiFuel.primitiveFieldRef();
 
         phiCellIn.rmap
         (
-            phiFuel.internalField()*fuelInternalFaceMask,
+            phiFuel.primitiveFieldRef()*fuelInternalFaceMask,
             fuelInternalFaceMap
         );
 
@@ -138,14 +138,14 @@
                 curFpm -= mesh.boundary()
                     [fuelPatchesMap[patchI]].patch().start();
 
-                phiCell.boundaryField()[fuelPatchesMap[patchI]].
+                phiCell.boundaryFieldRef()[fuelPatchesMap[patchI]].
                     scalarField::rmap
                     (
                         phiFuel.boundaryField()[patchI]*curMask,
                         curFpm
                     );
 
-                rhoCpPhiCell.boundaryField()[fuelPatchesMap[patchI]].
+                rhoCpPhiCell.boundaryFieldRef()[fuelPatchesMap[patchI]].
                     scalarField::rmap
                     (
                         (
diff --git a/src/appSrc/mapToCell.H b/src/appSrc/mapToCell.H
index 0e5881b..969871d 100755
--- a/src/appSrc/mapToCell.H
+++ b/src/appSrc/mapToCell.H
@@ -2,19 +2,19 @@
     // Map to cell.  In practice, fields for rho, Cp and k will come from
     // components (air, fuel, electrolyte) but currently they are constant
 
-    scalarField& rhoCpCellIn = rhoCpCell.internalField();
+    scalarField& rhoCpCellIn = rhoCpCell.primitiveFieldRef();
     rhoCpCellIn = 0.0;
 
-    scalarField& kCellIn = kCell.internalField();
+    scalarField& kCellIn = kCell.primitiveFieldRef();
     kCellIn = 0.0;
 
-    scalarField& TsourceCellIn = TsourceCell.internalField();
+    scalarField& TsourceCellIn = TsourceCell.primitiveFieldRef();
     TsourceCellIn = 0.0;
 
-    scalarField& phiCellIn = phiCell.internalField();
+    scalarField& phiCellIn = phiCell.primitiveFieldRef();
     phiCellIn = 0.0;
 
-    scalarField& rhoCpPhiCellIn = rhoCpPhiCell.internalField();
+    scalarField& rhoCpPhiCellIn = rhoCpPhiCell.primitiveFieldRef();
     rhoCpPhiCellIn = 0.0;
 
 #   include "mapAirToCell.H"
diff --git a/src/appSrc/muAir.H b/src/appSrc/muAir.H
index f682a5a..35b2340 100644
--- a/src/appSrc/muAir.H
+++ b/src/appSrc/muAir.H
@@ -1,7 +1,7 @@
 // obtain mu(air mixture) from mu(species components)
 
 {
-    scalarField& mu = muAir.internalField();
+    scalarField& mu = muAir.primitiveFieldRef();
     mu = 0;
 
     // initially using a simple linear combination weighted by mass fraction
diff --git a/src/appSrc/muFuel.H b/src/appSrc/muFuel.H
index 50f0562..d9ac631 100644
--- a/src/appSrc/muFuel.H
+++ b/src/appSrc/muFuel.H
@@ -1,7 +1,7 @@
 // obtain mu(fuel mixture) from mu(species components)
 
 {
-    scalarField& mu = muFuel.internalField();
+    scalarField& mu = muFuel.primitiveFieldRef();
     mu = 0;
 
     // initially using a simple linear combination weighted by mass fraction
diff --git a/src/appSrc/nuAir.H b/src/appSrc/nuAir.H
index 472d6e2..733dcea 100644
--- a/src/appSrc/nuAir.H
+++ b/src/appSrc/nuAir.H
@@ -1,6 +1,6 @@
 // obtain nu(air mixture) from mu(air mixture) and rho(air mixture)
 {
-    scalarField& nu = nuAir.internalField();
+    scalarField& nu = nuAir.primitiveFieldRef();
     nu = 0;
     nu = muAir/rhoAir;
     nuAir.correctBoundaryConditions();
diff --git a/src/appSrc/nuFuel.H b/src/appSrc/nuFuel.H
index 78612cb..0e8eb64 100644
--- a/src/appSrc/nuFuel.H
+++ b/src/appSrc/nuFuel.H
@@ -1,6 +1,6 @@
 // obtain nu(fuel mixture) from mu(fuel mixture) and rhoFuel(fuel mixture)
 {
-    scalarField& nu = nuFuel.internalField();
+    scalarField& nu = nuFuel.primitiveFieldRef();
     nu = 0;
     nu = muFuel/rhoFuel;
     nuFuel.correctBoundaryConditions();
diff --git a/src/appSrc/rhoAir.H b/src/appSrc/rhoAir.H
index ef05def..562ab4f 100644
--- a/src/appSrc/rhoAir.H
+++ b/src/appSrc/rhoAir.H
@@ -13,7 +13,7 @@
     {
         // sumAirYIonMI was calculated in getXair.H
 
-        scalarField& rho = rhoAir.internalField();
+        scalarField& rho = rhoAir.primitiveFieldRef();
         rho = pAir/(RR*Tair)/sumAirYIonMI;
         rhoAir.correctBoundaryConditions();
 
diff --git a/src/appSrc/rhoFuel.H b/src/appSrc/rhoFuel.H
index 80f6879..2ffbc28 100644
--- a/src/appSrc/rhoFuel.H
+++ b/src/appSrc/rhoFuel.H
@@ -13,7 +13,7 @@
     {
         // sumFuelYIonMI was calculated in getXfuel.H
 
-        scalarField& rho = rhoFuel.internalField();
+        scalarField& rho = rhoFuel.primitiveFieldRef();
         rho = pFuel/(RR*Tfuel)/sumFuelYIonMI;
         rhoFuel.correctBoundaryConditions();
 
diff --git a/src/appSrc/solveAir.H b/src/appSrc/solveAir.H
index 9de1e0e..7da2b67 100755
--- a/src/appSrc/solveAir.H
+++ b/src/appSrc/solveAir.H
@@ -24,7 +24,7 @@
 
     for (int corr = 0; corr < nCorr; corr++)
     {
-        pAir.boundaryField().updateCoeffs();
+        pAir.boundaryFieldRef().updateCoeffs();
         rUAair = 1.0/UEqn.A();
         Uair = UEqn.H()/UEqn.A();
         Uair.correctBoundaryConditions();
diff --git a/src/appSrc/solveElectrochemistry.H b/src/appSrc/solveElectrochemistry.H
index 52c2943..cf91f89 100755
--- a/src/appSrc/solveElectrochemistry.H
+++ b/src/appSrc/solveElectrochemistry.H
@@ -51,11 +51,11 @@
 
     forAll(Yair, j)
     {
-        sumAYjOnMj += Yair[j].boundaryField()[cathodeID]/airSpecies[j].MW();
+        sumAYjOnMj += Yair[j].boundaryFieldRef()[cathodeID]/airSpecies[j].MW();
     }
 
     PtrList<scalarField> xAir(airSpecies.size());
-    const scalarField& pAirPatch = pAir.boundaryField()[cathodeID]; //added Qing 04.06.2014
+    const scalarField& pAirPatch = pAir.boundaryFieldRef()[cathodeID]; //added Qing 04.06.2014
     forAll(airSpecies, s)
     {
         if (airSpecies[s].rSign() != 0)    //reactant or product (S.Keuler, corrected 30.01.2014)
@@ -63,13 +63,13 @@
             xAir.set
             (
                 s,
-                new scalarField(anodeT.size(), 0)
+                new scalarField(cathodeT.size(), 0)
             );
 
             xAir[s] =
                 cathodeToAnode.faceInterpolate
                 (
-                    Yair[s].boundaryField()[cathodeID]/airSpecies[s].MW()
+                    Yair[s].boundaryFieldRef()[cathodeID]/airSpecies[s].MW()
                     /sumAYjOnMj
                 );
 
@@ -94,13 +94,14 @@
 
     forAll(Yfuel, j)
     {
-      sumFYjOnMj += Yfuel[j].boundaryField()[anodeID]/fuelSpecies[j].MW();
+      sumFYjOnMj += Yfuel[j].boundaryFieldRef()[anodeID]/fuelSpecies[j].MW();
     }
 
     PtrList<scalarField> xFuel(fuelSpecies.size());
-    const scalarField& pFuelPatch = pFuel.boundaryField()[anodeID]; //added hj 03-09-2013
+    const scalarField& pFuelPatch = pFuel.boundaryFieldRef()[anodeID]; //added hj 03-09-2013
     forAll(fuelSpecies,s)
     {
+
         if(fuelSpecies[s].rSign() != 0)    //reactant or product
         {
             xFuel.set
@@ -110,10 +111,10 @@
             );
 
             xFuel[s] =
-            (
-                Yfuel[s].boundaryField()[cathodeID]/fuelSpecies[s].MW()
-                /sumFYjOnMj
-            );
+                (
+                    Yfuel[s].boundaryFieldRef()[anodeID]/fuelSpecies[s].MW()
+                    /sumFYjOnMj
+                );
 
             // ensure positivity
             xFuel[s] = Foam::max(xFuel[s], Foam::doubleScalarSMALL);
@@ -136,7 +137,7 @@
     // Takes voltage from cell center of cell adjacent to boundary face and applies it at the face   
     forAll(electrolyteAnodePatch, facei)
     {
-	voltage.boundaryField()[electrolyteAnodeID][facei] = voltage[electrolyteAnodePatch.faceCells()[facei]];
+	voltage.boundaryFieldRef()[electrolyteAnodeID][facei] = voltage[electrolyteAnodePatch.faceCells()[facei]];
     }
 
 
@@ -234,8 +235,8 @@
                 fixedGradientFvPatchScalarField& YsBC =
                     refCast<fixedGradientFvPatchScalarField>
                     (
-                        //Ys.boundaryField()[anodeID]
-                        Ys.boundaryField()[cathodeID]
+                        //Ys.boundaryFieldRef()[anodeID]
+                        Ys.boundaryFieldRef()[cathodeID]
                     );
 
                 // gradient boundary condition
@@ -266,7 +267,7 @@
 
 
         // Set the interface velocity condition
-        Uair.boundaryField()[cathodeID] ==
+        Uair.boundaryFieldRef()[cathodeID] ==
         (
             -anodeToCathode.faceInterpolate(mfluxSum)
             /rhoAir.boundaryField()[cathodeID]
@@ -327,7 +328,7 @@
                 fixedGradientFvPatchScalarField& YsBC =
                     refCast<fixedGradientFvPatchScalarField>
                     (
-                        Ys.boundaryField()[anodeID]
+                        Ys.boundaryFieldRef()[anodeID]
                     );
 
                 // gradient boundary condition
@@ -352,7 +353,7 @@
 
 
         // Set the interface velocity condition
-        Ufuel.boundaryField()[anodeID] ==
+        Ufuel.boundaryFieldRef()[anodeID] ==
         (
             -(mfluxSum)
             /rhoFuel.boundaryField()[anodeID]
diff --git a/src/appSrc/solveFuel.H b/src/appSrc/solveFuel.H
index ae77ece..ce5981b 100755
--- a/src/appSrc/solveFuel.H
+++ b/src/appSrc/solveFuel.H
@@ -28,7 +28,7 @@
 
     for (int corr = 0; corr < nCorr; corr++)
     {
-        pFuel.boundaryField().updateCoeffs();
+        pFuel.boundaryFieldRef().updateCoeffs();
         rUAfuel = 1.0/UEqn.A();
         Ufuel = UEqn.H()/UEqn.A();
         Ufuel.correctBoundaryConditions();
diff --git a/src/libSrc/Make/options b/src/libSrc/Make/options
index 9a986c6..b81a978 100644
--- a/src/libSrc/Make/options
+++ b/src/libSrc/Make/options
@@ -4,12 +4,13 @@ endif
 
 EXE_INC = -g \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/AABBTree \
     -I$(LIB_SRC)/meshTools/coordinateSystems \
     -I$(LIB_SRC)/meshTools/coordinateSystems/coordinateRotation \
     -I$(LIB_SRC)/meshTools/AMIInterpolation/patches/cyclic/cyclicAMILduInterfaceField \
     -I$(LIB_SRC)/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMILduInterfaceField \
     -I$(LIB_SRC)/meshTools/AMIInterpolation/AMIInterpolation \
-    -I$(LIB_SRC)/meshTools/searchableSurface \
+    -I$(LIB_SRC)/meshTools/searchableSurfaces/searchableSurface \
     -I$(LIB_SRC)/meshTools/AMIInterpolation/faceAreaIntersect \
     -I$(LIB_SRC)/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod \
     -I$(LIB_SRC)/meshTools/indexedOctree \
diff --git a/src/libSrc/MeshWave/FaceCellWave.C b/src/libSrc/MeshWave/FaceCellWave.C
index ba38cf3..e261e91 100644
--- a/src/libSrc/MeshWave/FaceCellWave.C
+++ b/src/libSrc/MeshWave/FaceCellWave.C
@@ -477,7 +477,7 @@ void Foam::FaceCellWave<Type>::sendPatchInfo
     const List<Type>& faceInfo
 ) const
 {
-    OPstream toNeighbour(Pstream::blocking, neighbour);
+    OPstream toNeighbour(Pstream::commsTypes::blocking, neighbour);
 
     writeFaces(nFaces, faceLabels, faceInfo, toNeighbour);
 }
@@ -492,7 +492,7 @@ Foam::label Foam::FaceCellWave<Type>::receivePatchInfo
     List<Type>& faceInfo
 ) const
 {
-    IPstream fromNeighbour(Pstream::blocking, neighbour);
+    IPstream fromNeighbour(Pstream::commsTypes::blocking, neighbour);
 
     label nFaces = 0;
     readFaces(nFaces, faceLabels, faceInfo, fromNeighbour);
diff --git a/src/libSrc/porosityModel/porosityModel/porousZone.C b/src/libSrc/porosityModel/porosityModel/porousZone.C
index 8f6b64a..383452b 100644
--- a/src/libSrc/porosityModel/porosityModel/porousZone.C
+++ b/src/libSrc/porosityModel/porosityModel/porousZone.C
@@ -179,7 +179,7 @@ Foam::tmp<Foam::vectorField> Foam::porousZone::porousZone::force
 
     if (!cellZoneIDs_.empty())
     {
-        this->calcForce(U, rho, mu, tforce());
+        this->calcForce(U, rho, mu, tforce.ref());
     }
 
     return tforce;
diff --git a/src/libSrc/smearPatchToMesh/smearPatchToMeshTemplates.C b/src/libSrc/smearPatchToMesh/smearPatchToMeshTemplates.C
index 2aa3a8a..cd4b7a9 100755
--- a/src/libSrc/smearPatchToMesh/smearPatchToMeshTemplates.C
+++ b/src/libSrc/smearPatchToMesh/smearPatchToMeshTemplates.C
@@ -60,7 +60,7 @@ Foam::tmp<Foam::Field<Type> > Foam::smearPatchToMesh::smear
 
     // Create result
     tmp<Field<Type> > tresult(new Field<Type>(mesh_.nCells()));
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const labelList& addr = addressing();
 
