diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
index 11053f31a952b0aabfcb4e228da801f80f0233dd..502938c53c027fcd59c4ae9bc1e7b7c7c84d52fd 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../rhoPorousMRFPimpleFoam \
     -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
index ce7e7a2d85c57b9304c970f0a7019cc8a110c9e7..68bcd9ef06f4228e90f1fa9e852783691e64630b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
@@ -10,5 +10,4 @@ EXE_INC = \
 LIB_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
-    -lcompressibleMultiphaseEulerianInterfacialModels \
     -lfiniteVolume
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
index e42c52451fce119c4084e9f8655da650a0700f8d..bc4a36351798496a553e918208b05dc85bc4271e 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
@@ -64,6 +64,7 @@ initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
 initialPointsMethod/faceCentredCubic/faceCentredCubic.C
 initialPointsMethod/pointFile/pointFile.C
 initialPointsMethod/autoDensity/autoDensity.C
+initialPointsMethod/rayShooting/rayShooting.C
 
 relaxationModel/relaxationModel/relaxationModel.C
 relaxationModel/adaptiveLinear/adaptiveLinear.C
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 6f137249c1c4b8ecaa5680b60aa19663b011e9f7..e580febdd8f90d45631759cbffd25b4310270d49 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1249,7 +1249,7 @@ void Foam::conformalVoronoiMesh::move()
                     if
                     (
                         (
-                            (vA->internalPoint() && vB->internalPoint())
+                            (vA->internalPoint() || vB->internalPoint())
                          && (!vA->referred() || !vB->referred())
 //                         ||
 //                            (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 20a8b278ebbee234255272bc079e4ad0fc4430ce..0360cf0d8316623e4dc7c0f4c29f9095d480ef62 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -670,6 +670,32 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 continue;
             }
 
+            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
+
+            if (!surface.hasVolumeType())
+            {
+                pointField sample(1, samplePts[i]);
+                scalarField nearestDistSqr(1, GREAT);
+                List<pointIndexHit> info;
+
+                surface.findNearest(sample, nearestDistSqr, info);
+
+                vector hitDir = info[0].rawPoint() - samplePts[i];
+                hitDir /= mag(hitDir) + SMALL;
+
+                if
+                (
+                    findSurfaceAnyIntersection
+                    (
+                        samplePts[i],
+                        info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir)
+                    )
+                )
+                {
+                    continue;
+                }
+            }
+
             if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
             {
                 if
@@ -694,44 +720,6 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                     break;
                 }
             }
-//            else
-//            {
-//                // Surface volume type is unknown
-//                Info<< "UNKNOWN" << endl;
-//                // Get nearest face normal
-//
-//                pointField sample(1, samplePts[i]);
-//                scalarField nearestDistSqr(1, GREAT);
-//                List<pointIndexHit> info;
-//                vectorField norms(1);
-//
-//                surface.findNearest(sample, nearestDistSqr, info);
-//                surface.getNormal(info, norms);
-//
-//                vector fN = norms[0];
-//                fN /= mag(fN);
-//
-//                vector hitDir = info[0].rawPoint() - samplePts[i];
-//                hitDir /= mag(hitDir);
-//
-//                if ((fN & hitDir) < 0)
-//                {
-//                    // Point is OUTSIDE
-//
-//                    if
-//                    (
-//                        normalVolumeTypes_[regionI]
-//                     == extendedFeatureEdgeMesh::OUTSIDE
-//                    )
-//                    {
-//                    }
-//                    else
-//                    {
-//                        insidePoint[i] = false;
-//                        break;
-//                    }
-//                }
-//            }
         }
     }
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
new file mode 100644
index 0000000000000000000000000000000000000000..ebecad88e168f7377c61993ef4942d1d7dc7d61e
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rayShooting.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(rayShooting, 0);
+addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+rayShooting::rayShooting
+(
+    const dictionary& initialPointsDict,
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
+)
+:
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
+    randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
+    randomPerturbationCoeff_
+    (
+        readScalar(detailsDict().lookup("randomPerturbationCoeff"))
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+List<Vb::Point> rayShooting::initialPoints() const
+{
+    // Loop over surface faces
+    const searchableSurfaces& surfaces = geometryToConformTo().geometry();
+    const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
+
+    const scalar maxRayLength = surfaces.bounds().mag();
+
+    // Initialise points list
+    label initialPointsSize = 0;
+    forAll(surfaces, surfI)
+    {
+        initialPointsSize += surfaces[surfI].size();
+    }
+
+    DynamicList<Vb::Point> initialPoints(initialPointsSize);
+
+    forAll(surfacesToConformTo, surfI)
+    {
+        const searchableSurface& s = surfaces[surfacesToConformTo[surfI]];
+
+        tmp<pointField> faceCentresTmp(s.coordinates());
+        const pointField& faceCentres = faceCentresTmp();
+
+        Info<< "    Shoot rays from " << s.name() << nl
+            << "    nRays = " << faceCentres.size() << endl;
+
+
+        forAll(faceCentres, fcI)
+        {
+            const Foam::point& fC = faceCentres[fcI];
+
+            if
+            (
+                Pstream::parRun()
+             && !decomposition().positionOnThisProcessor(fC)
+            )
+            {
+                continue;
+            }
+
+            const scalar pert =
+                randomPerturbationCoeff_
+               *cellShapeControls().cellSize(fC);
+
+            pointIndexHit surfHitStart;
+            label hitSurfaceStart;
+
+            // Face centres should be on the surface so search distance can be
+            // small
+            geometryToConformTo().findSurfaceNearest
+            (
+                 fC,
+                 sqr(pert),
+                 surfHitStart,
+                 hitSurfaceStart
+            );
+
+            vectorField normStart(1, vector::min);
+            geometryToConformTo().getNormal
+            (
+                hitSurfaceStart,
+                List<pointIndexHit>(1, surfHitStart),
+                normStart
+            );
+
+            pointIndexHit surfHitEnd;
+            label hitSurfaceEnd;
+
+            geometryToConformTo().findSurfaceNearestIntersection
+            (
+                fC - normStart[0]*pert,
+                fC - normStart[0]*maxRayLength,
+                surfHitEnd,
+                hitSurfaceEnd
+            );
+
+            if (surfHitEnd.hit())
+            {
+                vectorField normEnd(1, vector::min);
+                geometryToConformTo().getNormal
+                (
+                    hitSurfaceEnd,
+                    List<pointIndexHit>(1, surfHitEnd),
+                    normEnd
+                );
+
+                if ((normStart[0] & normEnd[0]) < 0)
+                {
+                    line<point, point> l(fC, surfHitEnd.hitPoint());
+
+                    if (Pstream::parRun())
+                    {
+                        // Clip the line in parallel
+                        pointIndexHit procIntersection =
+                            decomposition().findLine
+                            (
+                                l.start(),
+                                l.end()
+                            );
+
+                        if (procIntersection.hit())
+                        {
+                            l =
+                                line<point, point>
+                                (
+                                    l.start(),
+                                    procIntersection.hitPoint()
+                                );
+                        }
+                    }
+
+                    Foam::point midPoint(l.centre());
+
+                    const scalar minDistFromSurfaceSqr =
+                        minimumSurfaceDistanceCoeffSqr_
+                       *sqr(cellShapeControls().cellSize(midPoint));
+
+                    if (randomiseInitialGrid_)
+                    {
+                        midPoint.x() += pert*(rndGen().scalar01() - 0.5);
+                        midPoint.y() += pert*(rndGen().scalar01() - 0.5);
+                        midPoint.z() += pert*(rndGen().scalar01() - 0.5);
+                    }
+
+                    if
+                    (
+                        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
+                     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
+                    )
+                    {
+                        initialPoints.append(toPoint(midPoint));
+                    }
+                }
+            }
+        }
+    }
+
+    return initialPoints.shrink();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
new file mode 100644
index 0000000000000000000000000000000000000000..e66df62b7c0f9c6d6632469a31874472274fbd2f
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rayShooting
+
+Description
+
+SourceFiles
+    rayShooting.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rayShooting_H
+#define rayShooting_H
+
+#include "initialPointsMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class rayShooting Declaration
+\*---------------------------------------------------------------------------*/
+
+class rayShooting
+:
+    public initialPointsMethod
+{
+
+private:
+
+    // Private data
+
+        //- Should the initial positions be randomised
+        Switch randomiseInitialGrid_;
+
+        //- Randomise the initial positions by fraction of the initialCellSize_
+        scalar randomPerturbationCoeff_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("rayShooting");
+
+    // Constructors
+
+        //- Construct from components
+        rayShooting
+        (
+            const dictionary& initialPointsDict,
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
+        );
+
+
+    //- Destructor
+    virtual ~rayShooting()
+    {}
+
+
+    // Member Functions
+
+        //- Return the initial points for the conformalVoronoiMesh
+        virtual List<Vb::Point> initialPoints() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 10c8a4b4ac440af69bb2e30ead371a0a0c1bd4c4..36ff9f4cadf50be45bf5aa585d286495e496198e 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -61,6 +61,9 @@ void reactingOneDim::readReactingOneDimControls()
 
     coeffs().lookup("gasHSource") >> gasHSource_;
     coeffs().lookup("QrHSource") >> QrHSource_;
+    useChemistrySolvers_ =
+        coeffs().lookupOrDefault<bool>("useChemistrySolvers", true);
+
 }
 
 
@@ -321,6 +324,8 @@ void reactingOneDim::solveEnergy()
     (
         fvm::ddt(rho_, h_)
       - fvm::laplacian(alpha, h_)
+      + fvc::laplacian(alpha, h_)
+      - fvc::laplacian(kappa(), T())
      ==
         chemistrySh_
       - fvm::Sp(solidChemistry_->RRg(), h_)
@@ -462,7 +467,8 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     totalGasMassFlux_(0.0),
     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
     gasHSource_(false),
-    QrHSource_(false)
+    QrHSource_(false),
+    useChemistrySolvers_(true)
 {
     if (active_)
     {
@@ -560,7 +566,8 @@ reactingOneDim::reactingOneDim
     totalGasMassFlux_(0.0),
     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
     gasHSource_(false),
-    QrHSource_(false)
+    QrHSource_(false),
+    useChemistrySolvers_(true)
 {
     if (active_)
     {
@@ -681,11 +688,18 @@ void reactingOneDim::evolveRegion()
 {
     Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
 
-    solidChemistry_->solve
-    (
-        time().value() - time().deltaTValue(),
-        time().deltaTValue()
-    );
+    if (useChemistrySolvers_)
+    {
+        solidChemistry_->solve
+        (
+            time().value() - time().deltaTValue(),
+            time().deltaTValue()
+        );
+    }
+    else
+    {
+        solidChemistry_->calculate();
+    }
 
     solveContinuity();
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
index 4a92fd69677304bb2b4d4902e5942b14dabcef6e..a5950a9804b6f52532dd707603c4d04b4eb328e5 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
@@ -154,6 +154,9 @@ protected:
             //- Add in depth radiation source term
             bool  QrHSource_;
 
+            //- Use chemistry solvers (ode or sequential)
+            bool useChemistrySolvers_;
+
 
     // Protected member functions
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index 67704bb9df1b32bebb5e38357590260e4b6c5098..75ea14eb676ac1a4f84508c8108b7643b823c30a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -140,6 +140,12 @@ public:
                     const label i
                 ) const = 0;
 
+                //- Return access to chemical source terms [kg/m3/s]
+                virtual DimensionedField<scalar, volMesh>& RR
+                (
+                    const label i
+                ) = 0;
+
 
             // Chemistry solution
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 7b3d1c663a3314a5cff53df22d3a523f6d4823d8..475fc7c8fd77f89ec249d4878d651ba05d456413 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -211,6 +211,12 @@ public:
                 const label i
             ) const;
 
+            //- Return non const access to chemical source terms [kg/m3/s]
+            virtual DimensionedField<scalar, volMesh>& RR
+            (
+                const label i
+            );
+
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
             virtual scalar solve(const scalar t0, const scalar deltaT);
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
index d128443fe9e4e9ea9fcd831ca20be545f085c84f..c009a11e4695e0e4f3188092129b1475d3857d92 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,5 +78,15 @@ Foam::chemistryModel<CompType, ThermoType>::RR
     return RR_[i];
 }
 
+template<class CompType, class ThermoType>
+Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::chemistryModel<CompType, ThermoType>::RR
+(
+    const label i
+)
+{
+    return RR_[i];
+}
+
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
index cfff71b9edc4a8393dba6b582ecf546866a3baa8..f94db5946b3bf7bbf88b507e670dd7e7730bbbd9 100644
--- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
@@ -53,4 +53,35 @@ Foam::basicSolidChemistryModel::~basicSolidChemistryModel()
 {}
 
 
+const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicSolidChemistryModel::RR(const label i) const
+{
+    notImplemented
+    (
+        "const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&"
+        "basicSolidChemistryModel::RR(const label)"
+    );
+    return (DimensionedField<scalar, volMesh>::null());
+}
+
+
+Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicSolidChemistryModel::RR(const label i)
+{
+    notImplemented
+    (
+        "Foam::DimensionedField<Foam::scalar, Foam::volMesh>&"
+        "basicSolidChemistryModel::RR(const label)"
+    );
+
+    return dynamic_cast<DimensionedField<scalar, volMesh>&>
+    (
+        const_cast<DimensionedField<scalar, volMesh>& >
+        (
+            DimensionedField<scalar, volMesh>::null()
+        )
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
index 2134bffca1a1df8bdfa429f7223f5329c6cf77e8..50d6a2e96bec1dcbe88d4a8ecaafcb7d571ab0de 100644
--- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
@@ -149,6 +149,15 @@ public:
 
         //- Calculates the reaction rates
         virtual void calculate() = 0;
+
+        //- Return const access to the total source terms
+        virtual const DimensionedField<scalar, volMesh>& RR
+        (
+            const label i
+        ) const;
+
+        //- Return non-const access to the total source terms
+        virtual DimensionedField<scalar, volMesh>& RR(const label i);
 };
 
 
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
index e301a1b85e24a08666f369baa04f4bc1b4dbda8b..4962bf716b6b5a34079b44e0db1fdee361f30913 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
@@ -535,14 +535,21 @@ updateConcsInReactionI
         c[si] = max(0.0, c[si]);
     }
 
+    scalar sr = 0.0;
     forAll(R.rhs(), s)
     {
         label si = R.rhs()[s].index;
         const scalar rhoR = this->solidThermo_[si].rho(p, T);
-        const scalar sr = rhoR/rhoL;
+        sr = rhoR/rhoL;
         c[si] += dt*sr*omeg;
         c[si] = max(0.0, c[si]);
     }
+
+    forAll(R.grhs(), g)
+    {
+        label gi = R.grhs()[g].index;
+        c[gi + this->nSolids_] += (1.0 - sr)*omeg*dt;
+    }
 }
 
 
@@ -561,24 +568,11 @@ updateRRInReactionI
     simpleMatrix<scalar>& RR
 ) const
 {
-    const Reaction<SolidThermo>& R = this->reactions_[index];
-    scalar rhoL = 0.0;
-    forAll(R.lhs(), s)
-    {
-        label si = R.lhs()[s].index;
-        rhoL = this->solidThermo_[si].rho(p, T);
-        RR[si][rRef] -= pr*corr;
-        RR[si][lRef] += pf*corr;
-    }
-
-    forAll(R.rhs(), s)
-    {
-        label si = R.rhs()[s].index;
-        const scalar rhoR = this->solidThermo_[si].rho(p, T);
-        const scalar sr = rhoR/rhoL;
-        RR[si][lRef] -= sr*pf*corr;
-        RR[si][rRef] += sr*pr*corr;
-    }
+    notImplemented
+    (
+        "void Foam::pyrolysisChemistryModel<CompType, SolidThermo,GasThermo>::"
+        "updateRRInReactionI"
+    );
 }
 
 
@@ -666,7 +660,6 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
 
                 scalar newCp = 0.0;
                 scalar newhi = 0.0;
-                scalar invRho = 0.0;
                 scalarList dcdt = (c - c0)/dt;
 
                 for (label i=0; i<this->nSolids_; i++)
@@ -675,7 +668,6 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
                     scalar Yi = c[i]/cTot;
                     newCp += Yi*this->solidThermo_[i].Cp(pi, Ti);
                     newhi -= dYi*this->solidThermo_[i].Hc();
-                    invRho += Yi/this->solidThermo_[i].rho(pi, Ti);
                 }
 
                 scalar dTi = (newhi/newCp)*dt;
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
index 2d4a11cc90fea0485b625f1e0bfd4657d52151ad..7ae76c5233383221b14539d5b092439617b8c299 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
@@ -137,8 +137,9 @@ public:
             const bool updateC0 = false
         ) const;
 
-        //- Return the reaction rate for reaction r and the reference
-        //  species and charateristic times
+        //- Return the reaction rate for reaction r
+        // NOTE: Currently does not calculate reference specie and
+        // characteristic times (pf, cf,l Ref, etc.)
         virtual scalar omega
         (
             const Reaction<SolidThermo>& r,
@@ -153,8 +154,10 @@ public:
             label& rRef
         ) const;
 
-        //- Return the reaction rate for iReaction and the reference
-        //  species and charateristic times
+
+        //- Return the reaction rate for iReaction
+        // NOTE: Currently does not calculate reference specie and
+        // characteristic times (pf, cf,l Ref, etc.)
         virtual scalar omegaI
         (
             label iReaction,
@@ -169,6 +172,7 @@ public:
             label& rRef
         ) const;
 
+
         //- Calculates the reaction rates
         virtual void calculate();
 
@@ -186,6 +190,7 @@ public:
 
 
         //- Update matrix RR for reaction i. Used by EulerImplicit
+        // (Not implemented)
         virtual void updateRRInReactionI
         (
             const label i,
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 03076d79ceeb9a330601066d22d1d500c826d429..bf67f1ca4d27c9662c29369d68d724b1aa854366 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,9 @@ class solidChemistryModel
 {
     // Private Member Functions
 
+        //- Disallow copy constructor
+        solidChemistryModel(const solidChemistryModel&);
+
         //- Disallow default bitwise assignment
         void operator=(const solidChemistryModel&);
 
@@ -151,6 +154,7 @@ public:
             label& rRef
         ) const = 0;
 
+
         //- Return the reaction rate for iReaction and the reference
         //  species and charateristic times
         virtual scalar omegaI
@@ -167,6 +171,10 @@ public:
             label& rRef
         ) const = 0;
 
+        //- Calculates the reaction rates
+        virtual void calculate() = 0;
+
+
         //- Update concentrations in reaction i given dt and reaction rate
         // omega used by sequential solver
         virtual void updateConcsInReactionI
@@ -194,11 +202,8 @@ public:
             simpleMatrix<scalar>& RR
         ) const = 0;
 
-        //- Calculates the reaction rates
-        virtual void calculate() = 0;
-
 
-        // Chemistry model functions
+        // Solid Chemistry model functions
 
             //- Return const access to the chemical source terms for solids
             inline const DimensionedField<scalar, volMesh>& RRs
@@ -209,13 +214,6 @@ public:
             //- Return total solid source term
             inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
 
-            //- Return const access to the total source terms
-            inline const DimensionedField<scalar, volMesh>& RR
-            (
-                const label i
-            ) const;
-
-
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
             virtual scalar solve(const scalar t0, const scalar deltaT) = 0;
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
index 92b461e11ea93273757627dab8ed36a55a990cf7..f2d0dad79ab58e027031889248da0c99f4704f09 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,16 +95,4 @@ Foam::solidChemistryModel<CompType, SolidThermo>::RRs() const
 }
 
 
-template<class CompType, class SolidThermo>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::solidChemistryModel<CompType, SolidThermo>::RR
-(
-    const label i
-) const
-{
-    notImplemented("solidChemistryModel::RR(const label)");
-    return (DimensionedField<scalar, volMesh>::null());
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
index 12f842a054399d2a2d5197c4bc61ccc50e1bfe20..447020e2f7ded4490ab5b6b82d2f45d7846063eb 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
@@ -52,9 +52,8 @@ namespace Foam
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
         SS##Schem##Comp##SThermo##GThermo,                                    \
-        (#SS"<" + word(Schem::typeName_())                                    \
-      + "<"#Comp"," + SThermo::typeName()                                     \
-      + "," + GThermo::typeName() + ">>").c_str(),                            \
+        (#SS"<"#Schem"<"#Comp"," + SThermo::typeName() + ","                  \
+      + GThermo::typeName() + ">>").c_str(),                                  \
         0                                                                     \
     );                                                                        \
                                                                               \
@@ -77,14 +76,6 @@ namespace Foam
         GThermo                                                               \
     );                                                                        \
                                                                               \
-    makeSolidChemistrySolverType                                              \
-    (                                                                         \
-        EulerImplicit,                                                        \
-        SolidChem,                                                            \
-        Comp,                                                                 \
-        SThermo,                                                              \
-        GThermo                                                               \
-    );                                                                        \
                                                                               \
     makeSolidChemistrySolverType                                              \
     (                                                                         \
@@ -104,7 +95,6 @@ namespace Foam
         GThermo                                                               \
     );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
index 5193a4b9deef4e3d95b10f06d9231f625fb46e72..07b6ffb050242c50661b5adfcff50fe1e14db104 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
+++ b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeSolidThermo_H
-#define makesolidThermo_H
+#define makeSolidThermo_H
 
 #include "addToRunTimeSelectionTable.H"