diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index 5de143de2c..72b5e0869b 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@@ -103,11 +103,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p_rgh
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
index fa06a820d8..5298950425 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
@@ -65,12 +65,7 @@ reaction->setTurbulence(turbulence());
 
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
-
-Info<< "Calculating field g.h\n" << endl;
-dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-volScalarField gh("gh", (g & mesh.C()) - ghRef);
-surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
-
+#include "gh.H"
 
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
index 8fd4fda768..0f554ef78d 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
@@ -82,11 +82,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index 91d584eb85..c54cebc9f4 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -82,11 +82,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
index b22934c60f..847c67bcf2 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
@@ -83,11 +83,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p_rgh
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index f1be9b3a37..fecec7a80d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -67,11 +67,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     // Mass flux
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
index 3ef9f32477..d281df81e9 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
@@ -51,11 +51,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     // Construct compressible turbulence model
diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H
index fc66a17f8b..988aba07b9 100644
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@@ -93,11 +93,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 5413bd193e..58960256b1 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -76,13 +76,11 @@
         incompressible::turbulenceModel::New(U, phi, mixture)
     );
 
+
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
+    #include "gh.H"
 
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
 
     volScalarField p
     (
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
index 8b1be416de..21fde57907 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
@@ -81,11 +81,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
     volScalarField p
     (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index 05bb201cb8..abebd61f17 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -69,11 +69,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
index bc4f6b362b..17253bf0a9 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@@ -55,11 +55,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index ad8b08a45c..1a01871985 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -72,11 +72,7 @@
 
     #include "readGravitationalAcceleration.H"
     #include "readhRef.H"
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+    #include "gh.H"
 
 
     volScalarField p
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index 67d7c82baa..efb2310399 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -30,12 +30,7 @@
         fluid.lookup("pMin")
     );
 
-
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
-
+    #include "gh.H"
 
     rhoThermo& thermo1 = phase1.thermo();
     rhoThermo& thermo2 = phase2.thermo();
diff --git a/src/finiteVolume/cfdTools/general/include/gh.H b/src/finiteVolume/cfdTools/general/include/gh.H
new file mode 100644
index 0000000000..b538cf69a0
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/gh.H
@@ -0,0 +1,9 @@
+    Info<< "Calculating field g.h\n" << endl;
+    dimensionedScalar ghRef
+    (
+        mag(g.value()) > SMALL
+      ? g & (cmptMag(g.value())/mag(g.value()))*hRef
+      : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
+    );
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);