diff --git a/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.C b/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.C
index 9a2f655624042f151c2dc57c4e4a4ca30c287df2..ad392c506bdef7e6c6b91aabf992118a9d0683e6 100644
--- a/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.C
+++ b/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.C
@@ -146,6 +146,7 @@ Foam::scalarTransport::scalarTransport
     active_(true),
     phiName_("phi"),
     UName_("U"),
+    rhoName_("rho"),
     DT_(0.0),
     userDT_(false),
     resetOnStartUp_(false),
@@ -192,6 +193,7 @@ void Foam::scalarTransport::read(const dictionary& dict)
 
         phiName_ = dict.lookupOrDefault<word>("phiName", "phi");
         UName_ = dict.lookupOrDefault<word>("UName", "U");
+        rhoName_ = dict.lookupOrDefault<word>("rhoName", "rho");
 
         userDT_ = false;
         if (dict.readIfPresent("DT", DT_))
@@ -237,24 +239,59 @@ void Foam::scalarTransport::execute()
         relaxCoeff = mesh_.equationRelaxationFactor(schemeVar);
     }
 
-    // solve
-    for (label i = 0; i <= nCorr_; i++)
+    if (phi.dimensions() == dimMass/dimTime)
     {
-        fvScalarMatrix TEqn
-        (
-            fvm::ddt(T_)
-          + fvm::div(phi, T_, divScheme)
-          - fvm::laplacian(DT, T_, laplacianScheme)
-         ==
-            fvOptions_(T_)
-        );
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        // solve
+        for (label i = 0; i <= nCorr_; i++)
+        {
+            fvScalarMatrix TEqn
+            (
+                fvm::ddt(rho, T_)
+              + fvm::div(phi, T_, divScheme)
+              - fvm::laplacian(DT, T_, laplacianScheme)
+             ==
+                fvOptions_(rho, T_)
+            );
+
+            TEqn.relax(relaxCoeff);
+
+            fvOptions_.constrain(TEqn);
+
+            TEqn.solve(mesh_.solverDict(schemeVar));
+        }
+    }
+    else if (phi.dimensions() == dimVolume/dimTime)
+    {
+        // solve
+        for (label i = 0; i <= nCorr_; i++)
+        {
+            fvScalarMatrix TEqn
+            (
+                fvm::ddt(T_)
+              + fvm::div(phi, T_, divScheme)
+              - fvm::laplacian(DT, T_, laplacianScheme)
+             ==
+                fvOptions_(T_)
+            );
 
-        TEqn.relax(relaxCoeff);
+            TEqn.relax(relaxCoeff);
 
-        fvOptions_.constrain(TEqn);
+            fvOptions_.constrain(TEqn);
 
-        TEqn.solve(mesh_.solverDict(UName_));
+            TEqn.solve(mesh_.solverDict(schemeVar));
+        }
     }
+    else
+    {
+        FatalErrorIn("void Foam::scalarTransport::execute()")
+            << "Incompatible dimensions for phi: " << phi.dimensions() << nl
+            << "Dimensions should be " << dimMass/dimTime << " or "
+            << dimVolume/dimTime << endl;
+    }
+
 
     Info<< endl;
 }
diff --git a/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.H b/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.H
index fdeecdf3604e4ce6d6e5f5c035ed8a8e883abd9b..a58ed2268049379831ceaff4c1042434300d6f52 100644
--- a/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.H
+++ b/src/postProcessing/functionObjects/utilities/scalarTransport/scalarTransport.H
@@ -90,6 +90,9 @@ class scalarTransport
         //- Name of velocity field (optional)
         word UName_;
 
+        //- Name of density field (optional)
+        word rhoName_;
+
         //- Diffusion coefficient (optional)
         scalar DT_;