diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index cd254643b08b07f2055462014ca67dbeffebf8c3..aec50069be2100bab93b9afb222df649a1b0df06 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -51,7 +51,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 977fb1542a5f4df9f98af6e10107c16e7dcffdd2..8c6402f8f48366795ce1288ea6cb2aa0f1ffc5d8 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -50,9 +50,11 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Calculating field g.h\n" << endl;
-volScalarField gh("gh", g & mesh.C());
-surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
 
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 403fb6fecd1bdd5eea1de8e7eea749fb9ee6cc7f..8950122e32df910d596880a9a3088877dcce1115 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -48,7 +48,6 @@ int main(int argc, char *argv[])
 
     simpleControl simple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index 37a76ba7c6311ee63e7659c834484a9e9768c04c..d7cf3eb704d2008a6c5d671d24ff7c701977a57e 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -49,9 +49,10 @@ autoPtr<compressible::RASModel> turbulence
 );
 
 
-Info<< "Calculating field g.h\n" << endl;
-volScalarField gh("gh", g & mesh.C());
-surfaceScalarField ghf("ghf", g & mesh.Cf());
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
 
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index 1811878bd738453bf1a8f4eebc23cc7396b8d39d..f3b0655405ab9973ca2d4b74696031f8350b45a9 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -4,10 +4,11 @@ PtrList<volScalarField> rhoFluid(fluidRegions.size());
 PtrList<volVectorField> UFluid(fluidRegions.size());
 PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
 PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
-PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
 PtrList<volScalarField> ghFluid(fluidRegions.size());
 PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
+PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+PtrList<volScalarField> p_rghFluid(fluidRegions.size());
 PtrList<radiation::radiationModel> radiation(fluidRegions.size());
 
 List<scalar> initialMassFluid(fluidRegions.size());
@@ -107,31 +108,65 @@ forAll(fluidRegions, i)
         )
     );
 
-    Info<< "    Adding to turbulence\n" << endl;
-    turbulence.set
+    Info<< "    Adding to hRefFluid\n" << endl;
+    hRefFluid.set
     (
         i,
-        compressible::turbulenceModel::New
+        new uniformDimensionedScalarField
         (
-            rhoFluid[i],
-            UFluid[i],
-            phiFluid[i],
-            thermoFluid[i]
-        ).ptr()
+            IOobject
+            (
+                "hRef",
+                runTime.constant(),
+                fluidRegions[i],
+                IOobject::READ_IF_PRESENT,
+                IOobject::NO_WRITE
+            ),
+            dimensionedScalar("hRef", dimLength, 0)
+        )
+    );
+
+    dimensionedScalar ghRef
+    (
+        mag(gFluid[i].value()) > SMALL
+      ? gFluid[i]
+          & (cmptMag(gFluid[i].value())/mag(gFluid[i].value()))*hRefFluid[i]
+      : dimensionedScalar("ghRef", gFluid[i].dimensions()*dimLength, 0)
     );
 
     Info<< "    Adding to ghFluid\n" << endl;
     ghFluid.set
     (
         i,
-        new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        new volScalarField
+        (
+            "gh",
+            (gFluid[i] & fluidRegions[i].C()) - ghRef
+        )
     );
 
     Info<< "    Adding to ghfFluid\n" << endl;
     ghfFluid.set
     (
         i,
-        new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        new surfaceScalarField
+        (
+            "ghf",
+            (gFluid[i] & fluidRegions[i].Cf()) - ghRef
+        )
+    );
+
+    Info<< "    Adding to turbulence\n" << endl;
+    turbulence.set
+    (
+        i,
+        compressible::turbulenceModel::New
+        (
+            rhoFluid[i],
+            UFluid[i],
+            phiFluid[i],
+            thermoFluid[i]
+        ).ptr()
     );
 
     p_rghFluid.set
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 2e665464832c0bf9cd2ef5cbf0b343fc6e933d46..2097833ed9dcf2e4845a66ddc498110e650a6c5b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -4,10 +4,11 @@ PtrList<volScalarField> rhoFluid(fluidRegions.size());
 PtrList<volVectorField> UFluid(fluidRegions.size());
 PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
 PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
-PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
 PtrList<volScalarField> ghFluid(fluidRegions.size());
 PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
+PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+PtrList<volScalarField> p_rghFluid(fluidRegions.size());
 PtrList<radiation::radiationModel> radiation(fluidRegions.size());
 PtrList<volScalarField> KFluid(fluidRegions.size());
 PtrList<volScalarField> dpdtFluid(fluidRegions.size());
@@ -104,31 +105,65 @@ forAll(fluidRegions, i)
         )
     );
 
-    Info<< "    Adding to turbulence\n" << endl;
-    turbulence.set
+    Info<< "    Adding to hRefFluid\n" << endl;
+    hRefFluid.set
     (
         i,
-        compressible::turbulenceModel::New
+        new uniformDimensionedScalarField
         (
-            rhoFluid[i],
-            UFluid[i],
-            phiFluid[i],
-            thermoFluid[i]
-        ).ptr()
+            IOobject
+            (
+                "hRef",
+                runTime.constant(),
+                fluidRegions[i],
+                IOobject::READ_IF_PRESENT,
+                IOobject::NO_WRITE
+            ),
+            dimensionedScalar("hRef", dimLength, 0)
+        )
+    );
+
+    dimensionedScalar ghRef
+    (
+        mag(gFluid[i].value()) > SMALL
+      ? gFluid[i]
+          & (cmptMag(gFluid[i].value())/mag(gFluid[i].value()))*hRefFluid[i]
+      : dimensionedScalar("ghRef", gFluid[i].dimensions()*dimLength, 0)
     );
 
     Info<< "    Adding to ghFluid\n" << endl;
     ghFluid.set
     (
         i,
-        new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        new volScalarField
+        (
+            "gh",
+            (gFluid[i] & fluidRegions[i].C()) - ghRef
+        )
     );
 
     Info<< "    Adding to ghfFluid\n" << endl;
     ghfFluid.set
     (
         i,
-        new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        new surfaceScalarField
+        (
+            "ghf",
+            (gFluid[i] & fluidRegions[i].Cf()) - ghRef
+        )
+    );
+
+    Info<< "    Adding to turbulence\n" << endl;
+    turbulence.set
+    (
+        i,
+        compressible::turbulenceModel::New
+        (
+            rhoFluid[i],
+            UFluid[i],
+            phiFluid[i],
+            thermoFluid[i]
+        ).ptr()
     );
 
     p_rghFluid.set