diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index d99c35dcfb..f537f6283a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -26,6 +26,8 @@ License
 
 #include "physicoChemicalConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -229,7 +231,7 @@ Foam::ThermoCloud<ParcelType>::Ep() const
             Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
         }
 
-        Ep *= epsilon*constant::physicoChemical::sigma.value()/V;
+        Ep *= epsilon*physicoChemical::sigma.value()/V;
     }
 
     return tEp;
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index cb301d8210..a895e9df10 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -102,7 +102,7 @@ void Foam::KinematicParcel<ParcelType>::calc
     const scalar mass0 = mass();
 
     // Reynolds number
-    const scalar Re = this->Re(U0, d0, muc_);
+    const scalar Re = this->Re(U0, d0, rhoc_, muc_);
 
 
     // Sources
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 29979d106d..973c97a2d8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -389,12 +389,13 @@ public:
             //- Surface area for given diameter
             inline scalar areaS(const scalar d) const;
 
-            //- Reynolds number - particle properties input
+            //- Reynolds number
             inline scalar Re
             (
-                const vector& U,
-                const scalar d,
-                const scalar mu
+                const vector& U,        // particle velocity
+                const scalar d,         // particle diameter
+                const scalar rhoc,      // carrier density
+                const scalar muc        // carrier dynamic viscosity
             ) const;
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index d847ac38e0..84ccbdec5e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -26,6 +26,8 @@ License
 
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template <class ParcelType>
@@ -351,7 +353,7 @@ template <class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::volume(const scalar d) const
 {
-    return constant::math::pi/6.0*pow3(d);
+    return math::pi/6.0*pow3(d);
 }
 
 
@@ -381,7 +383,7 @@ template <class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::areaS(const scalar d) const
 {
-    return constant::math::pi*d*d;
+    return math::pi*d*d;
 }
 
 
@@ -391,10 +393,11 @@ Foam::KinematicParcel<ParcelType>::Re
 (
     const vector& U,
     const scalar d,
-    const scalar mu
+    const scalar rhoc,
+    const scalar muc
 ) const
 {
-    return rhoc_*mag(U - Uc_)*d/mu;
+    return rhoc*mag(U - Uc_)*d/muc;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 31d8e9b2d6..c035ecf2dc 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -27,6 +27,8 @@ License
 #include "ReactingMultiphaseParcel.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -214,7 +216,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
 
     // Reynolds number
-    scalar Re = this->Re(U0, d0, mus);
+    scalar Re = this->Re(U0, d0, rhos, mus);
 
 
     // Sources
@@ -495,7 +497,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         }
         else
         {
-            this->d_ = cbrt(mass1/this->rho_*6.0/constant::math::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/math::pi);
         }
     }
 }
@@ -642,7 +644,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     const scalar hRetentionCoeffMod =
         (1.0 - xsi*xsi)*td.constProps().hRetentionCoeff();
 
-    Sh += hRetentionCoeffMod *hReaction/dt;
+    Sh += hRetentionCoeffMod*hReaction/dt;
 
     dhsTrans += (1.0 - hRetentionCoeffMod)*hReaction;
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index fdc4d91a28..4ec58ce79e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -159,7 +159,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
     rhos = 0;
     mus = 0;
     kappa = 0;
-    scalar cps = 0.0;
+    scalar cps = 0;
     scalar sumYiSqrtW = 0;
     scalar sumYiCbrtW = 0;
 
@@ -233,11 +233,10 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Calc surface values
     // ~~~~~~~~~~~~~~~~~~~
     scalar Ts, rhos, mus, Pr, kappa;
-    ThermoParcel<ParcelType>::
-        calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
+    this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
 
     // Reynolds number
-    scalar Re = this->Re(U0, d0, mus);
+    scalar Re = this->Re(U0, d0, rhos, mus);
 
 
     // Sources
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 94b71f7d1a..7851841a74 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -133,7 +133,7 @@ void Foam::ThermoParcel<ParcelType>::calc
     calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
 
     // Reynolds number
-    scalar Re = this->Re(U0, d0, mus);
+    scalar Re = this->Re(U0, d0, rhos, mus);
 
 
     // Sources
@@ -231,37 +231,26 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     // Calc heat transfer coefficient
     scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
 
-    const scalar As = this->areaS(d);
-
     if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
     {
         return  T + dt*Sh/(this->volume(d)*rho*cp);
     }
 
-    scalar ap;
-    scalar bp;
-
+    const scalar As = this->areaS(d);
+    scalar ap = Tc_ + Sh/As/htc;
+    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
     if (td.cloud().radiation())
     {
         const scalarField& G =
             td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
+        const scalar Gc = G[cellI];
         const scalar sigma = physicoChemical::sigma.value();
         const scalar epsilon = td.constProps().epsilon0();
 
-        ap =
-            (Sh/As + htc*Tc_ + epsilon*G[cellI]/4.0)
-           /(htc + epsilon*sigma*pow3(T));
-
-        bp =
-            6.0
-           *(Sh/As + htc*(Tc_ - T) + epsilon*(G[cellI]/4.0 - sigma*pow4(T)))
-           /(rho*d*cp*(ap - T));
-    }
-    else
-    {
-        ap = Tc_ + Sh/As/htc;
-        bp = 6.0*(Sh/As + htc*(Tc_ - T))/(rho*d*cp*(ap - T));
+        ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T)/htc);
+        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T)));
     }
+    bp /= rho*d*cp*(ap - T);
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
index a6297af086..dd14fe9240 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C
@@ -28,6 +28,8 @@ License
 #include "DataEntry.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
@@ -201,7 +203,7 @@ void Foam::ConeInjection<CloudType>::setProperties
 )
 {
     // set particle velocity
-    const scalar deg2Rad = constant::math::pi/180.0;
+    const scalar deg2Rad = math::pi/180.0;
 
     scalar t = time - this->SOI_;
     scalar ti = thetaInner_().value(t);
@@ -211,7 +213,7 @@ void Foam::ConeInjection<CloudType>::setProperties
     coneAngle *= deg2Rad;
     scalar alpha = sin(coneAngle);
     scalar dcorr = cos(coneAngle);
-    scalar beta = constant::math::twoPi*this->owner().rndGen().scalar01();
+    scalar beta = math::twoPi*this->owner().rndGen().scalar01();
 
     vector normal = alpha*(tanVec1_*cos(beta) + tanVec2_*sin(beta));
     vector dirVec = dcorr*direction_;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjectionMP/ConeInjectionMP.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjectionMP/ConeInjectionMP.C
index 94a58b41c4..e8eb092a39 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjectionMP/ConeInjectionMP.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjectionMP/ConeInjectionMP.C
@@ -28,6 +28,8 @@ License
 #include "DataEntry.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
@@ -247,7 +249,7 @@ void Foam::ConeInjectionMP<CloudType>::setProperties
     // set particle velocity
     const label i = parcelI%positions_.size();
 
-    const scalar deg2Rad = constant::math::pi/180.0;
+    const scalar deg2Rad = math::pi/180.0;
 
     scalar t = time - this->SOI_;
     scalar ti = thetaInner_().value(t);
@@ -257,7 +259,7 @@ void Foam::ConeInjectionMP<CloudType>::setProperties
     coneAngle *= deg2Rad;
     scalar alpha = sin(coneAngle);
     scalar dcorr = cos(coneAngle);
-    scalar beta = constant::math::twoPi*this->owner().rndGen().scalar01();
+    scalar beta = math::twoPi*this->owner().rndGen().scalar01();
 
     vector normal = alpha*(tanVec1_[i]*cos(beta) + tanVec2_[i]*sin(beta));
     vector dirVec = dcorr*axes_[i];
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
index d10171de3d..7e8fc31e18 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
@@ -28,6 +28,8 @@ License
 #include "volFields.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
@@ -128,9 +130,7 @@ Foam::FieldActivatedInjection<CloudType>::FieldActivatedInjection
 
     // Determine total volume of particles to inject
     this->volumeTotal_ =
-         nParcelsPerInjector_
-        *sum(pow3(diameters_))
-        *constant::math::pi/6.0;
+        nParcelsPerInjector_*sum(pow3(diameters_))*math::pi/6.0;
 
     // Set/cache the injector cells
     forAll(positions_, i)
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index ccb56a4767..7f5f3e6c47 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -28,6 +28,8 @@ License
 #include "mathConstants.H"
 #include "meshTools.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -195,9 +197,10 @@ Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
     {
         case pbMass:
         {
-            nP = volume/volumeTotal_
-                *massTotal_/rho
-               /(parcels*constant::math::pi/6.0*pow3(diameter));
+            nP =
+                volume/volumeTotal_
+               *massTotal_/rho
+               /(parcels*math::pi/6.0*pow3(diameter));
             break;
         }
         case pbNumber:
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
index d0018efbd6..8ebf75af0b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C
@@ -27,6 +27,8 @@ License
 #include "ManualInjection.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
@@ -106,7 +108,7 @@ Foam::ManualInjection<CloudType>::ManualInjection
     }
 
     // Determine volume of particles to inject
-    this->volumeTotal_ = sum(pow3(diameters_))*constant::math::pi/6.0;
+    this->volumeTotal_ = sum(pow3(diameters_))*math::pi/6.0;
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index ae1ee89074..322ca30284 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -28,6 +28,8 @@ License
 #include "specie.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
@@ -148,7 +150,7 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     scalarField Xc = calcXc(cellI);
 
     // droplet surface area
-    scalar A = constant::math::pi*sqr(d);
+    scalar A = math::pi*sqr(d);
 
     // calculate mass transfer of each specie in liquid
     forAll(activeLiquids_, i)