diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
index df719b60..36db2d47 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C
@@ -160,6 +160,7 @@ cfdemCloud::cfdemCloud
             turbulenceModelType_
         )
     ),
+    particlePropertyTable(32),
     dataExchangeModel_
     (
         dataExchangeModel::New
@@ -402,6 +403,22 @@ cfdemCloud::~cfdemCloud()
     if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
     if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
     if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
+
+    for
+    (
+        HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+        iter != particlePropertyTable.end();
+        ++iter
+    )
+    {
+        if ((*(iter().ti)) == typeid(int**)) {
+            dataExchangeM().destroy(iter().ref<int**>(),-1);
+        } else if ((*(iter().ti)) == typeid(double**)) {
+            dataExchangeM().destroy(iter().ref<double**>(),-1);
+        } else {
+            FatalError << "Trying to destroy property of type " << iter().ti->name() << endl;
+        }
+    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -780,9 +797,60 @@ bool cfdemCloud::reAllocArrays()
         if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
         if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
         if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
+
+        for
+        (
+            HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+            iter != particlePropertyTable.end();
+            ++iter
+        )
+        {
+            if (iter().size > 0) {
+                ///Info << "!! about to realloc property of type " << iter().ti->name() << endl;
+                if ((*(iter().ti)) == typeid(int**)) {
+                    dataExchangeM().allocateArray(iter().ref<int**>(),iter().initVal,iter().size);
+                } else if ((*(iter().ti)) == typeid(double**)) {
+                    dataExchangeM().allocateArray(iter().ref<double**>(),iter().initVal,iter().size);
+                } else {
+                    FatalError << "Trying to realloc property of type " << iter().ti->name() << endl;
+                }
+            }
+        }
+
         arraysReallocated_ = true;
         return true;
     }
+    else
+    {
+        for
+        (
+            HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+            iter != particlePropertyTable.end();
+            ++iter
+        )
+        {
+            if (iter().size > 0 && iter().reset) {
+                if ((*(iter().ti)) == typeid(int**)) {
+                    int**& property = iter().ref<int**>();
+                    for (int index=0; index<numberOfParticles(); ++index) {
+                        for (int icomponent=0; icomponent<iter().size; ++icomponent) {
+                            property[index][icomponent] = iter().initVal;
+                        }
+                    }
+                } else if ((*(iter().ti)) == typeid(double**)) {
+                    double**& property = iter().ref<double**>();
+                    for (int index=0; index<numberOfParticles(); ++index) {
+                        for (int icomponent=0; icomponent<iter().size; ++icomponent) {
+                            property[index][icomponent] = iter().initVal;
+                        }
+                    }
+                } else {
+                    FatalError << "Trying to reset property of type " << iter().ti->name() << endl;
+                }
+            }
+        }
+    }
+
     return false;
 }
 
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
index 0b0d02ee..636c2c18 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H
@@ -45,6 +45,7 @@ SourceFiles
 // choose version
 #include "OFversion.H"
 #include <vector>
+#include <typeinfo>
 
 #include "fvCFD.H"
 #include "IFstream.H"
@@ -188,6 +189,21 @@ protected:
 
     const turbulenceModel& turbulence_;
 
+    struct particleProperty {
+        void** property;            // pointer to per-particle data; memory deallocation by cfdemCloud
+        const std::type_info* ti;   // type of particle data (int**, double**)
+        int size;                   // size of single particle data; memory allocation by cfdemCloud if size > 0
+        double initVal;             // initial property value
+        bool reset;                 // if true, data is reset to initial value every coupling step
+        template<typename T>
+        T& ref() {
+          if (*ti == typeid(T)) return *reinterpret_cast<T*>(&property);
+          else throw std::bad_cast();
+        }
+    };
+
+    HashTable<particleProperty> particlePropertyTable; // table of registered per-particle properties
+
     autoPtr<dataExchangeModel> dataExchangeModel_;
 
     PtrList<forceModel> forceModel_;
@@ -443,6 +459,14 @@ public:
         void otherForces(volVectorField&);
 
         bool checkPeriodicCells() const { return checkPeriodicCells_; }
+
+        template<typename T>
+        void registerParticleProperty(const word& property, int size=0, double initVal=0.0, bool reset=true);
+        template<typename T>
+        T& getParticlePropertyRef(const word& property);
+protected:
+        virtual int**& getParticlePropertyImpl(const word& property, int**);
+        virtual double**& getParticlePropertyImpl(const word& property, double**);
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
index 5c9d426f..dc3a4125 100644
--- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
+++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H
@@ -403,5 +403,28 @@ inline const turbulenceModel& cfdemCloud::turbulence() const
     return turbulence_;
 }
 
+
+template<typename T>
+void cfdemCloud::registerParticleProperty(const word& property, int size, double initVal, bool reset)
+{
+    particlePropertyTable.insert(property,{NULL,&typeid(T),size,initVal,reset});
+}
+
+template<typename T>
+T& cfdemCloud::getParticlePropertyRef(const word& property)
+{
+    return getParticlePropertyImpl(property, static_cast<T>(0));
+}
+
+inline int**& cfdemCloud::getParticlePropertyImpl(const word& property, int**)
+{
+    return particlePropertyTable[property].ref<int**>();
+}
+
+inline double**& cfdemCloud::getParticlePropertyImpl(const word& property, double**)
+{
+    return particlePropertyTable[property].ref<double**>();
+}
+
 }
 // ************************************************************************* //
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
index 9ce61c26..3ad1a38b 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
@@ -71,15 +71,17 @@ diffusionCoefficient::diffusionCoefficient
     pressureFieldName_(propsDict_.lookupOrDefault<word>("pressureFieldName","p")),
     P_(sm.mesh().lookupObject<volScalarField>(pressureFieldName_)),
     partPressureName_(propsDict_.lookupOrDefault<word>("partPressureName","partP")),
-    partPressure_(NULL),
     X_(speciesNames_.size()),
     diffusantGasNames_(propsDict_.lookup("diffusantGasNames")),
-    diffusionCoefficients_(diffusantGasNames_.size(),NULL),
     Xdiffusant_(diffusantGasNames_.size()),
     initialized_(false)
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partPressureName_,1);
+    for (int i=0; i<diffusantGasNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>(diffusantGasNames_[i],1);
+    }
     createCoeffs();
     molWeightTable();
 }
@@ -88,37 +90,10 @@ diffusionCoefficient::diffusionCoefficient
 
 diffusionCoefficient::~diffusionCoefficient()
 {
-    particleCloud_.dataExchangeM().destroy(partPressure_,1);
-    for (int i=0; i<diffusantGasNames_.size(); i++) particleCloud_.dataExchangeM().destroy(diffusionCoefficients_[i],1);
-
-    coeffs.clearStorage();
-    molWeight.clearStorage();
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void diffusionCoefficient::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
-        for (int i=0; i<diffusantGasNames_.size(); i++)
-        {
-            particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1,"nparticles");
-        }
-    }
-}
-
-void diffusionCoefficient::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
-    for (int i=0; i<diffusantGasNames_.size(); i++)
-    {
-        particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
-    }
-}
 
 void diffusionCoefficient::init()
 {
@@ -156,9 +131,6 @@ void diffusionCoefficient::execute()
         init();
     }
 
-    // realloc the arrays
-    reAllocMyArrays();
-
     label  cellI=0;
     scalar Tfluid(0);
     scalar Pfluid(0);
@@ -172,6 +144,8 @@ void diffusionCoefficient::execute()
     interpolationCellPoint <scalar> TInterpolator_(tempField_);
     interpolationCellPoint <scalar> PInterpolator_(P_);
 
+    double**& partPressure_ = particleCloud_.getParticlePropertyRef<double**>(partPressureName_);
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
@@ -207,6 +181,7 @@ void diffusionCoefficient::execute()
 
             for (int j=0; j<diffusantGasNames_.size(); j++)
             {
+                double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
                 TotalFraction_[j] = 0.0;
                 dBinary_ = 0.0;
 
@@ -243,9 +218,9 @@ void diffusionCoefficient::execute()
 
                             // pass on dCoeff values to array
                             if (TotalFraction_[j] < VSMALL)
-                                diffusionCoefficients_[j][index][0] = VSMALL;
+                                diffusionCoefficients_[index][0] = VSMALL;
                             else
-                                diffusionCoefficients_[j][index][0] = (1.0 - Xdiffusant_[j][cellI]) / TotalFraction_[j];
+                                diffusionCoefficients_[index][0] = (1.0 - Xdiffusant_[j][cellI]) / TotalFraction_[j];
                         }
                         else
                         {
@@ -258,7 +233,7 @@ void diffusionCoefficient::execute()
                 }
 
                 if(verbose_)
-                    Info << "diffusionCoefficient of species " << diffusantGasNames_[j] << " = " << diffusionCoefficients_[j][index][0] << endl;
+                    Info << "diffusionCoefficient of species " << diffusantGasNames_[j] << " = " << diffusionCoefficients_[index][0] << endl;
             }
         }
     }
@@ -268,7 +243,8 @@ void diffusionCoefficient::execute()
     for (int j=0; j<diffusantGasNames_.size(); j++)
     {
         word pushName = diffusantGasNames_[j] + "_diffCoeff";
-        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_[j]);
+        double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
+        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_);
     }
 
     Info << "give data done" << endl;
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
index c8e9b1cf..72ad8a68 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
@@ -76,19 +76,15 @@ private:
 
     word partPressureName_;
 
-    mutable double **partPressure_;
-
     UPtrList<volScalarField> X_;
 
     wordList diffusantGasNames_;
 
-    mutable List<double**> diffusionCoefficients_;
-
     UPtrList<volScalarField> Xdiffusant_;
 
-    HashTable<scalar, word> coeffs;
+    HashTable<scalar> coeffs;
 
-    HashTable<scalar, word> molWeight;
+    HashTable<scalar> molWeight;
 
     void createCoeffs();
 
@@ -100,8 +96,6 @@ private:
     // calculate denominator part diffusion volume equation
     double calcDiffVol(int, int);
 
-    void allocateMyArrays() const;
-
     bool initialized_;
 
     void init();
@@ -126,8 +120,6 @@ private:
 
     // Member Functions
         void execute();
-
-        void reAllocMyArrays() const;
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
index 2f49a729..2293e7bf 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
@@ -63,51 +63,26 @@ massTransferCoeff::massTransferCoeff
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partNuName_(propsDict_.lookupOrDefault<word>("partViscos","partNu")),
-    partNu_(NULL),
-    partReynolds_(propsDict_.lookupOrDefault<word>("partReynolds","partRe")),
-    partRe_(NULL),
+    partReName_(propsDict_.lookupOrDefault<word>("partReynolds","partRe")),
     scaleDia_(1)
 {
     particleCloud_.checkCG(true);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partNuName_,1);
+    particleCloud_.registerParticleProperty<double**>(partReName_,1);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 massTransferCoeff::~massTransferCoeff()
 {
-    int nP_ = particleCloud_.numberOfParticles();
-
-    particleCloud_.dataExchangeM().destroy(partNu_,nP_);
-    particleCloud_.dataExchangeM().destroy(partRe_,nP_);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void massTransferCoeff::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1,"nparticles");
-    }
-}
-
-void massTransferCoeff::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void massTransferCoeff::execute()
 {
-    // realloc the arrays
-    reAllocMyArrays();
-
     #ifdef compre
         const volScalarField nufField = particleCloud_.turbulence().mu()/rho_;
     #else
@@ -139,6 +114,9 @@ void massTransferCoeff::execute()
     interpolationCellPoint <vector> UluidInterpolator_(U_);
     interpolationCellPoint <scalar> voidfractionInterpolator_(voidfraction_);
 
+    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuName_);
+    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReName_);
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI=particleCloud_.cellIDs()[index][0];
@@ -195,7 +173,7 @@ void massTransferCoeff::execute()
 
     // give DEM data
     particleCloud_.dataExchangeM().giveData(partNuName_, "scalar-atom", partNu_);
-    particleCloud_.dataExchangeM().giveData(partReynolds_, "scalar-atom", partRe_);
+    particleCloud_.dataExchangeM().giveData(partReName_, "scalar-atom", partRe_);
 
     Info << "give data done" << endl;
 }
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
index 5021fbc8..927cac8b 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
@@ -56,27 +56,21 @@ private:
 
     const fvMesh& mesh_;
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partNuName_;
+    const word partNuName_;
 
-    mutable double **partNu_;
-
-    word partReynolds_;
-
-    mutable double **partRe_;
-
-    void allocateMyArrays() const;
+    const word partReName_;
 
     mutable scalar scaleDia_;
 
@@ -101,8 +95,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
index 0d2151c7..99673ede 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
@@ -57,7 +57,7 @@ reactantPerParticle::reactantPerParticle
     propsDict_(dict.subDict(typeName + "Props")),
     mesh_(sm.mesh()),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
-    reactantPerParticle_(NULL),
+    partReactantName_("reactantPerParticle"),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     particlesPerCell_
@@ -76,32 +76,16 @@ reactantPerParticle::reactantPerParticle
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1))
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partReactantName_,1);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 reactantPerParticle::~reactantPerParticle()
 {
-    particleCloud_.dataExchangeM().destroy(reactantPerParticle_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void reactantPerParticle::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1,"nparticles");
-    }
-}
-
-void reactantPerParticle::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
@@ -112,8 +96,6 @@ void reactantPerParticle::execute()
     {
         return;
     }
-    // realloc the arrays
-    reAllocMyArrays();
 
     particlesPerCell_ *= 0.0;
 
@@ -121,6 +103,7 @@ void reactantPerParticle::execute()
     scalar voidfraction(1);
     scalar cellvolume(0.0);
     scalar particlesPerCell(1.0);
+    double**& reactantPerParticle_ = particleCloud_.getParticlePropertyRef<double**>(partReactantName_);
 
     // first create particles per cell field
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
@@ -147,10 +130,10 @@ void reactantPerParticle::execute()
         if (verbose_) Info << "reactantPerParticle_" << reactantPerParticle_[index][0] << endl;
     }
 
-        // give DEM data
-        particleCloud_.dataExchangeM().giveData("reactantPerParticle", "scalar-atom", reactantPerParticle_);
+    // give DEM data
+    particleCloud_.dataExchangeM().giveData(partReactantName_, "scalar-atom", reactantPerParticle_);
 
-        Info << "give data done" << endl;
+    Info << "give data done" << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H
index 2e369755..a0ddbee0 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H
@@ -56,16 +56,14 @@ private:
 
     bool verbose_;
 
-    mutable double **reactantPerParticle_;
+    const word partReactantName_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
     mutable volScalarField particlesPerCell_;
 
-    void allocateMyArrays() const;
-
     label loopCounter_;
 
     label Nevery_;
@@ -91,8 +89,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
index c485a848..f9719eab 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
@@ -70,8 +70,6 @@ species::species
     speciesNames_(specDict_.lookup("species")),
     mod_spec_names_(speciesNames_.size()),
     X_(speciesNames_.size()),                           //volumeScalarFields of molarFractions
-    molarFractions_(speciesNames_.size(),NULL),         //the value of molar fractions for every species
-    changeOfSpeciesMass_(speciesNames_.size(),NULL),    //the values that are received from DEM with the name of Modified_+species name
     changeOfSpeciesMassFields_(speciesNames_.size()),   //the scalar fields generated with the values from Modified_+species names
     changeOfGasMassField_                               //the total change of Gas Mass field (when the Modified species
     (
@@ -89,18 +87,15 @@ species::species
     tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
     tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
     partTempName_(propsDict_.lookupOrDefault<word>("partTempName","partTemp")),
-    partTemp_(NULL),
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partRhoName_(propsDict_.lookupOrDefault<word>("partRhoName","partRho")),
-    partRho_(NULL),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     // total mole field
     molarConcFieldName_(propsDict_.lookupOrDefault<word>("totalMoleFieldName","molarConc")),
     molarConc_(sm.mesh().lookupObject<volScalarField>(molarConcFieldName_)),
     partMolarConcName_(propsDict_.lookupOrDefault<word>("partMoleName","partMolarConc")),
-    partMolarConc_(NULL),
     loopCounter_(-1),
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
     massSourceCurr_(0.0),
@@ -108,56 +103,25 @@ species::species
     initialized_(false)
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partRhoName_,1);
+    particleCloud_.registerParticleProperty<double**>(partMolarConcName_,1);
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>("X_"+speciesNames_[i],1);
+        particleCloud_.registerParticleProperty<double**>("Modified_"+speciesNames_[i],1);
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 species::~species()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partRho_,1);
-    particleCloud_.dataExchangeM().destroy(partMolarConc_,1);
-
-
-    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(molarFractions_[i],1);
-    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(changeOfSpeciesMass_[i],1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void species::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1,"nparticles");
-
-        for (int i=0; i<speciesNames_.size(); i++)
-        {
-            particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1,"nparticles");
-            particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1,"nparticles");
-        }
-    }
-}
-
-void species::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1);
-
-    for (int i=0; i<speciesNames_.size(); i++)
-    {
-        particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
-    }
-}
-
 void species::init()
 {
     if(verbose_)
@@ -219,8 +183,6 @@ void species::execute()
     {
         return;
     }
-    // realloc the arrays
-    reAllocMyArrays();
 
     // get X_i, T, rho at particle positions
     label  cellI = 0;
@@ -235,6 +197,9 @@ void species::execute()
     interpolationCellPoint <scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint <scalar> molarConcInterpolator_(molarConc_);
 
+    double**& partRho_ = particleCloud_.getParticlePropertyRef<double**>(partRhoName_);
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partMolarConc_ = particleCloud_.getParticlePropertyRef<double**>(partMolarConcName_);
 
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
@@ -264,7 +229,8 @@ void species::execute()
             for (int i=0; i<speciesNames_.size();i++)
             {
                 // attention for indices when not communicating all species
-                molarFractions_[i][index][0] = X_[i][cellI];
+                double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
+                molarFractions_[index][0] = X_[i][cellI];
             }
         }
     }
@@ -274,8 +240,9 @@ void species::execute()
     {
         for(int i =0; i<speciesNames_.size(); i++)
         {
+            double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
             Info << "X_i = " << X_[i].name() << endl;
-            Info << "molarFractions_= " << molarFractions_[i][0][0] << endl;
+            Info << "molarFractions_= " << molarFractions_[0][0] << endl;
             Info << "partRho_[index][0] = " << partRho_[0][0] << endl;
             Info << "rhofluid = " << rhofluid << endl;
             Info << "partTemp_[index][0] = " << partTemp_[0][0] << endl;
@@ -292,7 +259,8 @@ void species::execute()
 
         for (int i=0; i<speciesNames_.size();i++)
         {
-            particleCloud_.dataExchangeM().giveData("X_"+speciesNames_[i],"scalar-atom",molarFractions_[i]);
+            double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
+            particleCloud_.dataExchangeM().giveData("X_"+speciesNames_[i],"scalar-atom",molarFractions_);
         }
 
         if (verbose_) Info << "give data done" << endl;
@@ -305,17 +273,18 @@ void species::execute()
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
         {
+            double**& changeOfSpeciesMass_ = particleCloud_.getParticlePropertyRef<double**>("Modified_"+speciesNames_[i]);
             changeOfSpeciesMassFields_[i].primitiveFieldRef() = 0.0;
             changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
 
-            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i],particleCloud_.dataExchangeM().couplingInterval());
+            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_,particleCloud_.dataExchangeM().couplingInterval());
 
-            if (verbose_) Info << "changeOfSpeciesMass received from DEM = " << changeOfSpeciesMass_[i][0][0] << endl;
+            if (verbose_) Info << "changeOfSpeciesMass received from DEM = " << changeOfSpeciesMass_[0][0] << endl;
 
             particleCloud_.averagingM().setScalarSumCentre
             (
                 changeOfSpeciesMassFields_[i],
-                changeOfSpeciesMass_[i],
+                changeOfSpeciesMass_,
                 particleCloud_.particleWeights(),
                 NULL
             );
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
index 2b9529a0..53763a17 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
@@ -68,44 +68,32 @@ private:
 
     UPtrList<volScalarField> X_;
 
-    mutable List<double**> molarFractions_;
-
-    mutable List<double**> changeOfSpeciesMass_;
-
     PtrList<volScalarField> changeOfSpeciesMassFields_;
 
     volScalarField changeOfGasMassField_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;              // ref to gas temperature field
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;                    // gas temperature at particle positions
-
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partRhoName_;
+    const word partRhoName_;
 
-    mutable double **partRho_;                     // gas density at particle positions
-
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    void allocateMyArrays() const;
-
     // total mole field
-    word molarConcFieldName_;
+    const word molarConcFieldName_;
 
     const volScalarField& molarConc_;
 
-    word partMolarConcName_;
-
-    mutable double **partMolarConc_;
+    const word partMolarConcName_;
 
     label loopCounter_;
 
@@ -141,8 +129,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
         tmp <volScalarField> Smi(const label i) const;
 
         tmp <volScalarField> Sm() const;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
index cf87a4d7..3f11ca68 100755
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
@@ -129,31 +129,31 @@ public:
         template <typename T>
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             T ** const& field
         ) const { getData(name,type,field,couplingStep_-1); }
 
         virtual void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const = 0;
 
         virtual void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const = 0;
 
         virtual void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype="double"
         ) const = 0;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/noDataExchange/noDataExchange.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/noDataExchange/noDataExchange.H
index 1fda4c54..6d7a8e1f 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/noDataExchange/noDataExchange.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/noDataExchange/noDataExchange.H
@@ -78,24 +78,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const {}
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const {}
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const {}
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.C
index 327519ec..285b80db 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.C
@@ -84,8 +84,8 @@ oneWayVTK::~oneWayVTK()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 void oneWayVTK::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label step
 ) const
@@ -194,8 +194,8 @@ void oneWayVTK::getData
 
 void oneWayVTK::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.H
index 72731202..0f030c3f 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.H
@@ -87,16 +87,16 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const
@@ -104,8 +104,8 @@ public:
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.C
index 73aed9d6..fcc4e0bd 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.C
@@ -80,7 +80,7 @@ twoWayFiles::~twoWayFiles()
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-fileName twoWayFiles::getFilePath(word& name, bool in) const
+fileName twoWayFiles::getFilePath(const word& name, bool in) const
 {
     char timeStep[40];
 
@@ -100,7 +100,7 @@ fileName twoWayFiles::getFilePath(word& name, bool in) const
     return particleFilePathOld;
 }
 
-void twoWayFiles::renameFilePath(fileName& particleFilePathOld,word& name) const
+void twoWayFiles::renameFilePath(const fileName& particleFilePathOld, const word& name) const
 {
     char timeStep[40];
 
@@ -116,8 +116,8 @@ void twoWayFiles::renameFilePath(fileName& particleFilePathOld,word& name) const
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayFiles::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label step
 ) const
@@ -164,8 +164,8 @@ void twoWayFiles::getData
 
 void twoWayFiles::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.H
index e5f15814..6279ca5e 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.H
@@ -65,9 +65,9 @@ private:
 
   // private member functions
 
-    fileName getFilePath(word&, bool) const;
+    fileName getFilePath(const word&, bool) const;
 
-    void renameFilePath(fileName&,word&) const;
+    void renameFilePath(const fileName&, const word&) const;
 
 public:
 
@@ -92,24 +92,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const { FatalError<<"function not implemented !!! twoWayFiles::getData!!!\n" << abort(FatalError); }
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
index 45af8a5c..2b028010 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
@@ -138,8 +138,8 @@ double twoWayMPI::DEMVariableValue(word variablename)
 
 void twoWayMPI::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -149,8 +149,8 @@ void twoWayMPI::getData
 
 void twoWayMPI::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -160,8 +160,8 @@ void twoWayMPI::getData
 
 void twoWayMPI::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
index 9801ec52..1c647436 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
@@ -107,24 +107,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C
index da2df875..c1a7ef9d 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C
@@ -178,8 +178,8 @@ twoWayMany2Many::~twoWayMany2Many()
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayMany2Many::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -231,8 +231,8 @@ void twoWayMany2Many::getData
 
 void twoWayMany2Many::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -242,8 +242,8 @@ void twoWayMany2Many::getData
 
 void twoWayMany2Many::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* /*datatype*/
 ) const
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H
index a9697249..18b4a448 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H
@@ -157,24 +157,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
index b8c314f6..43a1c236 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
@@ -253,8 +253,8 @@ twoWayOne2One::~twoWayOne2One()
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayOne2One::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -361,8 +361,8 @@ void twoWayOne2One::getData
 
 void twoWayOne2One::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -383,8 +383,8 @@ void twoWayOne2One::getData
 
 void twoWayOne2One::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H
index a6ee34da..59f5f9f9 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.H
@@ -152,24 +152,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C
index 135850f7..cbd168d3 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C
@@ -49,7 +49,8 @@ heatTransferGranConduction::heatTransferGranConduction
     totalHeatFlux_(0.0),
     QPartPartName_(propsDict_.lookupOrDefault<word>("QPartPartName","QPartPart")),
     QPartPart_
-    (   IOobject
+    (
+        IOobject
         (
             QPartPartName_,
             sm.mesh().time().timeName(),
@@ -72,7 +73,7 @@ heatTransferGranConduction::heatTransferGranConduction
         ),
         sm.mesh(),
         dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     partThermCondField_
     (
@@ -86,7 +87,7 @@ heatTransferGranConduction::heatTransferGranConduction
         ),
         sm.mesh(),
         dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     partTempField_(sm.mesh().lookupObject<volScalarField>("partTemp")),
     prescribedVoidfractionFieldName_(propsDict_.lookupOrDefault<word>("prescribedVoidfractionFieldName","voidfraction")),
@@ -94,11 +95,11 @@ heatTransferGranConduction::heatTransferGranConduction
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","conductiveHeatFlux")),
-    partHeatFlux_(NULL),
     typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
-    partThermCond_(NULL)
+    partThermCondRegName_(typeName + "partThermCond")
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    particleCloud_.registerParticleProperty<double**>(partThermCondRegName_,1);
 
     if (typePartThermCond_[0] < 0.0)
     {
@@ -127,25 +128,13 @@ heatTransferGranConduction::heatTransferGranConduction
 
 heatTransferGranConduction::~heatTransferGranConduction()
 {
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partThermCond_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferGranConduction::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partThermCond_,initVal,1);
-}
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferGranConduction::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
-
     calcPartEffThermCond();
 
     QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField_);
@@ -156,6 +145,7 @@ void heatTransferGranConduction::calcEnergyContribution()
     scalar voidfraction(1);
 
     totalHeatFlux_ = 0.0;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
@@ -199,6 +189,8 @@ void heatTransferGranConduction::calcPartThermCond()
 {
     label cellI=0;
     label partType = 1;
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>(partThermCondRegName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -227,7 +219,8 @@ void heatTransferGranConduction::calcPartThermCond()
 
 void heatTransferGranConduction::heatFlux(label index, scalar vol, scalar voidfraction, scalar QPartPart)
 {
-        partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+    partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
 }
 
 void heatTransferGranConduction::giveData()
@@ -238,6 +231,7 @@ void heatTransferGranConduction::giveData()
         Info << "total conductive particle-particle heat flux [W] (Eulerian) = " << totalHeatFlux_ << endl;
     }
 
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
index 6ab1e180..c98256f9 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
@@ -54,7 +54,7 @@ protected:
 
     scalar totalHeatFlux_;
 
-    word QPartPartName_;
+    const word QPartPartName_;
 
     volScalarField QPartPart_;
 
@@ -64,23 +64,19 @@ protected:
 
     const volScalarField& partTempField_;
 
-    word prescribedVoidfractionFieldName_;
+    const word prescribedVoidfractionFieldName_;
 
     const volScalarField& prescribedVoidfraction_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    word partHeatFluxName_;
-
-    mutable double **partHeatFlux_;
+    const word partHeatFluxName_;
 
     scalarList typePartThermCond_;
 
-    mutable double **partThermCond_;
-
-    void allocateMyArrays() const;
+    const word partThermCondRegName_;
 
     void calcPartEffThermCond();
 
@@ -118,9 +114,9 @@ public:
 
         void addEnergyCoefficient(volScalarField&) const {}
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 
 };
 
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
index ece81a57..11064fa0 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
@@ -143,17 +143,25 @@ heatTransferGunn::heatTransferGunn
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","convectiveHeatFlux")),
-    partHeatFlux_(NULL),
-    partHeatFluxCoeff_(NULL),
-    partRe_(NULL),
-    partNu_(NULL),
+    partHeatFluxCoeffRegName_(typeName + "partHeatFluxCoeff"),
+    partReRegName_(typeName + "partRe"),
+    partNuRegName_(typeName + "partNu"),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partHeatFluxCoeffRegName_,1);
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partReRegName_,1);
+        particleCloud_.registerParticleProperty<double**>(partNuRegName_,1);
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -226,41 +234,16 @@ heatTransferGunn::heatTransferGunn
 
 heatTransferGunn::~heatTransferGunn()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partRe_,1);
-    particleCloud_.dataExchangeM().destroy(partNu_,1);
-    if (implicit_)
-    {
-        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferGunn::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    if(implicit_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
-    }
-
-    if(verbose_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferGunn::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     QPartFluid_.primitiveFieldRef() = 0.0;
@@ -387,6 +370,8 @@ void heatTransferGunn::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -433,6 +418,7 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -448,6 +434,8 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -515,6 +503,8 @@ scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) con
 void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid, scalar cg3)
 {
     scalar hAs = h * As * cg3;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     if (particleCloud_.getParticleEffVolFactors())
     {
@@ -529,6 +519,7 @@ void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid,
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -541,6 +532,7 @@ void heatTransferGunn::giveData()
         reduce(totalHeatFlux_, sumOp<scalar>());
         Info << "total convective particle-fluid heat flux [W] = " << totalHeatFlux_ << endl;
     }
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 
@@ -552,6 +544,9 @@ void heatTransferGunn::postFlow()
         scalar Tfluid(0.0);
         scalar Tpart(0.0);
         interpolationCellPoint<scalar> TInterpolator_(tempField_);
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+        double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
 
         totalHeatFlux_ = 0.0;
 
@@ -585,6 +580,7 @@ void heatTransferGunn::postFlow()
 
 void heatTransferGunn::partTempField()
 {
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
     partTempField_.primitiveFieldRef() = 0.0;
     particleCloud_.averagingM().resetWeightFields();
     particleCloud_.averagingM().setScalarAverage
@@ -607,6 +603,8 @@ void heatTransferGunn::initPartTemp()
 {
     label cellI = 0;
     scalar T = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
index 06b9a5ac..88f56189 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
@@ -68,11 +68,11 @@ protected:
 
     scalar NusseltScalingFactor_;
 
-    word QPartFluidName_;
+    const word QPartFluidName_;
 
     volScalarField QPartFluid_;
 
-    word QPartFluidCoeffName_;
+    const word QPartFluidCoeffName_;
 
     volScalarField QPartFluidCoeff_;
 
@@ -90,37 +90,33 @@ protected:
 
     dimensionedScalar partTempAve_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;           // ref to temperature field
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;        // ref to voidfraction field
 
     scalar maxSource_;                          // max (limited) value of src field
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;
+    const word partHeatFluxName_;
 
-    word partHeatFluxName_;
+    const word partHeatFluxCoeffRegName_;
 
-    mutable double **partHeatFlux_;
+    const word partReRegName_;
 
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
+    const word partNuRegName_;
 
     mutable scalar scaleDia_;
 
@@ -128,8 +124,6 @@ protected:
 
     const label maxTypeCG_;
 
-    void allocateMyArrays() const;
-
     void partTempField();
 
     scalar Nusselt(scalar, scalar, scalar) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
index fdfabc13..2b5a4e7c 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
@@ -142,17 +142,25 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","convectiveHeatFlux")),
-    partHeatFlux_(NULL),
-    partHeatFluxCoeff_(NULL),
-    partRe_(NULL),
-    partNu_(NULL),
+    partHeatFluxCoeffRegName_(typeName + "partHeatFluxCoeff"),
+    partReRegName_(typeName + "partRe"),
+    partNuRegName_(typeName + "partNu"),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partHeatFluxCoeffRegName_,1);
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partReRegName_,1);
+        particleCloud_.registerParticleProperty<double**>(partNuRegName_,1);
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -225,41 +233,16 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
 
 heatTransferRanzMarshall::~heatTransferRanzMarshall()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partRe_,1);
-    particleCloud_.dataExchangeM().destroy(partNu_,1);
-    if (implicit_)
-    {
-        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferRanzMarshall::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    if(implicit_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
-    }
-
-    if(verbose_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferRanzMarshall::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     QPartFluid_.primitiveFieldRef() = 0.0;
@@ -387,6 +370,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -427,6 +412,7 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -442,6 +428,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -503,8 +491,10 @@ scalar heatTransferRanzMarshall::Nusselt(scalar voidfraction, scalar Rep, scalar
 void heatTransferRanzMarshall::heatFlux(label index, scalar h, scalar As, scalar Tfluid, scalar cg3)
 {
     scalar hAs = h * As * cg3;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
-    if (particleCloud_.getParticleEffVolFactors()) 
+    if (particleCloud_.getParticleEffVolFactors())
     {
         scalar effVolFac = particleCloud_.particleEffVolFactor(index);
         hAs *= effVolFac;
@@ -517,6 +507,7 @@ void heatTransferRanzMarshall::heatFlux(label index, scalar h, scalar As, scalar
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -529,6 +520,7 @@ void heatTransferRanzMarshall::giveData()
         reduce(totalHeatFlux_, sumOp<scalar>());
         Info << "total convective particle-fluid heat flux [W] = " << totalHeatFlux_ << endl;
     }
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 
@@ -540,6 +532,9 @@ void heatTransferRanzMarshall::postFlow()
         scalar Tfluid(0.0);
         scalar Tpart(0.0);
         interpolationCellPoint<scalar> TInterpolator_(tempField_);
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+        double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
 
         totalHeatFlux_ = 0.0;
 
@@ -573,6 +568,7 @@ void heatTransferRanzMarshall::postFlow()
 
 void heatTransferRanzMarshall::partTempField()
 {
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
         partTempField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
@@ -595,6 +591,8 @@ void heatTransferRanzMarshall::initPartTemp()
 {
     label cellI = 0;
     scalar T = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H
index 17baedb8..8d9faa19 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H
@@ -67,11 +67,11 @@ protected:
 
     scalar NusseltScalingFactor_;
 
-    word QPartFluidName_;
+    const word QPartFluidName_;
 
     volScalarField QPartFluid_;
 
-    word QPartFluidCoeffName_;
+    const word QPartFluidCoeffName_;
 
     volScalarField QPartFluidCoeff_;
 
@@ -89,37 +89,33 @@ protected:
 
     dimensionedScalar partTempAve_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;            // ref to temperature field
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;    // ref to voidfraction field
 
     scalar maxSource_;                           // max (limited) value of src field
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;
+    const word partHeatFluxName_;
 
-    word partHeatFluxName_;
+    const word partHeatFluxCoeffRegName_;
 
-    mutable double **partHeatFlux_;
+    const word partReRegName_;
 
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
+    const word partNuRegName_;
 
     mutable scalar scaleDia_;
 
@@ -127,8 +123,6 @@ protected:
 
     const label maxTypeCG_;
 
-    void allocateMyArrays() const;
-
     void partTempField();
 
     scalar Nusselt(scalar, scalar, scalar) const;
@@ -163,11 +157,11 @@ public:
 
         void addEnergyContribution(volScalarField&) const;
 
-	void addEnergyCoefficient(volScalarField&) const;
+        void addEnergyCoefficient(volScalarField&) const;
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
index a5f3e2e1..71c30154 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
@@ -50,7 +50,6 @@ reactionHeat::reactionHeat
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
-    reactionHeat_(NULL),
     reactionHeatField_
     (
         IOobject
@@ -65,7 +64,7 @@ reactionHeat::reactionHeat
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0),0.0)
     )
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(reactionHeatName_,1);
 
     if(propsDict_.found("maxsource"))
     {
@@ -79,23 +78,15 @@ reactionHeat::reactionHeat
 
 reactionHeat::~reactionHeat()
 {
-    particleCloud_.dataExchangeM().destroy(reactionHeat_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void reactionHeat::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(reactionHeat_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void reactionHeat::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& reactionHeat_ = particleCloud_.getParticlePropertyRef<double**>(reactionHeatName_);
 
     particleCloud_.dataExchangeM().getData(reactionHeatName_,"scalar-atom",reactionHeat_);
 
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
index 4ebd7a02..87aff767 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
@@ -56,12 +56,8 @@ protected:
 
     word reactionHeatName_;
 
-    mutable double **reactionHeat_;
-
     volScalarField reactionHeatField_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C
index 496556d2..1c7a1963 100755
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C
@@ -74,8 +74,8 @@ KochHillRWDrag::KochHillRWDrag
     interpolation_(propsDict_.found("interpolation")),
     scale_(1.),
     randomTauE_(propsDict_.found("randomTauE")),
-    partTime_(NULL),
-    partUfluct_(NULL),
+    partTimeRegName_(typeName + "partTime"),
+    partUfluctRegName_(typeName + "partUfluct"),
     RanGen_(label(0))
 {
 
@@ -99,16 +99,8 @@ KochHillRWDrag::KochHillRWDrag
     if (propsDict_.found("rhoP"))
         rhoP_= readScalar(propsDict_.lookup("rhoP"));
 
-
-//    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-//    {
-        //allocate memory
-        particleCloud_.dataExchangeM().allocateArray(partTime_,0.,1);
-        particleCloud_.dataExchangeM().allocateArray(partUfluct_,0.,3);
-//    }
-
-    //Pout << "RW-TEST: maxNumberOfParticles() == " << particleCloud_.dataExchangeM().maxNumberOfParticles() << endl; // TEST-Output
-
+    particleCloud_.registerParticleProperty<double**>(partTimeRegName_,1,0.0,false);
+    particleCloud_.registerParticleProperty<double**>(partUfluctRegName_,3,0.0,false);
 }
 
 
@@ -116,18 +108,12 @@ KochHillRWDrag::KochHillRWDrag
 
 KochHillRWDrag::~KochHillRWDrag()
 {
-    particleCloud_.dataExchangeM().destroy(partTime_, 1);
-    particleCloud_.dataExchangeM().destroy(partUfluct_, 3);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void KochHillRWDrag::setForce() const
 {
-
-    // realloc the arrays
-    reAllocArrays();
-
     if (scale_ > 1.0)
     {
         Info << "KochHillRW using scale = " << scale_ << endl;
@@ -189,6 +175,9 @@ void KochHillRWDrag::setForce() const
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
 
+    double**& partTime_ = particleCloud_.getParticlePropertyRef<double**>(partTimeRegName_);
+    double**& partUfluct_ = particleCloud_.getParticlePropertyRef<double**>(partUfluctRegName_);
+
     //Info << "RW-TEST: We are in setForce() at t = " << t << endl; // TEST-Output
 
     for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
@@ -382,16 +371,6 @@ void KochHillRWDrag::setForce() const
 }
 
 
-void KochHillRWDrag::reAllocArrays() const
-{
-    if (particleCloud_.numberOfParticlesChanged())
-    {
-        particleCloud_.dataExchangeM().allocateArray(partTime_,0.0,1);  // field/initVal/with/lenghtFromLigghts
-        particleCloud_.dataExchangeM().allocateArray(partUfluct_,0.0,3);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H
index 5a8b762d..3aeffc37 100755
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H
@@ -89,9 +89,9 @@ private:
 
     const bool randomTauE_;
 
-    mutable double **partTime_;                // Lagrangian array
+    const word partTimeRegName_;
 
-    mutable double **partUfluct_;              // Lagrangian array
+    const word partUfluctRegName_;
 
     mutable Random RanGen_;
 
@@ -117,8 +117,6 @@ public:
 
     // Member Functions
         void setForce() const;
-
-        void reAllocArrays() const;
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
index bb8813fc..84604074 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C
@@ -71,14 +71,10 @@ LaEuScalarTemp::LaEuScalarTemp
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
-    partHeatFlux_(NULL),
     lambda_(readScalar(propsDict_.lookup("lambda"))),
     Cp_(readScalar(propsDict_.lookup("Cp")))
 {
-    allocateMyArrays();
-
     if (propsDict_.found("maxSource"))
     {
         maxSource_=readScalar(propsDict_.lookup ("maxSource"));
@@ -96,8 +92,10 @@ LaEuScalarTemp::LaEuScalarTemp
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
 
-
     particleCloud_.checkCG(false);
+
+    particleCloud_.registerParticleProperty<double**>(partTempName_);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_);
 }
 
 
@@ -105,8 +103,6 @@ LaEuScalarTemp::LaEuScalarTemp
 
 LaEuScalarTemp::~LaEuScalarTemp()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -114,6 +110,9 @@ void LaEuScalarTemp::allocateMyArrays() const
 {
     // get memory for 2d arrays
     double initVal = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
 }
@@ -128,6 +127,8 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
 {
     // realloc the arrays
     allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     EuField.primitiveFieldRef() = 0.0;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H
index c9de04f9..98fa42e2 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H
@@ -81,12 +81,8 @@ private:
 
     word partTempName_;
 
-    mutable double **partTemp_;                  // Lagrangian array
-
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;              // Lagrangian array
-
     scalar lambda_;                              // fluid thermal conductivity [W/(m*K)]
 
     scalar Cp_;                                  // specific heat capacity [W*s/(kg*K)]
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C
index 88fed5d4..f202af0f 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C
@@ -54,8 +54,8 @@ dSauter::dSauter
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     multiTypes_(false),
-    d2_(NULL),
-    d3_(NULL),
+    d2RegName_(typeName + "d2"),
+    d3RegName_(typeName + "d3"),
     maxTypeCG_(1),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     d2Field_
@@ -101,9 +101,11 @@ dSauter::dSauter
         multiTypes_ = true;
         maxTypeCG_ = typeCG_.size();
     }
-    allocateMyArrays();
-    dSauter_.write();
 
+    particleCloud_.registerParticleProperty<double**>(d2RegName_,1);
+    particleCloud_.registerParticleProperty<double**>(d3RegName_,1);
+
+    dSauter_.write();
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -114,19 +116,10 @@ dSauter::dSauter
 
 dSauter::~dSauter()
 {
-    particleCloud_.dataExchangeM().destroy(d2_,1);
-    particleCloud_.dataExchangeM().destroy(d3_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void dSauter::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(d2_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(d3_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
@@ -137,7 +130,8 @@ void dSauter::setForce() const
         Info << "dSauter using CG factor(s) = " << typeCG_ << endl;
     }
 
-    allocateMyArrays();
+    double**& d2_ = particleCloud_.getParticlePropertyRef<double**>(d2RegName_);
+    double**& d3_ = particleCloud_.getParticlePropertyRef<double**>(d3RegName_);
 
     label cellI = 0;
     label partType = 1;
@@ -151,11 +145,11 @@ void dSauter::setForce() const
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI >= 0)
         {
-            if (particleCloud_.getParticleEffVolFactors()) 
+            if (particleCloud_.getParticleEffVolFactors())
             {
                 effVolFac = particleCloud_.particleEffVolFactor(index);
             }
-            if (multiTypes_) 
+            if (multiTypes_)
             {
                 partType = particleCloud_.particleType(index);
                 if (partType > maxTypeCG_)
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H
index eb875ca5..843553d0 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H
@@ -45,21 +45,19 @@ private:
 
     bool multiTypes_;
 
-    mutable double **d2_;
-    
-    mutable double **d3_;
+    const word d2RegName_;
+
+    const word d3RegName_;
 
     label maxTypeCG_;
 
     scalarList typeCG_;
-    
-    mutable volScalarField d2Field_;
-    
-    mutable volScalarField d3Field_;
-    
-    mutable volScalarField dSauter_;
 
-    void allocateMyArrays() const;
+    mutable volScalarField d2Field_;
+
+    mutable volScalarField d3Field_;
+
+    mutable volScalarField dSauter_;
 
 public:
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C
index 072ea2de..88a5d2ee 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C
@@ -53,7 +53,7 @@ granKineticEnergy::granKineticEnergy
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    vfluc_(NULL),
+    vflucRegName_(typeName + "vfluc_mag"),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     granKineticEnergy_
@@ -70,7 +70,7 @@ granKineticEnergy::granKineticEnergy
         "zeroGradient"
     )
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(vflucRegName_,1);
     granKineticEnergy_.write();
 
 
@@ -82,21 +82,14 @@ granKineticEnergy::granKineticEnergy
 
 granKineticEnergy::~granKineticEnergy()
 {
-    particleCloud_.dataExchangeM().destroy(vfluc_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void granKineticEnergy::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,1);
-}
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
 void granKineticEnergy::setForce() const
 {
-    allocateMyArrays();
+    double**& vfluc_ = particleCloud_.getParticlePropertyRef<double**>(vflucRegName_);
 
     label cellI = 0;
     vector velfluc(0,0,0);
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
index 46a8a7ff..c16a4ffd 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
@@ -43,16 +43,14 @@ private:
 
     dictionary propsDict_;
 
-    mutable double **vfluc_;
+    const word vflucRegName_;
 
-    word UsFieldName_;
+    const word UsFieldName_;
 
     const volVectorField& UsField_;
 
     mutable volScalarField granKineticEnergy_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.C b/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.C
index ed1444c7..b2481866 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.C
@@ -60,9 +60,9 @@ particleDeformation::particleDeformation
     partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
     lowerBounds_(propsDict_.lookupOrDefault<scalarList>("lowerBounds",scalarList(1,-1.0))),
     upperBounds_(propsDict_.lookupOrDefault<scalarList>("upperBounds",scalarList(1,-1.0))),
-    partDeformations_(NULL)
+    partDeformationsName_("partDeformations")
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partDeformationsName_,1);
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -120,16 +120,9 @@ particleDeformation::particleDeformation
 
 particleDeformation::~particleDeformation()
 {
-    particleCloud_.dataExchangeM().destroy(partDeformations_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void particleDeformation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(partDeformations_,initVal,1);
-}
 
 bool particleDeformation::defaultDeformCell(label cell) const
 {
@@ -145,8 +138,8 @@ void particleDeformation::setForce() const
         init();
         initialExec_ = false;
     }
-    // realloc the arrays
-    allocateMyArrays();
+
+    double**& partDeformations_ = particleCloud_.getParticlePropertyRef<double**>(partDeformationsName_);
 
     label cellI = 0;
     label partType = -1;
@@ -204,7 +197,7 @@ void particleDeformation::setForce() const
     }
 
     // give DEM data
-    particleCloud_.dataExchangeM().giveData("partDeformations","scalar-atom", partDeformations_);
+    particleCloud_.dataExchangeM().giveData(partDeformationsName_,"scalar-atom", partDeformations_);
 }
 
 void particleDeformation::init() const
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.H b/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.H
index ce5718cd..84120c8f 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/particleDeformation/particleDeformation.H
@@ -55,12 +55,12 @@ private:
 
     mutable bool initialExec_;
 
-    word refFieldName_;
+    const word refFieldName_;
 
     mutable autoPtr<volScalarField> refField_;
 
     // default deformation in region
-    word defaultDeformCellsName_;
+    const word defaultDeformCellsName_;
 
     autoPtr<cellSet> defaultDeformCells_;
 
@@ -74,12 +74,10 @@ private:
 
     scalarList upperBounds_;
 
-    mutable double **partDeformations_;
+    const word partDeformationsName_;
 
     label getListIndex(label) const;
 
-    void allocateMyArrays() const;
-
     void init() const;
 
     bool defaultDeformCell(label) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C
index d79f03d2..65d64548 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C
@@ -18,7 +18,7 @@ SourceFiles
     pdCorrelation.C
 
 Contributing Author
-    2018    Paul Kieckhefen, TUHH 
+    2018    Paul Kieckhefen, TUHH
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
@@ -55,11 +55,11 @@ pdCorrelation::pdCorrelation
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    d_(nullptr),
-    p_(nullptr),
-    d2_(nullptr),
-    pd_(nullptr),
-    cg3_(nullptr),
+    dRegName_(typeName + "d"),
+    pRegName_(typeName + "p"),
+    d2RegName_(typeName + "d2"),
+    pdRegName_(typeName + "pd"),
+    cg3RegName_(typeName + "cg3"),
     dField_
     (   IOobject
         (
@@ -151,7 +151,11 @@ pdCorrelation::pdCorrelation
                   << abort(FatalError);
     }
 
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(dRegName_,1);
+    particleCloud_.registerParticleProperty<double**>(pRegName_,3);
+    particleCloud_.registerParticleProperty<double**>(d2RegName_,1);
+    particleCloud_.registerParticleProperty<double**>(pdRegName_,3);
+    particleCloud_.registerParticleProperty<double**>(cg3RegName_,1);
 
     dField_.write();
     pdField_.write();
@@ -165,24 +169,9 @@ pdCorrelation::pdCorrelation
 
 pdCorrelation::~pdCorrelation()
 {
-    particleCloud_.dataExchangeM().destroy(cg3_, 1);
-    particleCloud_.dataExchangeM().destroy(d_,  1);
-    particleCloud_.dataExchangeM().destroy(p_,  3);
-    particleCloud_.dataExchangeM().destroy(d2_, 1);
-    particleCloud_.dataExchangeM().destroy(pd_, 3);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void pdCorrelation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(d_,  initVal, 1);
-    particleCloud_.dataExchangeM().allocateArray(p_,  initVal, 3);
-    particleCloud_.dataExchangeM().allocateArray(d2_, initVal, 1);
-    particleCloud_.dataExchangeM().allocateArray(pd_, initVal, 3);
-    particleCloud_.dataExchangeM().allocateArray(cg3_, initVal, 1);
-}
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
 void pdCorrelation::setForce() const
@@ -191,7 +180,11 @@ void pdCorrelation::setForce() const
 
     if (runOnWriteOnly_ && !mesh.write()) return; // skip if it's not write time
 
-    allocateMyArrays();
+    double**& d_ = particleCloud_.getParticlePropertyRef<double**>(dRegName_);
+    double**& p_ = particleCloud_.getParticlePropertyRef<double**>(pRegName_);
+    double**& d2_ = particleCloud_.getParticlePropertyRef<double**>(d2RegName_);
+    double**& pd_ = particleCloud_.getParticlePropertyRef<double**>(pdRegName_);
+    double**& cg3_ = particleCloud_.getParticlePropertyRef<double**>(cg3RegName_);
 
     const Switch densityFromList
     (
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H
index 04fb2bf0..8c7d4d7d 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H
@@ -45,11 +45,11 @@ private:
 
     dictionary propsDict_;
 
-    mutable double **d_;
-    mutable double **p_;
-    mutable double **d2_;
-    mutable double **pd_;
-    mutable double **cg3_;
+    const word dRegName_;
+    const word pRegName_;
+    const word d2RegName_;
+    const word pdRegName_;
+    const word cg3RegName_;
 
     mutable volScalarField dField_;
     mutable volVectorField pField_;
@@ -65,8 +65,6 @@ private:
     const Switch     CG_;
     const Switch     runOnWriteOnly_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.C b/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.C
index ff8abc89..bfe83895 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.C
@@ -1,6 +1,6 @@
 /*---------------------------------------------------------------------------*\
     CFDEMcoupling academic - Open Source CFD-DEM coupling
-    
+
     Contributing authors:
     Thomas Lichtenegger
     Copyright (C) 2015- Johannes Kepler University, Linz
@@ -83,10 +83,10 @@ potentialRelaxation::potentialRelaxation
     ignoreReg_(propsDict_.lookupOrDefault<bool>("ignoreRegion",false)),
     ignoreDirection_(propsDict_.lookupOrDefault<vector>("ignoreDirection",vector::zero)),
     ignorePoint_(propsDict_.lookupOrDefault<vector>("ignorePoint",vector::zero)),
-    vfluc_(NULL)
+    vflucName_("vfluc")
 {
-    allocateMyArrays();
-    
+    particleCloud_.registerParticleProperty<double**>(vflucName_,3);
+
     if(ignoreReg_)
     {
         if(mag(ignoreDirection_) < SMALL)
@@ -103,48 +103,39 @@ potentialRelaxation::potentialRelaxation
 
 potentialRelaxation::~potentialRelaxation()
 {
-    delete vfluc_;
 }
 
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void potentialRelaxation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,3); 
-}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void potentialRelaxation::setForce() const
 {
-
     relax(D0_,D1_);
-    
+
     volVectorField relaxStream = -fvc::grad(correctedField_);
-    
+
   // volVectorField relaxStream = DField_ * fvc::grad(voidfraction_ - voidfractionRec_);
-    
-     // realloc the arrays
-    allocateMyArrays();
-    
+
+    double**& vfluc_ = particleCloud_.getParticlePropertyRef<double**>(vflucName_);
+
     vector position(0,0,0);
     scalar voidfraction(0.0);
     vector flucU(0,0,0);
     label cellI=0;
-   
+
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> relaxStreamInterpolator_(relaxStream);
-    
+
     scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
     scalar dtCFD = voidfraction_.mesh().time().deltaTValue();
-    
+
     // if DEM time step > CFD time step, scale velocity down
     scalar timeFac = 1.0;
     if (dtDEM > dtCFD) timeFac = dtCFD / dtDEM;
 
-    
+
     for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -164,18 +155,18 @@ void potentialRelaxation::setForce() const
                     {
                       position = particleCloud_.position(index);
                       voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
-                      flucU = relaxStreamInterpolator_.interpolate(position,cellI);      
+                      flucU = relaxStreamInterpolator_.interpolate(position,cellI);
                     }
                     else
                     {
                         voidfraction = voidfraction_[cellI];
                         flucU = relaxStream[cellI];
                     }
-                    
+
                     if (voidfraction > 1.0-SMALL) voidfraction = 1.0 - SMALL;
                     flucU /= (1-voidfraction);
-		    flucU *= timeFac;
-                    // write particle based data to global array 
+                    flucU *= timeFac;
+                    // write particle based data to global array
                     for(int i = 0; i < 3; i++)
                     {
                         vfluc_[index][i]=flucU[i];
@@ -183,13 +174,13 @@ void potentialRelaxation::setForce() const
                 }
             }
     }
-    
-    particleCloud_.dataExchangeM().giveData("vfluc","vector-atom", vfluc_);
-    
+
+    particleCloud_.dataExchangeM().giveData(vflucName_,"vector-atom", vfluc_);
+
     if (measureDiff_)
     {
         dimensionedScalar diff( fvc::domainIntegrate( sqr( voidfraction_ - voidfractionRec_ ) ) );
-        scalar t = particleCloud_.mesh().time().timeOutputValue(); 
+        scalar t = particleCloud_.mesh().time().timeOutputValue();
         recErrorFile_ << t << "\t" << diff.value() << endl;
     }
 }
@@ -198,12 +189,12 @@ void potentialRelaxation::relax(scalar D0, scalar D1) const
 {
     volScalarField src0 = voidfraction_ - voidfractionRec_;
     volScalarField src1 = voidfraction_ - voidfractionRec_;
-    
+
     forAll(src1, cellI)
     {
         if(src1[cellI] > 0.0) src1[cellI] = 0.0;
     }
-    
+
     solve
     (
         fvm::laplacian(correctedField_)
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.H b/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.H
index 81cebc5a..2bc51dfb 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/potentialRelaxation/potentialRelaxation.H
@@ -1,6 +1,6 @@
 /*---------------------------------------------------------------------------*\
     CFDEMcoupling academic - Open Source CFD-DEM coupling
-    
+
     Contributing authors:
     Thomas Lichtenegger
     Copyright (C) 2015- Johannes Kepler University, Linz
@@ -44,44 +44,42 @@ class potentialRelaxation
 {
 private:
     dictionary propsDict_;
-    
+
     bool interpolate_;
-    
+
     bool measureDiff_;
-    
+
     mutable OFstream recErrorFile_;
-    
-    word voidfractionFieldName_;
-    
+
+    const word voidfractionFieldName_;
+
     const volScalarField& voidfraction_;
-    
-    word voidfractionRecFieldName_;
-    
+
+    const word voidfractionRecFieldName_;
+
     const volScalarField& voidfractionRec_;
-    
+
     scalar critVoidfraction_;
-    
+
     scalar D0_;
-    
+
     scalar D1_;
-    
+
     mutable volScalarField correctedField_;
-    
+
     const scalar dt_;
-    
+
     // ignore particles in cells below plane given with ref point and normal vector
     // normal vector points towards region where fluctuations are permitted
-    
-    bool ignoreReg_;
-    
-    vector ignoreDirection_;
-    
-    vector ignorePoint_;
-    
-    mutable double **vfluc_;              // Lagrangian array
 
-    void allocateMyArrays() const;
-    
+    bool ignoreReg_;
+
+    vector ignoreDirection_;
+
+    vector ignorePoint_;
+
+    const word vflucName_;
+
     void relax(scalar, scalar) const;
 
 public:
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C
index 122b50e7..39dff7b2 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C
@@ -69,16 +69,11 @@ virtualMassForce::virtualMassForce
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     phiFieldName_(propsDict_.lookup("phiFieldName")),
     phi_(sm.mesh().lookupObject<surfaceScalarField> (phiFieldName_)),
-    UrelOld_(NULL),
+    UrelOldRegName_(typeName + "UrelOld"),
     splitUrelCalculation_(propsDict_.lookupOrDefault<bool>("splitUrelCalculation",false)),
     Cadd_(0.5)
 {
-
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d array
-        particleCloud_.dataExchangeM().allocateArray(UrelOld_,NOTONCPU,3);
-    }
+    particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false);
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -117,15 +112,13 @@ virtualMassForce::virtualMassForce
 
 virtualMassForce::~virtualMassForce()
 {
-    particleCloud_.dataExchangeM().destroy(UrelOld_,3);
 }
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void virtualMassForce::setForce() const
 {
-    reAllocArrays();
+    double**& UrelOld_ = particleCloud_.getParticlePropertyRef<double**>(UrelOldRegName_);
 
     scalar dt = U_.mesh().time().deltaT().value();
 
@@ -233,17 +226,6 @@ void virtualMassForce::setForce() const
     }
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::virtualMassForce::reAllocArrays() const
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        Pout << "virtualMassForce::reAllocArrays..." << endl;
-        particleCloud_.dataExchangeM().allocateArray(UrelOld_,NOTONCPU,3);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H
index 28597d0b..9a25ed2e 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H
@@ -71,7 +71,7 @@ private:
 
     const surfaceScalarField& phi_;
 
-    mutable double **UrelOld_;
+    const word UrelOldRegName_;
 
     const bool splitUrelCalculation_;      //indicator to split calculation of Urel between CFDEM and LIGGGHTS
                                            //requires the integration fix to take dv/dt into account!
@@ -100,8 +100,6 @@ public:
 
     // Member Functions
         void setForce() const;
-
-        void reAllocArrays() const;
 };
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/regionModel/allRegion/allRegion.C b/src/lagrangian/cfdemParticle/subModels/regionModel/allRegion/allRegion.C
index 6fe64953..8ed7f2c2 100644
--- a/src/lagrangian/cfdemParticle/subModels/regionModel/allRegion/allRegion.C
+++ b/src/lagrangian/cfdemParticle/subModels/regionModel/allRegion/allRegion.C
@@ -78,7 +78,6 @@ allRegion::~allRegion()
 
 void allRegion::defineRegion() const
 {
-    reAllocArrays();
     // do nothing
 }
 
diff --git a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C
index b1770d07..41397f8d 100644
--- a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C
+++ b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C
@@ -46,15 +46,6 @@ defineRunTimeSelectionTable(regionModel, dictionary);
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void regionModel::reAllocArrays() const
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        // get arrays of new length
-        particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1);
-        particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -66,16 +57,10 @@ regionModel::regionModel
 )
 :
     dict_(dict),
-    particleCloud_(sm),
-    inRegion_(NULL),
-    outRegion_(NULL)
+    particleCloud_(sm)
 {
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1);
-        particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1);
-    }
+    particleCloud_.registerParticleProperty<double**>("inRegion",1,1.0,false);
+    particleCloud_.registerParticleProperty<double**>("outRegion",1,1.0,false);
 }
 
 
@@ -83,8 +68,6 @@ regionModel::regionModel
 
 regionModel::~regionModel()
 {
-    particleCloud_.dataExchangeM().destroy(inRegion_,1);
-    particleCloud_.dataExchangeM().destroy(outRegion_,1);
 }
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.H b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.H
index c087cf1f..735329af 100644
--- a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.H
+++ b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.H
@@ -60,10 +60,6 @@ protected:
 
         cfdemCloud& particleCloud_;
 
-        mutable double **inRegion_;
-
-        mutable double **outRegion_;
-
 public:
 
     friend class voidFractionModel;
@@ -117,13 +113,6 @@ public:
 
         virtual void resetVolFields(volVectorField&) const = 0;
 
-        void reAllocArrays() const;
-
-    // Access
-        inline double ** const& inRegion()const{ return inRegion_; };
-
-        inline double ** const& outRegion()const { return outRegion_; };
-
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
index 7a5b1ce9..3c6e2949 100644
--- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
+++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
@@ -55,14 +55,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     typeKs_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
-    partKs_(NULL)
+    partKsRegName_(typeName + "partKs")
 {
     if (typeKs_[0] < 0.0)
     {
         FatalError << "ZehnerSchluenderThermCond: provide list of thermal conductivities." << abort(FatalError);
     }
 
-    if (typeKs_.size() > 1) allocateMyArrays();
+    if (typeKs_.size() > 1)
+    {
+        particleCloud_.registerParticleProperty<double**>(partKsRegName_,1);
+    }
     else
     {
         partKsField_ *= typeKs_[0];
@@ -74,14 +77,6 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
 
 ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond()
 {
-    if (typeKs_.size() > 1) particleCloud_.dataExchangeM().destroy(partKs_,1);
-}
-
-
-void ZehnerSchluenderThermCond::allocateMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partKs_,initVal,1);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -133,7 +128,7 @@ void ZehnerSchluenderThermCond::calcPartKsField() const
         FatalError << "ZehnerSchluenderThermCond needs data for more than one type, but types are not communicated." << abort(FatalError);
     }
 
-    allocateMyArrays();
+    double**& partKs_ = particleCloud_.getParticlePropertyRef<double**>(partKsRegName_);
     label cellI=0;
     label partType = 0;
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H
index b681483e..6d10ccc9 100644
--- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H
+++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H
@@ -54,19 +54,17 @@ private:
 
     dictionary propsDict_;
 
-    word partKsFieldName_;
+    const word partKsFieldName_;
 
     volScalarField& partKsField_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
     scalarList typeKs_;
 
-    mutable double **partKs_;
-
-    void allocateMyArrays() const;
+    const word partKsRegName_;
 
     void calcPartKsField() const;
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
index 9a8c5a22..eedd4f8a 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C
@@ -89,8 +89,6 @@ GaussVoidFraction::~GaussVoidFraction()
 
 void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     voidfractionNext_.ref()=1;
 
     scalar radius(-1);
@@ -103,12 +101,12 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
         //if(mask[index][0])
         //{
             //reset
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for(int subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
             {
                 particleWeights[index][subcell]=0;
                 particleVolumes[index][subcell]=0;
             }
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -142,7 +140,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
 
                     //making sure that the cell containing the center is the first subcell
                     particleCloud_.cellIDs()[index][0]=particleCenterCellID;
@@ -200,7 +198,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
     //bringing eulerian field to particle array
     for(label index=0; index< particleCloud_.numberOfParticles(); index++)
     {
-        for(label subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+        for(label subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
index fcc1947b..7f0172d0 100755
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C
@@ -91,8 +91,6 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 {
     const boundBox& globalBb = particleCloud_.mesh().bounds();
 
-    reAllocArrays();
-
     voidfractionNext_.ref() = 1.0;
 
     for (int index=0; index < particleCloud_.numberOfParticles(); index++)
@@ -100,13 +98,13 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
         //if(mask[index][0])
         //{
             //reset
-            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+            for (int subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
             {
                 particleWeights[index][subcell] = 0.0;
                 particleVolumes[index][subcell] = 0.0;
             }
 
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -266,7 +264,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
                     hashSett.erase(particleCenterCellID);
 
                     for (label i=0; i < hashSetLength-1; i++)
@@ -281,7 +279,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     for (label index=0; index < particleCloud_.numberOfParticles(); index++)
     {
-        for (label subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+        for (label subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
index 4d254b42..44ab10e9 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C
@@ -88,8 +88,6 @@ bigParticleVoidFraction::~bigParticleVoidFraction()
 
 void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     voidfractionNext_.ref()=1;
 
     scalar radius(-1);
@@ -102,12 +100,12 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
         //if(mask[index][0])
         //{
             //reset
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for(int subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
             {
                 particleWeights[index][subcell]=0;
                 particleVolumes[index][subcell]=0;
             }
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -137,7 +135,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
 
                     //making sure that the cell containing the center is the first subcell
                     particleCloud_.cellIDs()[index][0]=particleCenterCellID;
@@ -190,7 +188,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
     //bringing eulerian field to particle array
     for(label index=0; index< particleCloud_.numberOfParticles(); index++)
     {
-        for(label subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+        for(label subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C
index 355211cc..1c57f970 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C
@@ -80,8 +80,6 @@ centreVoidFraction::~centreVoidFraction()
 
 void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     scalar radius(-1);
     scalar volume(0);
     scalar cellVol(0);
@@ -93,7 +91,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
         //{
             // reset
             particleWeights[index][0]=0;
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
 
             label cellI = particleCloud_.cellIDs()[index][0];
 
@@ -128,9 +126,9 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
                 if (index==0)
                 {
                     Info << "centre cellI = " << cellI << endl;
-                    Info << "cellsPerParticle_=" << cellsPerParticle_[index][0] << endl;
+                    Info << "cellsPerParticle =" << cellsPerParticle()[index][0] << endl;
 
-                    for(int i=0;i<cellsPerParticle_[index][0];i++)
+                    for(int i=0;i<cellsPerParticle()[index][0];i++)
                     {
                        if(i==0)Info << "cellids, voidfractions, particleWeights, : \n";
                        Info << particleCloud_.cellIDs()[index][i] << " ," << endl;
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C
index 979bf064..d345f03d 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C
@@ -69,8 +69,7 @@ dividedVoidFraction::dividedVoidFraction
     alphaMin_(readScalar(propsDict_.lookup("alphaMin"))),
     alphaLimited_(0),
     tooMuch_(0.0),
-    interpolation_(propsDict_.found("interpolation")),
-    cfdemUseOnly_(propsDict_.lookupOrDefault<bool>("cfdemUseOnly", false))
+    interpolation_(propsDict_.found("interpolation"))
 {
     maxCellsPerParticle_ = numberOfMarkerPoints;
 
@@ -158,10 +157,6 @@ dividedVoidFraction::~dividedVoidFraction()
 
 void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes, double**& particleV)
 {
-    if (cfdemUseOnly_)
-        reAllocArrays(particleCloud_.numberOfParticles());
-    else
-        reAllocArrays();
 
     vector position(0.,0.,0.);
     label cellID = -1;
@@ -180,14 +175,14 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
         //{
             // reset
 
-            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+            for (int subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
             {
                 particleWeights[index][subcell] = 0.;
                 particleVolumes[index][subcell] = 0.;
             }
             particleV[index][0] = 0.;
 
-            cellsPerParticle_[index][0] = 1;
+            cellsPerParticle()[index][0] = 1;
             position = particleCloud_.position(index);
             cellID = particleCloud_.cellIDs()[index][0];
             radius = particleCloud_.radius(index);
@@ -270,7 +265,7 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
     for(int index=0; index < particleCloud_.numberOfParticles(); index++)
     {
         {
-            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+            for (int subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
             {
                 label cellID = particleCloud_.cellIDs()[index][subcell];
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H
index de8c5409..11fac124 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H
@@ -74,8 +74,6 @@ private:
 
         const bool interpolation_;
 
-        const bool cfdemUseOnly_;
-
         vector offsets[numberOfMarkerPoints];
 
         virtual inline scalar Vp(int index, scalar radius, scalar scaleVol) const
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H
index d8a418db..b0e0dabc 100755
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H
@@ -16,7 +16,7 @@
     }
 
     label partCellId = particleCloud_.locateM().findSingleCell(subPosition,cellID);
-               
+
     //NP fprintf(lmp->screen,"cellID=%d, partCellId=%d\n",static_cast<int>(cellID),static_cast<int>(partCellId));
 
     if(partCellId!=cellID)
@@ -35,7 +35,7 @@
         scalar partCellVol = particleCloud_.mesh().V()[partCellId];
         scalar particleVolume = volume/static_cast<scalar>(nPoints);
         scalar newAlpha = voidfractionNext_[partCellId]- particleVolume / partCellVol;
-        
+
         if(newAlpha > alphaMin_) voidfractionNext_[partCellId] = newAlpha;
         else
         {
@@ -49,9 +49,9 @@
         // add sub particle representation
         bool createNew = true;
         label storeInIndex=0;
-        for(int i=0; i < cellsPerParticle_[index][0] ; i++)
+        for(int i=0; i < cellsPerParticle()[index][0] ; i++)
         {
-            if(partCellId == particleCloud_.cellIDs()[index][i]) 
+            if(partCellId == particleCloud_.cellIDs()[index][i])
             {
                 storeInIndex = i;
                 createNew = false;
@@ -61,8 +61,8 @@
 
         if(createNew)
         {
-            cellsPerParticle_[index][0] ++;
-            storeInIndex = cellsPerParticle_[index][0]-1;
+            cellsPerParticle()[index][0] ++;
+            storeInIndex = cellsPerParticle()[index][0]-1;
             particleCloud_.cellIDs()[index][storeInIndex] = partCellId;
 
         }
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C
index 7cbc8bbd..18d12a79 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C
@@ -90,8 +90,6 @@ trilinearVoidFraction::~trilinearVoidFraction()
 
 void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     scalar radius(-1.);
     scalar volume(0.);
     scalar scaleVol = weight();
@@ -131,7 +129,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         // reset
-        cellsPerParticle_[index][0] = 8;
+        cellsPerParticle()[index][0] = 8;
         //TODO do we need to set particleVolumes, particleV?
         // ===
 
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
index b07d9667..4ac07188 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C
@@ -84,12 +84,12 @@ voidFractionModel::voidFractionModel
         /*sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 1)*/
     ),
-    cellsPerParticle_(NULL),
+    partCellsRegName_("cellsPerParticle"),
     maxCellsPerParticle_(1),
     weight_(1.),
     porosity_(1.)
 {
-    particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1);
+    particleCloud_.registerParticleProperty<int**>(partCellsRegName_,1,1.0,false);
     if (particleCloud_.getParticleEffVolFactors()) multiWeights_ = true;
 }
 
@@ -128,19 +128,18 @@ voidFractionModel::voidFractionModel
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0), initVoidfraction)
     ),
-    cellsPerParticle_(NULL),
+    partCellsRegName_("cellsPerParticle"),
     maxCellsPerParticle_(1),
     weight_(1.),
     porosity_(1.)
 {
-    particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1);
+    particleCloud_.registerParticleProperty<int**>(partCellsRegName_,1,1.0,false);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 voidFractionModel::~voidFractionModel()
 {
-    particleCloud_.dataExchangeM().destroy(cellsPerParticle_,1);
 }
 
 // * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
@@ -162,7 +161,7 @@ void voidFractionModel::resetVoidFractions()
 
 int** const& voidFractionModel::cellsPerParticle() const
 {
-    return cellsPerParticle_;
+    return particleCloud_.getParticlePropertyRef<int**>(partCellsRegName_);
 }
 
 int voidFractionModel::maxCellsPerParticle() const
@@ -170,24 +169,6 @@ int voidFractionModel::maxCellsPerParticle() const
     return maxCellsPerParticle_;
 }
 
-void voidFractionModel::reAllocArrays()
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        // get arrays of new length
-        particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1);
-    }
-}
-
-void voidFractionModel::reAllocArrays(int nP)
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        // get arrays of new length
-        particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1,nP);
-    }
-}
-
 scalar voidFractionModel::pointInParticle(int index, const vector& positionCenter, const vector& point, double scale) const
 {
     const scalar radius = particleCloud_.radius(index);
diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H
index 5e148fb2..ebc4dabc 100644
--- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H
+++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H
@@ -69,7 +69,7 @@ protected:
 
         volScalarField voidfractionNext_;
 
-        int ** cellsPerParticle_;
+        const word partCellsRegName_;
 
         int maxCellsPerParticle_;
 
@@ -165,10 +165,6 @@ public:
 
         int maxCellsPerParticle() const;
 
-        void reAllocArrays();
-
-        void reAllocArrays(int nP);  //force number of particles during reallocation, for CFD offline-use
-
         virtual void setParticleType(label type) const {}
 
         virtual bool checkParticleType(label) const { return true; } //consider all particles by default