diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
index 9ce61c26..f5a8a329 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**>("partPressure");
+    for (int i=0; i<diffusantGasNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>(diffusantGasNames_[i]);
+    }
     createCoeffs();
     molWeightTable();
 }
@@ -88,35 +90,17 @@ 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");
+    particleCloud_.dataExchangeM().allocateArray(particleCloud_.getParticlePropertyRef<double**>("partPressure"),initVal,1);
     for (int i=0; i<diffusantGasNames_.size(); i++)
     {
-        particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[i]),initVal,1);
     }
 }
 
@@ -172,6 +156,8 @@ void diffusionCoefficient::execute()
     interpolationCellPoint <scalar> TInterpolator_(tempField_);
     interpolationCellPoint <scalar> PInterpolator_(P_);
 
+    double**& partPressure_ = particleCloud_.getParticlePropertyRef<double**>("partPressure");
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
@@ -207,6 +193,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 +230,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 +245,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 +255,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..add7be45 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();
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
index 2f49a729..053aeb4c 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
@@ -63,42 +63,27 @@ 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),
     scaleDia_(1)
 {
     particleCloud_.checkCG(true);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partNuName_);
+    particleCloud_.registerParticleProperty<double**>(partReynolds_);
 }
 
 // * * * * * * * * * * * * * * * * 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);
+    particleCloud_.dataExchangeM().allocateArray(particleCloud_.getParticlePropertyRef<double**>(partNuName_),initVal,1);
+    particleCloud_.dataExchangeM().allocateArray(particleCloud_.getParticlePropertyRef<double**>(partReynolds_),initVal,1);
 }
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
@@ -139,6 +124,9 @@ void massTransferCoeff::execute()
     interpolationCellPoint <vector> UluidInterpolator_(U_);
     interpolationCellPoint <scalar> voidfractionInterpolator_(voidfraction_);
 
+    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuName_);
+    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReynolds_);
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI=particleCloud_.cellIDs()[index][0];
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
index 5021fbc8..8fc14ce1 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.H
@@ -70,14 +70,8 @@ private:
 
     word partNuName_;
 
-    mutable double **partNu_;
-
     word partReynolds_;
 
-    mutable double **partRe_;
-
-    void allocateMyArrays() const;
-
     mutable scalar scaleDia_;
 
 public:
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
index 0d2151c7..eabe40a7 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.C
@@ -57,7 +57,6 @@ reactantPerParticle::reactantPerParticle
     propsDict_(dict.subDict(typeName + "Props")),
     mesh_(sm.mesh()),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
-    reactantPerParticle_(NULL),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     particlesPerCell_
@@ -76,30 +75,21 @@ reactantPerParticle::reactantPerParticle
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1))
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>("reactantPerParticle");
 }
 
 // * * * * * * * * * * * * * * * * 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;
+    double**& reactantPerParticle_ = particleCloud_.getParticlePropertyRef<double**>("reactantPerParticle");
     particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1);
 }
 
@@ -121,6 +111,7 @@ void reactantPerParticle::execute()
     scalar voidfraction(1);
     scalar cellvolume(0.0);
     scalar particlesPerCell(1.0);
+    double**& reactantPerParticle_ = particleCloud_.getParticlePropertyRef<double**>("reactantPerParticle");
 
     // first create particles per cell field
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
@@ -147,10 +138,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("reactantPerParticle", "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..163dbb19 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/reactantPerParticle/reactantPerParticle.H
@@ -56,16 +56,12 @@ private:
 
     bool verbose_;
 
-    mutable double **reactantPerParticle_;
-
     word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
     mutable volScalarField particlesPerCell_;
 
-    void allocateMyArrays() const;
-
     label loopCounter_;
 
     label Nevery_;
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
index c485a848..d37773af 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,53 +103,42 @@ species::species
     initialized_(false)
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_);
+    particleCloud_.registerParticleProperty<double**>(partRhoName_);
+    particleCloud_.registerParticleProperty<double**>(partMolarConcName_);
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>("X_"+speciesNames_[i]);
+        particleCloud_.registerParticleProperty<double**>("Modified_"+speciesNames_[i]);
+    }
 }
 
 // * * * * * * * * * * * * * * * * 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;
+    double**& partRho_ = particleCloud_.getParticlePropertyRef<double**>(partRhoName_);
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partMolarConc_ = particleCloud_.getParticlePropertyRef<double**>(partMolarConcName_);
+
     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);
+        double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
+        double**& changeOfSpeciesMass_ = particleCloud_.getParticlePropertyRef<double**>("Modified_"+speciesNames_[i]);
+        particleCloud_.dataExchangeM().allocateArray(molarFractions_,initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_,initVal,1);
     }
 }
 
@@ -235,6 +219,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 +251,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 +262,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 +281,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 +295,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..a1119945 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
@@ -68,10 +68,6 @@ private:
 
     UPtrList<volScalarField> X_;
 
-    mutable List<double**> molarFractions_;
-
-    mutable List<double**> changeOfSpeciesMass_;
-
     PtrList<volScalarField> changeOfSpeciesMassFields_;
 
     volScalarField changeOfGasMassField_;
@@ -82,22 +78,16 @@ private:
 
     word partTempName_;
 
-    mutable double **partTemp_;                    // gas temperature at particle positions
-
     word densityFieldName_;
 
     const volScalarField& rho_;
 
     word partRhoName_;
 
-    mutable double **partRho_;                     // gas density at particle positions
-
     word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    void allocateMyArrays() const;
-
     // total mole field
     word molarConcFieldName_;
 
@@ -105,8 +95,6 @@ private:
 
     word partMolarConcName_;
 
-    mutable double **partMolarConc_;
-
     label loopCounter_;
 
     label Nevery_;
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C
index 135850f7..5c8ed8cd 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,10 @@ 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)
+    typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0)))
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_);
+    particleCloud_.registerParticleProperty<double**>("partThermCond");
 
     if (typePartThermCond_[0] < 0.0)
     {
@@ -127,8 +127,6 @@ heatTransferGranConduction::heatTransferGranConduction
 
 heatTransferGranConduction::~heatTransferGranConduction()
 {
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partThermCond_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -136,6 +134,9 @@ void heatTransferGranConduction::allocateMyArrays() const
 {
     // get memory for 2d arrays
     double initVal=0.0;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>("partThermCond");
+
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
     particleCloud_.dataExchangeM().allocateArray(partThermCond_,initVal,1);
 }
@@ -156,6 +157,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 +201,8 @@ void heatTransferGranConduction::calcPartThermCond()
 {
     label cellI=0;
     label partType = 1;
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>("partThermCond");
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -227,7 +231,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 +243,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..f5d28d8e 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
@@ -74,12 +74,8 @@ protected:
 
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;
-
     scalarList typePartThermCond_;
 
-    mutable double **partThermCond_;
-
     void allocateMyArrays() const;
 
     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..08cea07f 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
@@ -143,17 +143,22 @@ 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),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>("partHeatFluxCoeff");
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>("partRe");
+        particleCloud_.registerParticleProperty<double**>("partNu");
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -226,14 +231,6 @@ 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  * * * * * * * * * * * * * //
@@ -241,15 +238,22 @@ void heatTransferGunn::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);
+
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
     }
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
         particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
         particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
     }
@@ -261,6 +265,8 @@ 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 +393,8 @@ void heatTransferGunn::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -433,6 +441,7 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -448,6 +457,8 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -515,6 +526,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 +542,7 @@ void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid,
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -541,6 +555,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 +567,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**>("partHeatFluxCoeff");
 
         totalHeatFlux_ = 0.0;
 
@@ -585,6 +603,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 +626,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..037c2711 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
@@ -110,18 +110,8 @@ protected:
 
     word partTempName_;
 
-    mutable double **partTemp_;
-
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;
-
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
-
     mutable scalar scaleDia_;
 
     scalarList typeCG_;
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
index fdfabc13..a08c16e5 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C
@@ -142,17 +142,22 @@ 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),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>("partHeatFluxCoeff");
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>("partRe");
+        particleCloud_.registerParticleProperty<double**>("partNu");
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -225,14 +230,6 @@ 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  * * * * * * * * * * * * * //
@@ -240,15 +237,22 @@ void heatTransferRanzMarshall::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);
+
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
     }
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
         particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
         particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
     }
@@ -260,6 +264,8 @@ 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 +393,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -427,6 +435,7 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -442,6 +451,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>("partRe");
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>("partNu");
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -503,8 +514,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 +530,7 @@ void heatTransferRanzMarshall::heatFlux(label index, scalar h, scalar As, scalar
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>("partHeatFluxCoeff");
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -529,6 +543,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 +555,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**>("partHeatFluxCoeff");
 
         totalHeatFlux_ = 0.0;
 
@@ -573,6 +591,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 +614,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..10ab27f9 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H
@@ -109,18 +109,8 @@ protected:
 
     word partTempName_;
 
-    mutable double **partTemp_;
-
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;
-
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
-
     mutable scalar scaleDia_;
 
     scalarList typeCG_;
@@ -163,11 +153,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..52a10b5d 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_);
 
     if(propsDict_.found("maxsource"))
     {
@@ -79,7 +78,6 @@ reactionHeat::reactionHeat
 
 reactionHeat::~reactionHeat()
 {
-    particleCloud_.dataExchangeM().destroy(reactionHeat_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -87,6 +85,7 @@ void reactionHeat::allocateMyArrays() const
 {
     // get memory for 2d arrays
     double initVal=0.0;
+    double**& reactionHeat_ = particleCloud_.getParticlePropertyRef<double**>(reactionHeatName_);
     particleCloud_.dataExchangeM().allocateArray(reactionHeat_,initVal,1);
 }
 
@@ -96,6 +95,7 @@ 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..bdd3d196 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
@@ -56,8 +56,6 @@ protected:
 
     word reactionHeatName_;
 
-    mutable double **reactionHeat_;
-
     volScalarField reactionHeatField_;
 
     void allocateMyArrays() const;