diff --git a/applications/solvers/dsmc/dsmcFoam/createFields.H b/applications/solvers/dsmc/dsmcFoam/createFields.H
index 74c1d5c055f64033f0decfb1bc486609e10fe5f8..d024bd2017ae346c53e9adbc86a775a9ba496ba5 100644
--- a/applications/solvers/dsmc/dsmcFoam/createFields.H
+++ b/applications/solvers/dsmc/dsmcFoam/createFields.H
@@ -1,10 +1,10 @@
 
-    Info<< nl << "Reading field T" << endl;
-    volScalarField T
+    Info<< nl << "Reading field boundaryT" << endl;
+    volScalarField boundaryT
     (
         IOobject
         (
-            "T",
+            "boundaryT",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -13,12 +13,54 @@
         mesh
     );
 
-    Info<< nl << "Reading field U" << endl;
-    volVectorField U
+    Info<< nl << "Reading field boundaryU" << endl;
+    volVectorField boundaryU
     (
         IOobject
         (
-            "U",
+            "boundaryU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoN (number density)" << endl;
+    volScalarField rhoN
+    (
+        IOobject
+        (
+            "rhoN",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoM (mass density)" << endl;
+    volScalarField rhoM
+    (
+        IOobject
+        (
+            "rhoM",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl;
+    volScalarField dsmcRhoN
+    (
+        IOobject
+        (
+            "dsmcRhoN",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -41,12 +83,14 @@
         mesh
     );
 
-    Info<< nl << "Reading field rhoN" << endl;
-    volScalarField rhoN
+    Info<< nl << "Reading field linearKE (linear kinetic energy density)"
+        << endl;
+
+    volScalarField linearKE
     (
         IOobject
         (
-            "rhoN",
+            "linearKE",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -55,12 +99,12 @@
         mesh
     );
 
-    Info<< nl << "Reading field rhoM" << endl;
-    volScalarField rhoM
+    Info<< nl << "Reading field internalE (internal energy density)" << endl;
+    volScalarField internalE
     (
         IOobject
         (
-            "rhoM",
+            "internalE",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -69,12 +113,14 @@
         mesh
     );
 
-    Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl;
-    volScalarField dsmcRhoN
+    Info<< nl << "Reading field iDof (internal degree of freedom density)"
+        << endl;
+
+    volScalarField iDof
     (
         IOobject
         (
-            "dsmcRhoN",
+            "iDof",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -83,10 +129,6 @@
         mesh
     );
 
-    Info<< nl
-        << "Fields required to record solid surface forces and heat flux:"
-        << endl;
-
     Info<< nl << "Reading field q (surface heat transfer)" << endl;
     volScalarField q
     (
@@ -101,7 +143,7 @@
         mesh
     );
 
-    Info<< nl << "Reading field fD (surface force density" << endl;
+    Info<< nl << "Reading field fD (surface force density)" << endl;
     volVectorField fD
     (
         IOobject
@@ -115,6 +157,6 @@
         mesh
     );
 
-    Info<< "Constructing dsmcCloud " << endl;
+    Info<< nl << "Constructing dsmcCloud " << endl;
 
-dsmcCloud dsmc("dsmc", T, U);
+    dsmcCloud dsmc("dsmc", boundaryT, boundaryU);
diff --git a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
index a9d30e9ff038c7aa38e401b6b2ff695bb1563269..920bdc5f204c98beecab9853adaa12ec1ff9f163 100644
--- a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
+++ b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
@@ -66,6 +66,12 @@ int main(int argc, char *argv[])
 
         momentum = dsmc.momentum();
 
+        linearKE = dsmc.linearKE();
+
+        internalE = dsmc.internalE();
+
+        iDof = dsmc.iDof();
+
         q = dsmc.q();
 
         fD = dsmc.fD();
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/files b/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/files
new file mode 100755
index 0000000000000000000000000000000000000000..c667194d886e9844d31d125b2fa8502b6bb04f82
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/files
@@ -0,0 +1,3 @@
+dsmcFields.C
+
+EXE = $(FOAM_APPBIN)/dsmcFields
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/options b/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/options
new file mode 100755
index 0000000000000000000000000000000000000000..4eb433b999c3b4d223c83a2d62216a3936415b7c
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/postProcessing/postCalc \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    $(FOAM_LIBBIN)/postCalc.o \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -ldsmc
+
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFields/dsmcFields.C b/applications/utilities/postProcessing/miscellaneous/dsmcFields/dsmcFields.C
new file mode 100644
index 0000000000000000000000000000000000000000..6d7450835553226b243e6b7499b3058aff5d20c5
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFields/dsmcFields.C
@@ -0,0 +1,231 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    dsmcFields
+
+Description
+    Calculate intensive fields (U and T) from averaged extensive fields from a
+    DSMC calculation.
+
+\*---------------------------------------------------------------------------*/
+
+#include "calc.H"
+#include "fvc.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
+{
+    bool writeResults = !args.options().found("noWrite");
+
+    IOobject rhoNMeanheader
+    (
+        "rhoNMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject rhoMMeanheader
+    (
+        "rhoMMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject momentumMeanheader
+    (
+        "momentumMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject linearKEMeanheader
+    (
+        "linearKEMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject internalEMeanheader
+    (
+        "internalEMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject iDofMeanheader
+    (
+        "iDofMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    if (!rhoNMeanheader.headerOk())
+    {
+        Info<< "    No rhoNMean" << endl;
+    }
+    else if (!rhoMMeanheader.headerOk())
+    {
+        Info<< "    No rhoMMean" << endl;
+    }
+    else if (!momentumMeanheader.headerOk())
+    {
+        Info<< "    No momentumMean" << endl;
+    }
+    else if (!linearKEMeanheader.headerOk())
+    {
+        Info<< "    No linearKEMean" << endl;
+    }
+    else if (!internalEMeanheader.headerOk())
+    {
+        Info<< "    No internalEMean" << endl;
+    }
+    else if (!iDofMeanheader.headerOk())
+    {
+        Info<< "    No iDofMean" << endl;
+    }
+    else
+    {
+        Info<< "Reading field rhoNMean" << endl;
+        volScalarField rhoNMean(rhoNMeanheader, mesh);
+
+        Info<< "Reading field rhoMMean" << endl;
+        volScalarField rhoMMean(rhoMMeanheader, mesh);
+
+        Info<< "Reading field momentumMean" << endl;
+        volVectorField momentumMean(momentumMeanheader, mesh);
+
+        Info<< "Reading field linearKEMean" << endl;
+        volScalarField linearKEMean(linearKEMeanheader, mesh);
+
+        Info<< "Reading field internalEMean" << endl;
+        volScalarField internalEMean(internalEMeanheader, mesh);
+
+        Info<< "Reading field iDofMean" << endl;
+        volScalarField iDofMean(iDofMeanheader, mesh);
+
+        // Check if there are any zero values in the density fields
+
+        if (min(rhoNMean).value() > VSMALL)
+        {
+            // TODO Sort out boundary field values if required
+
+            Info<< nl << "Calculating UMean field." << endl;
+
+            volVectorField UMean
+            (
+                IOobject
+                (
+                    "UMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                momentumMean/rhoMMean
+            );
+
+            Info<< nl << "Calculating TMean fields." << endl;
+
+            volScalarField translationalTMean
+            (
+                IOobject
+                (
+                    "translationalTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(3.0*dsmcCloud::kb*rhoNMean)
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
+            );
+
+            volScalarField internalTMean
+            (
+                IOobject
+                (
+                    "internalTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*iDofMean)*internalEMean
+            );
+
+            volScalarField overallTMean
+            (
+                IOobject
+                (
+                    "overallTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
+            );
+
+            Info<< nl << "magUMean max/min : "
+                << max(mag(UMean)).value() << " "
+                << min(mag(UMean)).value() << endl;
+
+            Info<< nl << "translationalTMean max/min : "
+                << max(translationalTMean).value() << " "
+                << min(translationalTMean).value() << endl;
+
+            Info<< nl << "internalTMean max/min : "
+                << max(internalTMean).value() << " "
+                << min(internalTMean).value() << endl;
+
+            Info<< nl << "overallTMean max/min : "
+                << max(overallTMean).value() << " "
+                << min(overallTMean).value() << endl;
+
+            if (writeResults)
+            {
+                UMean.write();
+                translationalTMean.write();
+                internalTMean.write();
+                overallTMean.write();
+            }
+        }
+        else
+        {
+            Info<< "Small or negative value (" << min(rhoNMean)
+                << ") found in rhoNMean field. "
+                << "Not calculating fields to avoid division by zero "
+                << "or invalid results."
+                << endl;
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index d5dda1f680345e3793d894823753fa9563ad59da..44ce59af501888027737060667f81be25a817734 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -568,12 +568,11 @@ Foam::DsmcCloud<ParcelType>::iDof() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimless, 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
         )
     );
 
     scalarField& iDof = tiDof().internalField();
-    scalarField nDsmcParticles(iDof.size(),0);
 
     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
     {
@@ -581,12 +580,9 @@ Foam::DsmcCloud<ParcelType>::iDof() const
         const label cellI = p.cell();
 
         iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
-
-        // Avoiding divide by zero for nDsmcParticles where a cell is empty
-        nDsmcParticles[cellI] = max(1, ++nDsmcParticles[cellI]);
     }
 
-    iDof /= nDsmcParticles;
+    iDof *= nParticle_/mesh().cellVolumes();
 
     return tiDof;
 }