diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
index b04cd4f6..b9cf2fcb 100644
--- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
+++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
@@ -51,7 +51,7 @@ Description
 int main(int argc, char *argv[])
 {
     #include "postProcess.H"
-  
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -60,13 +60,14 @@ int main(int argc, char *argv[])
     #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createFieldRefs.H"
     #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
     cfdemCloudEnergy particleCloud(mesh);
     #include "checkModelType.H"
-    
+
     turbulence->validate();
   //        #include "compressibleCourantNo.H"
   //  #include "setInitialDeltaT.H"
@@ -82,12 +83,12 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         runTime++;
-	
-	particleCloud.clockM().start(1,"Global");
+
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-	// do particle stuff
+        // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
@@ -95,7 +96,7 @@ int main(int argc, char *argv[])
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
         }
-    
+
         Info << "update Ksl.internalField()" << endl;
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
@@ -111,12 +112,12 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-	
+
         if (pimple.nCorrPIMPLE() <= 1)
         {
             #include "rhoEqn.H"
         }
-        
+
         volScalarField rhoeps("rhoeps",rho*voidfraction);
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
@@ -127,9 +128,9 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-	        // besides this pEqn, OF offers a "pimple consistent"-option
+                // besides this pEqn, OF offers a "pimple consistent"-option
                 #include "pEqn.H"
-	      	rhoeps=rho*voidfraction;	
+                rhoeps=rho*voidfraction;
             }
 
             if (pimple.turbCorr())
@@ -140,7 +141,7 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-       
+
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
diff --git a/applications/solvers/cfdemSolverRhoPimple/createFieldRefs.H b/applications/solvers/cfdemSolverRhoPimple/createFieldRefs.H
new file mode 100644
index 00000000..5842906a
--- /dev/null
+++ b/applications/solvers/cfdemSolverRhoPimple/createFieldRefs.H
@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/cfdemSolverRhoPimple/createFields.H b/applications/solvers/cfdemSolverRhoPimple/createFields.H
index eac9a5a7..366cc1d0 100644
--- a/applications/solvers/cfdemSolverRhoPimple/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/createFields.H
@@ -7,8 +7,6 @@ Info<< "Reading thermophysical properties\n" << endl;
     psiThermo& thermo = pThermo();
     thermo.validate(args.executable(), "h", "e");
     volScalarField& p = thermo.p();
-    volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field rho\n" << endl;
     volScalarField rho
@@ -25,7 +23,6 @@ Info<< "Reading thermophysical properties\n" << endl;
     );
 
 
-    
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -39,7 +36,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         ),
         mesh
     );
-    
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
@@ -81,7 +78,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh,
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
-    
+
     Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
     volScalarField QCoeff
     (
@@ -96,7 +93,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh,
         dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
     );
-    
+
   /*  Info<< "\nCreating thermal diffusivity field\n" << endl;
     volScalarField thDiff
     (
@@ -141,7 +138,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh,
         dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
     );
-    
+
     Info<< "\nCreating body force field\n" << endl;
     volVectorField fOther
     (
@@ -204,7 +201,7 @@ Info<< "Reading thermophysical properties\n" << endl;
             thermo
         )
     );
-    
+
     mesh.setFluxRequired(p.name());
 
     Info<< "Creating field dpdt\n" << endl;
@@ -222,7 +219,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
-    
+
     Info<< "\nReading momentum exchange field Ksl\n" << endl;
     volScalarField Ksl
     (
diff --git a/src/lagrangian/cfdemParticle/etc/solver-list.txt b/src/lagrangian/cfdemParticle/etc/solver-list.txt
index 3ebcf8c5..9d9911fa 100644
--- a/src/lagrangian/cfdemParticle/etc/solver-list.txt
+++ b/src/lagrangian/cfdemParticle/etc/solver-list.txt
@@ -1,6 +1,5 @@
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
 cfdemSolverRhoPimple/dir
-cfdemSolverRhoPimpleChem/dir
 cfdemSolverIB/dir
-cfdemSolverPisoScalar/dir
\ No newline at end of file
+cfdemSolverPisoScalar/dir
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
index 8890bf8b..a3cf1635 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
@@ -97,8 +97,8 @@ void reactionHeat::calcEnergyContribution()
 
     particleCloud_.dataExchangeM().getData(reactionHeatName_,"scalar-atom",reactionHeat_);
 
-    reactionHeatField_.internalField() = 0.0;
-    reactionHeatField_.boundaryField() = 0.0;
+    reactionHeatField_.primitiveFieldRef() = 0.0;
+    reactionHeatField_.boundaryFieldRef() = 0.0;
 
     particleCloud_.averagingM().setScalarSum
     (
@@ -108,7 +108,7 @@ void reactionHeat::calcEnergyContribution()
         NULL
     );
 
-    reactionHeatField_.internalField() /= (reactionHeatField_.mesh().V());
+    reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V());
 
     forAll(reactionHeatField_,cellI)
     {
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C
index c36fafbb..17fffd4e 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C
@@ -107,7 +107,7 @@ ErgunStatFines::~ErgunStatFines()
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-const scalar ErgunStatFines::dSauter(label cellI) const
+scalar ErgunStatFines::dSauter(label cellI) const
 {
     // Sauter mean diameter without influence of medium-scale fines
     scalar dS = dSauter_[cellI] / scaleDist_;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H
index 108e0841..1b33e30e 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H
@@ -73,7 +73,7 @@ private:
 
     mutable scalar switchingVoidfraction_; //voidfraction above which dilute formulation will be used
 
-    inline const scalar dSauter(label) const;
+    inline scalar dSauter(label) const;
 
 public:
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
index 8cc189b6..46a8a7ff 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.H
@@ -44,12 +44,12 @@ private:
     dictionary propsDict_;
 
     mutable double **vfluc_;
-    
-    mutable volScalarField granKineticEnergy_;
-    
-    word UsFieldName_;   
 
-    const volVectorField& UsField_; 
+    word UsFieldName_;
+
+    const volVectorField& UsField_;
+
+    mutable volScalarField granKineticEnergy_;
 
     void allocateMyArrays() const;
 
diff --git a/src/lagrangian/cfdemParticleComp/Make/options b/src/lagrangian/cfdemParticleComp/Make/options
index e92ad1aa..1bbb8cc5 100644
--- a/src/lagrangian/cfdemParticleComp/Make/options
+++ b/src/lagrangian/cfdemParticleComp/Make/options
@@ -12,7 +12,7 @@ PFLAGS+= -Dcompre
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
-    -I ../cfdemParticle/lnInclude \
+    -I../cfdemParticle/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \