diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index 8159ee73aa67906e0c54937a76f12cb3cea3562e..5de143de2c2d71b3020ebefc6f0f309927b9da2d 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@@ -101,10 +101,14 @@
     volScalarField K("K", 0.5*magSqr(U));
 
 
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
+    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
 
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     volScalarField p_rgh
     (
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index d055be0de5d2285817150b849063329bdbb8a54a..d7ec701ad88b49863e31b7e01075695ca2cbb16c 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
index 4b901672bd0a63bbdedbef68cb190817376d4e62..fa06a820d89c5cddec9fa6fa9618e6848e44315e 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
@@ -63,9 +63,14 @@ autoPtr<compressible::turbulenceModel> turbulence
 reaction->setTurbulence(turbulence());
 
 
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+
 Info<< "Calculating field g.h\n" << endl;
-volScalarField gh("gh", g & mesh.C());
-surfaceScalarField ghf("ghf", g & mesh.Cf());
+dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+volScalarField gh("gh", (g & mesh.C()) - ghRef);
+surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+
 
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
index 4d6ca92e56c8cc8d38b6c61f4839dfceeba9fd5c..d8e68441c9d18e8ca6d29923d5843935355342c6 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
@@ -48,7 +48,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index c5be6738678f1668f9cd0bd86a76631ce43726da..f107467510800030b0040d7621d882e5de22d015 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createIncompressibleRadiationModel.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
index a7243f28179f9858f8dae1d4cfefe3e5013be39c..8fd4fda768a384337d2bccfd4a2a4f5745fdc69e 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
@@ -79,9 +79,15 @@
         mesh
     );
 
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index 42a843365d179b7ff240a32f613afdf4c52eb446..fbf754f6de591f1cfb2566e4ec46f53f56ee1a66 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -62,7 +62,6 @@ int main(int argc, char *argv[])
 
     simpleControl simple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index 242350fe64b7479551b0ff7874b02ce81751819d..91d584eb85c3b6a836dd25feb05f27ac760e21a5 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -79,9 +79,15 @@
         mesh
     );
 
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
index d3fffd024401a3d643095b4e5f149b88ad246793..b22934c60f4681d2e04bf15b56176b7e1560d57c 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
@@ -80,10 +80,15 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
 
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
+    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
 
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     volScalarField p_rgh
     (
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 9c2260696bca9371058fb0a148a4443edbed7cb5..757ede7fee48bad6196c3ae39a1d49c5990b76e7 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -52,7 +52,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 7f225f4607a1dd10269e96ef9aaa248008a3c705..9d56a2334f9c6dc8522c0b34690e1e89b55ac99e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -56,7 +56,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
@@ -98,8 +97,8 @@ int main(int argc, char *argv[])
                     << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                     << " s" << endl;
 
-                gh = g & mesh.C();
-                ghf = g & mesh.Cf();
+                gh = (g & mesh.C()) - ghRef;
+                ghf = (g & mesh.Cf()) - ghRef;
             }
 
             if (mesh.changing() && correctPhi)
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 0559985b54617a9e4990aacb03e1df997d84eead..de2c07e1488d9ad314d5d1add80f12823e712ea7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
 
     pimpleControl pimple(mesh);
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 5f7b4dddadf88dc0621b0c3e24567fad7327d1d6..f1be9b3a3719104112c4eb9009f8c90e21fa0989 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -57,11 +57,22 @@
     );
 
 
-    dimensionedScalar pMin(mixture.lookup("pMin"));
+    dimensionedScalar pMin
+    (
+        "pMin",
+        dimPressure,
+        mixture.lookup("pMin")
+    );
+
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
 
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+
 
     // Mass flux
     // Initialisation does not matter because rhoPhi is reset after the
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
index 24ba56e3e5ef257b86c00c884800eb91dc68e346..21ef8128899fb6c42ccda833cec442bdf34d6dc4 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
@@ -48,7 +48,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
 
     pimpleControl pimple(mesh);
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
index 7c864cc569a3a5ead4267ea8f1f823484d32b46f..3ef9f32477d7203130f32d7abefa904ad10a75ed 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
@@ -49,9 +49,14 @@
     dimensionedScalar pMin(mixture.lookup("pMin"));
 
 
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+
 
     // Construct compressible turbulence model
     autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H
index bfd899dd98649b66adade837c95cfc2fd951ff5e..fc66a17f8b3584fb22e3741aa265c3b8d7a5d550 100644
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@@ -90,9 +90,15 @@
         ::New(rho, U, rhoPhi, mixture)
     );
 
-    Info<< "Calculating field (g.h)f\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    #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);
+
 
     volScalarField p
     (
diff --git a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
index fc366873fa356facb3472466c543ddda620eabc8..35d4b441d0e9aeb08347977464a4f5e26f3264c2 100644
--- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
+++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
@@ -57,7 +57,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 93465ed299e0beea84f42b3f8693b2333d3624fb..5413bd193e9512729b02eb164136017a45e9c38e 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -77,10 +77,12 @@
     );
 
     #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
 
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 4df1c23dfc3b722bd55ec6abc77adc7ea8cd5ac0..84d1283d04a60a188ade032d27f8a5807eee811e 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -110,8 +110,8 @@ int main(int argc, char *argv[])
                         << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                         << " s" << endl;
 
-                    gh = g & mesh.C();
-                    ghf = g & mesh.Cf();
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
                 }
 
                 if (mesh.changing() && correctPhi)
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
index 284d9ed035e1df879482de5b8b5e71e3d3f38c0f..8b1be416de750cc97fdea3537296fee4b81a0466 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
@@ -79,10 +79,13 @@
         incompressible::turbulenceModel::New(U, phi, mixture)
     );
 
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
 
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index d51f5a83ab2823a4dd0ebfc4dcc6e997eda55ab0..abb18c4ba308877cbc1852b9045cbc3834dd77e1 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -50,7 +50,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index d8b06da38d1bd189a04a8b2abe06d48651e6a2cf..05bb201cb87bdcde700f68f826e07ca4d110c775 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -67,9 +67,14 @@
     );
 
 
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 5c05875e35811c63f0a837fe2bc1f58fca785c95..a1a416f08630c7162868f57798a1dd3a96819d1f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"
@@ -122,8 +121,8 @@ int main(int argc, char *argv[])
                         << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                         << " s" << endl;
 
-                    gh = g & mesh.C();
-                    ghf = g & mesh.Cf();
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
                 }
 
                 if (mesh.changing() && correctPhi)
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index 2364208b5a760aa1dbbe403f9048ed0b111fb064..8f5c97522fbc7d4553190c089fd844a241a91001 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -60,7 +60,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
index 26d7a17dc019b24565be7476609dfdeabe6d23c4..bc4f6b362b8230c6c90e6d0ae5a9f721415fbd2b 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@@ -52,11 +52,15 @@
         incompressible::turbulenceModel::New(U, phi, mixture)
     );
 
+
     #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
 
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.C
index f75e6d9c4bd227a2112e828160690e9d4c847097..fff99c0b77c9bb049b8edd45dc905bb7bbef8a6d 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.C
@@ -106,8 +106,8 @@ int main(int argc, char *argv[])
                         << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                         << " s" << endl;
 
-                    gh = g & mesh.C();
-                    ghf = g & mesh.Cf();
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
                 }
 
                 if (mesh.changing() && correctPhi)
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index 442bb9083b6dabe655615a262aeffe775ffd144a..ad8b08a45c6a34e926094b08995b102e9a125b54 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -69,9 +69,15 @@
         incompressible::turbulenceModel::New(U, phi, mixture)
     );
 
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    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/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 283724286b7d7b992790a119a1857bdf7371b12f..dab88b930aff6dccc07dd32d83d8002e20abcc58 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -49,7 +49,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
diff --git a/src/finiteVolume/cfdTools/general/include/readhRef.H b/src/finiteVolume/cfdTools/general/include/readhRef.H
new file mode 100644
index 0000000000000000000000000000000000000000..d34c3e7dd9d121a6042298fd6bebdb919b087f22
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/readhRef.H
@@ -0,0 +1,13 @@
+    Info<< "\nReading hRef" << endl;
+    uniformDimensionedScalarField hRef
+    (
+        IOobject
+        (
+            "hRef",
+            runTime.constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        dimensionedScalar("hRef", dimLength, 0)
+    );
diff --git a/tutorials/multiphase/LTSInterFoam/DTCHull/constant/hRef b/tutorials/multiphase/LTSInterFoam/DTCHull/constant/hRef
new file mode 100644
index 0000000000000000000000000000000000000000..bda49ba88369c5f285ac5508afde615fc9835de3
--- /dev/null
+++ b/tutorials/multiphase/LTSInterFoam/DTCHull/constant/hRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedScalarField;
+    location    "constant";
+    object      hRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+value           0.244;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/hRef b/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/hRef
new file mode 100644
index 0000000000000000000000000000000000000000..bda49ba88369c5f285ac5508afde615fc9835de3
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/hRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedScalarField;
+    location    "constant";
+    object      hRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+value           0.244;
+
+// ************************************************************************* //