diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/Make/options b/applications/solvers/Lagrangian/reactingParcelFoam/Make/options
index b0e474c308f5b7208b8f099789ca3247bae306fa..7d9fba61432e98b3f2cf2b390e4dcc547d488f8a 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/Make/options
@@ -4,7 +4,6 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
@@ -19,7 +18,6 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN) \
     -lfiniteVolume \
     -lmeshTools \
     -lcompressibleRASModels \
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/Lagrangian/reactingParcelFoam/UEqn.H
index ce4292dec9c933252515495d35babaec14835a7e..3c4a927091a474e5b8d0d4e50b2794a26202cc8a 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/UEqn.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/UEqn.H
@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + reactingParcels.SU1()
+      + parcels.SU()
     );
 
     UEqn.relax();
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/Lagrangian/reactingParcelFoam/YEqn.H
index d7dd43e6d70a22d354677234e06eaaee1c60b578..c7ade57fbeaf33effb8d4fef4abd4bda4b0e0551 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/YEqn.H
@@ -1,4 +1,3 @@
-
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
@@ -26,7 +25,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
-                reactingParcels.Srho1(i)
+                parcels.Srho(i)
               + kappa*chemistry.RR(i)().dimensionedInternalField()
             );
 
@@ -41,5 +40,4 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
-
 }
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/Lagrangian/reactingParcelFoam/createClouds.H
index 238ee3b18b7e84f73ccd9c21c6b11038669b4880..b4ec5e40569b90994c4575e1faae10cd304d1c78 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/createClouds.H
@@ -1,5 +1,3 @@
-Info<< "\nConstructing interpolation" << endl;
-
 Info << "\nConstructing gas properties" << endl;
 /*
 PtrList<specieConstProperties> gasProperties(Y.size());
@@ -31,7 +29,7 @@ forAll(gasProperties, i)
 }
 
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud reactingParcels
+basicReactingCloud parcels
 (
     "reactingCloud1",
     rho,
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/Lagrangian/reactingParcelFoam/createFields.H
index 079a4f2e1d1ba7b48486909490f1abcd7a7f0352..8f0d59e77d4d2e88e6080b7b029dcb1facefa450 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/createFields.H
@@ -28,51 +28,6 @@
         thermo->rho()
     );
 
-// lagrangian coal density field
-/*    volScalarField rholc
-    (
-        IOobject
-        (
-            "rholc",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );
-
-// lagrangian limestone density field
-    volScalarField rhols
-    (
-        IOobject
-        (
-            "rhols",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );
-
-// lagrangian total density field
-    volScalarField rhol
-    (
-        IOobject
-        (
-            "rhol",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );*/
-
     Info<< "\nReading field U\n" << endl;
     volVectorField U
     (
@@ -133,9 +88,3 @@
         fields.add(Y[i]);
     }
     fields.add(h);
-
-    Info<< "Creating radiation model\n" << endl;
-    autoPtr<radiation::radiationModel> radiation
-    (
-        radiation::radiationModel::New(T)
-    );
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/hEqn.H b/applications/solvers/Lagrangian/reactingParcelFoam/hEqn.H
index e4fc046097431ef7335ffd1737382d6c82ec4a89..7909d8c67d2b3962b6da9a75da8637ccd4673ed0 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/hEqn.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/hEqn.H
@@ -6,7 +6,7 @@
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
         DpDt
-     +  reactingParcels.Sh1()
+     +  parcels.Sh()
      +  radiation->Sh(thermo())
     );
 
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/Lagrangian/reactingParcelFoam/pEqn.H
index 5516afe5159f294059599c2f3456b2ed2dab18d8..63370f0cd41ea582ddaa646e9861f7d362293f00 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/pEqn.H
@@ -23,7 +23,7 @@ if (transonic)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rUA, p)
          ==
-            reactingParcels.Srho1()
+            parcels.Srho()
         );
 
         pEqn.solve();
@@ -51,7 +51,7 @@ else
           + fvc::div(phi)
           - fvm::laplacian(rho*rUA, p)
          ==
-            reactingParcels.Srho1()
+            parcels.Srho()
         );
 
         pEqn.solve();
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/Lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 122f7dffb74303ec90d0a375459852692763eb79..eb170cb7c4a189f308976f7a5c6f56f95123b56e 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -23,8 +23,11 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
+    reactingParcelFoam
 
 Description
+    Transient PISO solver for compressible, laminar or turbulent flow with
+    reacting Lagrangian parcels.
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,55 +45,54 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readChemistryProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "createClouds.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createClouds.H"
+    #include "createRadiationModel.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        Info << "Evolving reacting cloud" << endl;
-
-        reactingParcels.evolve();
+        parcels.evolve();
 
-        reactingParcels.info();
+        parcels.info();
 
-#       include "chemistry.H"
-#       include "rhoEqn.H"
+        #include "chemistry.H"
+        #include "rhoEqn.H"
 
         // --- PIMPLE loop
         for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
         {
-#           include "UEqn.H"
-#           include "YEqn.H"
+            #include "UEqn.H"
+            #include "YEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "hEqn.H"
-#               include "pEqn.H"
+                #include "hEqn.H"
+                #include "pEqn.H"
             }
 
             Info<< "T gas min/max   = " << min(T).value() << ", "
@@ -103,7 +105,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-#           include "additionalOutput.H"
+            #include "additionalOutput.H"
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/Lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/Lagrangian/reactingParcelFoam/rhoEqn.H
index d6b129452e7dd07680c66c4e8567e14ccccdbc38..e30e721004d7e251eb23e74182c475149f5bd2b4 100644
--- a/applications/solvers/Lagrangian/reactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/Lagrangian/reactingParcelFoam/rhoEqn.H
@@ -36,7 +36,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        reactingParcels.Srho1()
+        parcels.Srho()
     );
 }