diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index 6bec827a50946bb2e08ba0c79755ae13e84dd9df..8f4e9b8725911054742ab6b1c132b2689e78708f 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -172,7 +172,7 @@ updateCoeffs()
     radiativeIntensityRay& ray =
         const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
 
-    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
+    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());
 
     forAll(Iw, faceI)
     {
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
index 6dfffad991beb131b1ae5774b1aead76d426d276..d97689cb21e35856ea3026bff881f57d5c5b4c64 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -167,7 +167,7 @@ updateCoeffs()
     radiativeIntensityRay& ray =
         const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
 
-    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
+    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());
 
     const scalarField Eb =
         dom.blackBody().bLambda(lambdaId).boundaryField()[patchI];
diff --git a/src/thermophysicalModels/radiation/radiationModel/P1/P1.C b/src/thermophysicalModels/radiation/radiationModel/P1/P1.C
index b9fc90f74c08d836aa0573927614589389901dc6..74cf2165360c5702c563246f0db3b8a93d23feff 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_*physicoChemical::sigma*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 5300d3718ce94bdf88e19109b2ce060d6cc5d1b3..7c03d5313e0d1d07977944d253fd856be0e1bf0a 100644
--- a/src/thermophysicalModels/radiation/radiationModel/P1/P1.H
+++ b/src/thermophysicalModels/radiation/radiationModel/P1/P1.H
@@ -63,6 +63,9 @@ class P1
         //- Incident radiation / [W/m2]
         volScalarField G_;
 
+        //- Total radiative heat flux [W/m2]
+        volScalarField Qr_;
+
         //- Absorption coefficient
         volScalarField a_;