diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
index f1044030bda04737c4d76177ed0d286b3151286e..15f448ff5735ac57a8377d8999c8b6ca2a40c5b0 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2020 ENERCON GmbH
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 
 #include "atmPlantCanopyUSource.H"
 #include "addToRunTimeSelectionTable.H"
+#include "fvmSup.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -99,7 +100,7 @@ void Foam::fv::atmPlantCanopyUSource::addSup
     if (V_ > VSMALL)
     {
         // (SP:Eq. 42)
-        eqn -= (plantCd_*leafAreaDensity_*mag(U))*U;
+        eqn -= fvm::Sp(plantCd_*leafAreaDensity_*mag(U), U);
     }
 }
 
@@ -115,7 +116,7 @@ void Foam::fv::atmPlantCanopyUSource::addSup
 
     if (V_ > VSMALL)
     {
-        eqn -= rho*(plantCd_*leafAreaDensity_*mag(U))*U;
+        eqn -= fvm::Sp(rho*plantCd_*leafAreaDensity_*mag(U), U);
     }
 }
 
@@ -132,7 +133,7 @@ void Foam::fv::atmPlantCanopyUSource::addSup
 
     if (V_ > VSMALL)
     {
-        eqn -= alpha*rho*(plantCd_*leafAreaDensity_*mag(U))*U;
+        eqn -= fvm::Sp(alpha*rho*plantCd_*leafAreaDensity_*mag(U), U);
     }
 }