diff --git a/src/thermophysicalModels/radiation/radiationModel/P1/P1.C b/src/thermophysicalModels/radiation/radiationModel/P1/P1.C
index 0ff2685f5..7c3261bb7 100644
--- a/src/thermophysicalModels/radiation/radiationModel/P1/P1.C
+++ b/src/thermophysicalModels/radiation/radiationModel/P1/P1.C
@@ -68,6 +68,19 @@ Foam::radiation::P1::P1(const volScalarField& T)
         ),
         mesh_
     ),
+    Qr_
+    (
+        IOobject
+        (
+            "Qr",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+    ),
     a_
     (
         IOobject
@@ -162,6 +175,13 @@ void Foam::radiation::P1::calculate()
      ==
       - 4.0*(e_*radiation::sigmaSB*pow4(T_) + E_)
     );
+
+    // Calculate radiative heat flux on boundaries.
+    forAll(mesh_.boundaryMesh(), patchI)
+    {
+        Qr_.boundaryField()[patchI] =
+            -gamma.boundaryField()[patchI]*G_.boundaryField()[patchI].snGrad();
+    }
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/P1/P1.H b/src/thermophysicalModels/radiation/radiationModel/P1/P1.H
index 46b9398a3..6fc6ae8cd 100644
--- a/src/thermophysicalModels/radiation/radiationModel/P1/P1.H
+++ b/src/thermophysicalModels/radiation/radiationModel/P1/P1.H
@@ -64,6 +64,9 @@ class P1
         //- Incident radiation / [W/m2]
         volScalarField G_;
 
+        //- Total radiative heat flux [W/m2]
+        volScalarField Qr_;
+
         //- Absorption coefficient
         volScalarField a_;