diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index ea6760482a7bcbee806b9df86972e59283f53dd4..2c89e817fe94e1f7b4b89fbc4f7f9fca76f73008 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -113,15 +113,19 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
-                #include "setRegionFluidFields.H"
+                fvMesh& mesh = fluidRegions[i];
+
                 #include "readFluidMultiRegionPIMPLEControls.H"
+                #include "setRegionFluidFields.H"
                 #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
-                #include "setRegionSolidFields.H"
+                fvMesh& mesh = solidRegions[i];
+
                 #include "readSolidMultiRegionPIMPLEControls.H"
+                #include "setRegionSolidFields.H"
                 #include "solveSolid.H"
             }
 
@@ -133,8 +137,10 @@ int main(int argc, char *argv[])
 
                 forAll(fluidRegions, i)
                 {
-                    #include "setRegionFluidFields.H"
+                    fvMesh& mesh = fluidRegions[i];
+
                     #include "readFluidMultiRegionPIMPLEControls.H"
+                    #include "setRegionFluidFields.H"
                     if (!frozenFlow)
                     {
                         Info<< "\nSolving for fluid region "
@@ -166,20 +172,24 @@ int main(int argc, char *argv[])
 
                     forAll(fluidRegions, i)
                     {
+                        fvMesh& mesh = fluidRegions[i];
+
                         Info<< "\nSolving for fluid region "
                             << fluidRegions[i].name() << endl;
-                       #include "setRegionFluidFields.H"
-                       #include "readFluidMultiRegionPIMPLEControls.H"
-                       frozenFlow = true;
-                       #include "solveFluid.H"
+                        #include "readFluidMultiRegionPIMPLEControls.H"
+                        #include "setRegionFluidFields.H"
+                        frozenFlow = true;
+                        #include "solveFluid.H"
                     }
 
                     forAll(solidRegions, i)
                     {
+                        fvMesh& mesh = solidRegions[i];
+
                         Info<< "\nSolving for solid region "
                             << solidRegions[i].name() << endl;
-                        #include "setRegionSolidFields.H"
                         #include "readSolidMultiRegionPIMPLEControls.H"
+                        #include "setRegionSolidFields.H"
                         #include "solveSolid.H"
                     }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
index 21a452df9f839cf5cec07ff7dfd29cb05843dc5a..bc802e0265d244b05ba693d320b2fa78323fc542 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@@ -76,17 +76,21 @@ int main(int argc, char *argv[])
 
         forAll(fluidRegions, i)
         {
+            fvMesh& mesh = fluidRegions[i];
+
             Info<< "\nSolving for fluid region "
                 << fluidRegions[i].name() << endl;
-            #include "setRegionFluidFields.H"
             #include "readFluidMultiRegionSIMPLEControls.H"
+            #include "setRegionFluidFields.H"
             #include "solveFluid.H"
         }
 
         forAll(solidRegions, i)
         {
-            #include "setRegionSolidFields.H"
+            fvMesh& mesh = solidRegions[i];
+
             #include "readSolidMultiRegionSIMPLEControls.H"
+            #include "setRegionSolidFields.H"
             #include "solveSolid.H"
         }
 
@@ -99,8 +103,10 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
-                #include "setRegionFluidFields.H"
+                fvMesh& mesh = fluidRegions[i];
+
                 #include "readSolidMultiRegionSIMPLEControls.H"
+                #include "setRegionFluidFields.H"
                 if (!frozenFlow)
                 {
                     #include "pEqn.H"
@@ -121,20 +127,24 @@ int main(int argc, char *argv[])
 
                 forAll(fluidRegions, i)
                 {
+                    fvMesh& mesh = fluidRegions[i];
+
                     Info<< "\nSolving for fluid region "
                         << fluidRegions[i].name() << endl;
-                   #include "setRegionFluidFields.H"
-                   #include "readFluidMultiRegionSIMPLEControls.H"
-                   frozenFlow = true;
-                   #include "solveFluid.H"
+                    #include "readFluidMultiRegionSIMPLEControls.H"
+                    #include "setRegionFluidFields.H"
+                    frozenFlow = true;
+                    #include "solveFluid.H"
                 }
 
                 forAll(solidRegions, i)
                 {
+                    fvMesh& mesh = solidRegions[i];
+
                     Info<< "\nSolving for solid region "
                         << solidRegions[i].name() << endl;
-                    #include "setRegionSolidFields.H"
                     #include "readSolidMultiRegionSIMPLEControls.H"
+                    #include "setRegionSolidFields.H"
                     #include "solveSolid.H"
                 }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
index 332001c2945549089beeb3dd0a2d32f0914fa4bc..9455864173810658015e99f2519c3b5326a530cd 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
@@ -5,3 +5,5 @@
 
     const bool momentumPredictor =
         simple.getOrDefault("momentumPredictor", true);
+
+    simple.readIfPresent("frozenFlow", frozenFlowFluid[i]);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 28c4418da82aeb7aaa4185901d6fca48ae526ce2..783ebfd1da3d0fc03320eaba527ea8b5e84fc346 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -1,5 +1,3 @@
-    const fvMesh& mesh = fluidRegions[i];
-
     rhoThermo& thermo = thermoFluid[i];
     thermo.validate(args.executable(), "h", "e");
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C
index 0da321d188621f6dfd672c571ee125f67046a3fe..abc0f5f6ac240e330fa910da9ee5501742c7012e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -108,19 +108,23 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
+                fvMesh& mesh = fluidRegions[i];
+
                 Info<< "\nSolving for fluid region "
                     << fluidRegions[i].name() << endl;
-                #include "setRegionFluidFields.H"
                 #include "readFluidMultiRegionPIMPLEControls.H"
+                #include "setRegionFluidFields.H"
                 #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
+                fvMesh& mesh = solidRegions[i];
+
                 Info<< "\nSolving for solid region "
                     << solidRegions[i].name() << endl;
-                #include "setRegionSolidFields.H"
                 #include "readSolidMultiRegionPIMPLEControls.H"
+                #include "setRegionSolidFields.H"
                 #include "solveSolid.H"
             }
 
@@ -135,20 +139,24 @@ int main(int argc, char *argv[])
 
                     forAll(fluidRegions, i)
                     {
+                        fvMesh& mesh = fluidRegions[i];
+
                         Info<< "\nSolving for fluid region "
                             << fluidRegions[i].name() << endl;
-                       #include "setRegionFluidFields.H"
-                       #include "readFluidMultiRegionPIMPLEControls.H"
-                       frozenFlow = true;
-                       #include "solveFluid.H"
+                        #include "readFluidMultiRegionPIMPLEControls.H"
+                        #include "setRegionFluidFields.H"
+                        frozenFlow = true;
+                        #include "solveFluid.H"
                     }
 
                     forAll(solidRegions, i)
                     {
+                        fvMesh& mesh = solidRegions[i];
+
                         Info<< "\nSolving for solid region "
                             << solidRegions[i].name() << endl;
-                        #include "setRegionSolidFields.H"
                         #include "readSolidMultiRegionPIMPLEControls.H"
+                        #include "setRegionSolidFields.H"
                         #include "solveSolid.H"
                     }
                 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
index 7238b5a67a9477e0d98362396bb6d3f70c7863ad..9e76ce95743fe5a4c54e1f7398b59cc089b64804 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@@ -9,3 +9,5 @@
     (
         pimpleDict.getOrDefault<int>("nEnergyCorrectors", 1)
     );
+
+    pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H
index 754329349201a36ef3e8dc1f5222b4172bf5603a..c11535179be7ea8eab7dd5a3dbdb3954907be6c8 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H
@@ -1,5 +1,3 @@
-    fvMesh& mesh = fluidRegions[i];
-
     twoPhaseSystem& fluid  = phaseSystemFluid[i];
 
     phaseModel& phase1 = fluid.phase1();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
index af3cbe03c59ccf1637349a3a338e27d768ed1038..d477e4992fef4d442e4ffdf487ef23cac09649dc 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@@ -8,3 +8,5 @@
 
     const bool momentumPredictor =
         pimple.getOrDefault("momentumPredictor", true);
+
+    pimple.readIfPresent("frozenFlow", frozenFlowFluid[i]);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index a35a63d52b9e5e294b982a71119b1328b782c1b9..7647af9bd7d6bf4e7710679ed62310349b9baa3d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -1,5 +1,3 @@
-    fvMesh& mesh = fluidRegions[i];
-
     CombustionModel<rhoReactionThermo>& reaction = reactionFluid[i];
 
     rhoReactionThermo& thermo = reaction.thermo();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
index 7790ac9659ae6fe83cf9e9ab96edbd1656590bed..550325708512310e0fb6bfd88acaf06b64a0d8ef 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@@ -1,4 +1,3 @@
-fvMesh& mesh = solidRegions[i];
 solidThermo& thermo = thermos[i];
 
 tmp<volScalarField> trho = thermo.rho();