diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C
index ae1b8c20..9787583a 100755
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C
@@ -88,7 +88,7 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
 	pName_("p"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H
index a7abc6e7..d86a92f2 100755
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H
@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueTubeFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
 		word pName_; //JOKER pressure
 		
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C
index ed490709..e696fe31 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C
@@ -73,7 +73,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H
index 02c230e7..96f4e20c 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H
@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueVoidfractionFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
         word voidfractionFieldName_;
 
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
index e36e9b54..bc29645e 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
@@ -125,15 +125,7 @@ cfdemCloud::cfdemCloud
     ),
     turbulence_
     (
-        #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                mesh.lookupObject<compressible::turbulenceModel>
-            #else
-                mesh.lookupObject<incompressible::turbulenceModel>
-            #endif
-        #elif defined(version15)
-            mesh.lookupObject<incompressible::RASModel>
-        #endif
+        mesh.lookupObject<turbulenceModel>
         (
             turbulenceModelType_
         )
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
index ffb3a589..1eadd20a 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
@@ -49,11 +49,7 @@ SourceFiles
 #include "fvCFD.H"
 #include "IFstream.H"
 
-#if defined(version21) || defined(version16ext)
-    #include <turbulenceModel.H>
-#elif defined(version15)
-    #include <RASModel.H>
-#endif
+#include <turbulenceModel.H>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -157,15 +153,7 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
-    #if defined(version21) || defined(version16ext)
-        #ifdef compre
-            const compressible::turbulenceModel& turbulence_;
-        #else
-            const incompressible::turbulenceModel& turbulence_;
-        #endif
-    #elif defined(version15)
-        const incompressible::RASModel& turbulence_;
-    #endif
+    const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
 
@@ -364,16 +352,8 @@ public:
 
        inline autoPtr<liggghtsCommandModel>* liggghtsCommand() const;
 
-       #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                inline const compressible::turbulenceModel& turbulence() const;
-            #else
-                inline const incompressible::turbulenceModel& turbulence() const;
-            #endif
-       #elif defined(version15)
-           inline const incompressible::RASModel& turbulence() const;
-       #endif
-
+       inline const turbulenceModel& turbulence() const;
+       
     // Write
 
       // write cfdemCloud internal data
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
index 29558e29..84cd29c5 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
@@ -302,15 +302,7 @@ inline autoPtr<liggghtsCommandModel>* cfdemCloud::liggghtsCommand() const
     return liggghtsCommand_;
 }
 
-#if defined(version21) || defined(version16ext)
-    #ifdef compre
-        inline const compressible::turbulenceModel& cfdemCloud::turbulence() const
-    #else
-        inline const incompressible::turbulenceModel& cfdemCloud::turbulence() const
-    #endif
-#elif defined(version15)
-    inline const incompressible::RASModel& cfdemCloud::turbulence() const
-#endif
+inline const turbulenceModel& cfdemCloud::turbulence() const
 {
     return turbulence_;
 }
diff --git a/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C b/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C
index b186db91..0f65f783 100644
--- a/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C
+++ b/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C
@@ -59,7 +59,7 @@ void averagingModel::undoVectorAverage
 {
 // WARNING - not sure if this is valid for dilute model!!!
 
-    if(!single) fieldPrev.internalField() = fieldNext.internalField();
+    if(!single) fieldPrev.ref() = fieldNext.ref();
 
     label cellI;
     vector valueVec;
@@ -303,13 +303,13 @@ void averagingModel::setDSauter
 
 void averagingModel::resetVectorAverage(volVectorField& prev,volVectorField& next,bool single) const
 {
-    if(!single) prev.internalField() = next.internalField();
-    next.internalField() = vector::zero;
+    if(!single) prev.ref() = next.ref();
+    next.primitiveFieldRef() = vector::zero;
 }
 
 void averagingModel::resetWeightFields() const
 {
-    UsWeightField_.internalField() = 0;
+    UsWeightField_.ref() = 0;
 }
 
 
@@ -352,12 +352,12 @@ tmp<volVectorField> Foam::averagingModel::UsInterp() const
 
     if (particleCloud_.dataExchangeM().couplingStep() > 1)
     {
-        tsource() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
+        tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
                     + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
     }
     else
     {
-        tsource() = UsNext_;
+        tsource.ref() = UsNext_;
     }
 
     return tsource;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C
index f0becad6..50f5a7e4 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C
@@ -64,11 +64,7 @@ Archimedes::Archimedes
     propsDict_(dict.subDict(typeName + "Props")),
     twoDimensional_(false),
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
 
     //Append the field names to be probed
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.H
index 41d3c51d..edb09ff4 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.H
@@ -65,11 +65,7 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
+    const uniformDimensionedVectorField& g_;    // ref to gravity
 
 
 public:
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C
index 846f43aa..46c016b6 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C
@@ -66,11 +66,7 @@ ArchimedesIB::ArchimedesIB
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), //mod by alice
     voidfractions_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),//mod by alice
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif  defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, "archimedesIBF.logDat");
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.H b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.H
index 7b1ff8e6..c80b3900 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.H
@@ -71,12 +71,8 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
-
+    const uniformDimensionedVectorField& g_;    // ref to gravity
+    
 
 public:
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
index 7a94f570..09cb8c46 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
@@ -130,7 +130,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     allocateMyArrays();
 
     // reset Scalar field
-    EuField.internalField() = 0.0;
+    EuField.primitiveFieldRef() = 0.0;
 
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
@@ -232,7 +232,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     );
 
     // scale with -1/(Vcell*rho*Cp)
-    EuField.internalField() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
+    EuField.primitiveFieldRef() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
 
     // limit source term
     forAll(EuField,cellI)
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/particleCellVolume/particleCellVolume.C b/src/lagrangian/cfdemParticle/subModels/forceModel/particleCellVolume/particleCellVolume.C
index ee406c6e..54c956b1 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/particleCellVolume/particleCellVolume.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/particleCellVolume/particleCellVolume.C
@@ -120,7 +120,7 @@ void particleCellVolume::setForce() const
     {
         if(verbose_) Info << "particleCellVolume.C - setForce()" << endl;
 
-        scalarField_.internalField()=0.;
+        scalarField_.ref()=0.;
 
         // get reference to actual field
         const volScalarField& field = mesh_.lookupObject<volScalarField>(scalarFieldName_);
@@ -143,8 +143,8 @@ void particleCellVolume::setForce() const
                 scalarField2_[cellI] = 0.;
             }
         }
-        scalarField_.internalField() = gSum(scalarField_);
-        scalarField2_.internalField() = gSum(scalarField2_);
+        scalarField_.ref() = gSum(scalarField_);
+        scalarField2_.ref() = gSum(scalarField2_);
 
         if(verbose_)
         {
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/volWeightedAverage/volWeightedAverage.C b/src/lagrangian/cfdemParticle/subModels/forceModel/volWeightedAverage/volWeightedAverage.C
index e4054c75..d4d3d603 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/volWeightedAverage/volWeightedAverage.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/volWeightedAverage/volWeightedAverage.C
@@ -179,7 +179,7 @@ void volWeightedAverage::setForce() const
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             integralValue = gSum(scalarFields_[i]);
             volWeightedAverage = integralValue / (totVol_all+SMALL);
-            scalarFields_[i].internalField() = volWeightedAverage;
+            scalarFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {
@@ -224,7 +224,7 @@ void volWeightedAverage::setForce() const
 
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             volWeightedAverage = gSum(vectorFields_[i]) / (totVol_all+SMALL);
-            vectorFields_[i].internalField() = volWeightedAverage;
+            vectorFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {
diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C
index 224cd364..861ca33e 100644
--- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C
+++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C
@@ -61,14 +61,8 @@ engineSearch::engineSearch
 :
     locateModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    //faceDecomp_(propsDict_.lookup("faceDecomp")),
     treeSearch_(propsDict_.lookup("treeSearch")),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
-    //searchEngine_(particleCloud_.mesh(),faceDecomp_) // only 2.0.x
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H
index e855b438..33070599 100644
--- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H
+++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H
@@ -61,8 +61,6 @@ private:
 
     dictionary propsDict_;
 
-    //Switch faceDecomp_;
-
     Switch treeSearch_;
 
 protected:
diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C b/src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C
index 1ee8b044..afbdf504 100644
--- a/src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C
+++ b/src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C
@@ -89,11 +89,7 @@ label standardSearch::findCell
             for(int i=0;i<3;i++) position[i] = positions[index][i];
 
             // find cell
-            #ifdef version16ext
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position);
-            #elif defined(version21)
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-            #endif
+            cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
         }
         else cellIDs[index][0] = -1;
     }
@@ -108,11 +104,7 @@ label standardSearch::findSingleCell
 ) const
 {
     // find cell
-    #ifdef version16ext
-        return particleCloud_.mesh().findCell(position);
-    #elif defined(version21)
-        return particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-    #endif
+    return particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
 }
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C b/src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C
index 9975c0af..cc8137f4 100644
--- a/src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C
+++ b/src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C
@@ -64,11 +64,7 @@ turboEngineSearch::turboEngineSearch
     propsDict_(dict.subDict(typeName + "Props")),
     treeSearch_(propsDict_.lookup("treeSearch")),
     bb_(particleCloud_.mesh().points(),false),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C
index fc4572ed..8afed35b 100644
--- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C
+++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C
@@ -145,18 +145,18 @@ tmp<volVectorField> explicitCouple::expMomSource() const
                 if (fNext_[cellI][i] > fLimit_[i]) fNext_[cellI][i] = fLimit_[i];
             }
         }
-        tsource() = fPrev_;
+        tsource.ref() = fPrev_;
     }else
     {
-        tsource() = (1 - tsf) * fPrev_ + tsf * fNext_;
+        tsource.ref() = (1 - tsf) * fPrev_ + tsf * fNext_;
     }
     return tsource;
 }
 
 void Foam::explicitCouple::resetMomSourceField() const
 {
-    fPrev_.internalField() = fNext_.internalField();
-    fNext_.internalField() = vector::zero;
+    fPrev_.ref() = fNext_.ref();
+    fNext_.primitiveFieldRef() = vector::zero;
 }
 
 inline vector Foam::explicitCouple::arrayToField(label cellI) const
diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C
index 88b9631f..377fc137 100644
--- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C
+++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C
@@ -166,10 +166,10 @@ tmp<volScalarField> implicitCouple::impMomSource() const
             // limiter
             if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_;
         }
-        tsource() = KslPrev_;
+        tsource.ref() = KslPrev_;
     }else
     {
-        tsource() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
+        tsource.ref() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
     }
 
     return tsource;
@@ -177,8 +177,8 @@ tmp<volScalarField> implicitCouple::impMomSource() const
 
 void Foam::implicitCouple::resetMomSourceField() const
 {
-    KslPrev_.internalField() = KslNext_.internalField();
-    KslNext_.internalField() = 0;
+    KslPrev_.ref() = KslNext_.ref();
+    KslNext_.primitiveFieldRef() = 0;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
index 695e5154..1af57abe 100644
--- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
+++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
@@ -100,7 +100,7 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
     volScalarField sSmoothField = sSmoothField_;
 
     sSmoothField.dimensions().reset(fieldSrc.dimensions());
-    sSmoothField.internalField()=fieldSrc.internalField();
+    sSmoothField.ref()=fieldSrc.internalField();
     sSmoothField.correctBoundaryConditions();
     sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     sSmoothField.oldTime()=fieldSrc;
@@ -141,7 +141,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -176,7 +176,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -210,7 +210,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     forAll(vSmoothField,cellI)
     {
         if ( mag(vSmoothField.oldTime().internalField()[cellI]) > 0.0f)  // have a vector in the OLD vSmoothField, so keep it!
-            NLarge()[cellI] = sourceStrength;
+            NLarge.ref()[cellI] = sourceStrength;
     }
 
     // do the smoothing
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
index 1c90f8c8..9a929d67 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
@@ -91,7 +91,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
index 27bc9423..2b60cfac 100755
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
@@ -96,7 +96,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
index b8d5bbad..179a3b6a 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
@@ -90,7 +90,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
index 54182eb7..b5ddbddf 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
@@ -127,19 +127,19 @@ tmp<volScalarField> Foam::voidFractionModel::voidFractionInterp() const
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
     if(1-tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
     {
-        tsource() = voidfractionPrev_;
+        tsource.ref() = voidfractionPrev_;
     }
     else
     {
-        tsource() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
+        tsource.ref() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
     }
     return tsource;
 }
 
 void Foam::voidFractionModel::resetVoidFractions() const
 {
-    voidfractionPrev_.internalField() = voidfractionNext_.internalField();
-    voidfractionNext_.internalField() = 1;
+    voidfractionPrev_.ref() = voidfractionNext_.ref();
+    voidfractionNext_.ref() = 1;
 }
 
 /*void Foam::voidFractionModel::undoVoidFractions(double**const& mask) const