diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index bd4722a01b2be7cdf31b09255ea3f86245644bb7..bb8858fb01b422ca977562927e4d13c5d38bbf0b 100644
--- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -1045,7 +1045,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                             mDotL_[i] = dmdt_[i]*L[i];
 
                             // No quenching flux
-                            //qq_[i] = 0.0;
+                            qq_[i] = 0.0;
 
                             this->operator[](i) =
                             (
@@ -1104,6 +1104,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         A1[i] = 1.0;
                         qq_[i] = 0.0;
                         mDotL_[i] = 0.0;
+                        dmdt_[i] = 0.0;
 
                         // Turbulente thermal diffusivity for single phase.
                         this->operator[](i) =
@@ -1141,12 +1142,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     Info<< "  dmdt: " << gMin((dmdt_)) << " - "
                         << gMax((dmdt_)) << endl;
 
-                     Info<< "  alphatlEff: " << gMin(liquidw*(*this + alphaw))
+                    Info<< "  alphatlEff: " << gMin(liquidw*(*this + alphaw))
                         << " - " << gMax(liquidw*(*this + alphaw)) << endl;
 
                     scalar Qeff = gSum(qEff*patch().magSf());
-                    Info<< " Effective heat transfer rate to liquid:" << Qeff
-                        << endl;
+                    Info<< " Effective heat transfer rate to liquid: " << Qeff
+                        << endl << nl;
 
                     if (debug & 2)
                     {
@@ -1186,11 +1187,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                             }
                         }
 
-                        Info<< "Sub Cool faces : " << nSubCool << endl;
-                        Info<< "Transient faces : " << nTransient << endl;
-                        Info<< "Film faces : " << nFilm << endl;
-                        Info<< "Non Boiling faces : " << nNonBoiling << endl;
-                        Info<< "Total faces : " << this->size() << endl;
+                        Info<< "Faces regime :  " <<  nl << endl;
+
+                        Info<< "    sub Cool faces : " << nSubCool << endl;
+                        Info<< "    transient faces : " << nTransient << endl;
+                        Info<< "    film faces : " << nFilm << endl;
+                        Info<< "    non-Boiling faces : " << nNonBoiling << endl;
+                        Info<< "    total faces : " << this->size() << endl << nl;
 
                         const scalarField qc
                         (
@@ -1199,7 +1202,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         );
 
                         scalar Qc = gSum(qc*patch().magSf());
-                        Info<< " Convective heat transfer:" << Qc << endl;
+                        Info<< " Convective heat transfer: " << Qc << endl;
 
                         const scalarField qFilm
                         (
@@ -1219,14 +1222,21 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         Info<< " Transient boiling heat transfer:" << Qtbtot
                             << endl;
 
-                        Info<< " tDNB: " << gMin(tDNB) << " - " << gMax(tDNB)
+                        Info<< " TDNB: " << gMin(tDNB) << " - " << gMax(tDNB)
                             << endl;
 
-                        scalar QsubCool = gSum
+                        const scalarField qSubCool
                         (
-                            fLiquid*nSubCools*(qq_ + qe())*patch().magSf()
+                            fLiquid*nSubCools*
+                            (
+                                A1*alphatConv_*hew.snGrad()
+                              + qe() + qq()
+                            )
                         );
-                        Info<< " Sub Cool boiling heat transfer:" << QsubCool
+
+                        scalar QsubCool = gSum(qSubCool*patch().magSf());
+
+                        Info<< " Sub Cool boiling heat transfer: " << QsubCool
                             << endl;
 
                         Info<< "  N: " << gMin(nSubCools*N) << " - "
diff --git a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/water/fvSchemes b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/water/fvSchemes
index 18e87173202554ea0b3ae3ccb145e58920dc4a03..908a2f4570487101d3935d3c68777cb0ae17c82e 100644
--- a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/water/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/water/fvSchemes
@@ -22,58 +22,27 @@ ddtSchemes
 gradSchemes
 {
     default             Gauss linear;
-    /*
-    grad((1-alpha.gas)) leastSquares;//Gauss linear;
-    grad(alpha.gas)     leastSquares;//Gauss linear;
-    grad(U.gas)         leastSquares;//Gauss linear;
-    grad(U.liquid)      leastSquares;//   Gauss linear;
-
-    grad(h.gas)         leastSquares;
-    grad(h.liquid)      leastSquares;
-
-    grad(alpha.liquid) leastSquares;
-    grad(alpha.gas) leastSquares;
-
-    grad(rho) leastSquares;
-
-    grad(p_rgh) leastSquares;
-
-    grad(epsilon.liquid) leastSquares;
-    grad(k.liquid) leastSquares;
-    */
-
 }
 
 divSchemes
 {
-    default                         none;//Gauss upwind;
+    default                         none;
 
     "div\(phi,alpha.*\)"            Gauss vanLeer;
     "div\(phir,alpha.*\)"           Gauss vanLeer;
 
-    "div\(alphaRhoPhi.*,U.*\)"      Gauss upwind;//limitedLinearV 1;
-    "div\(phi.*,U.*\)"              Gauss upwind;//limitedLinearV 1;
+    "div\(alphaRhoPhi.*,U.*\)"      Gauss upwind;
+    "div\(phi.*,U.*\)"              Gauss upwind;
 
-    "div\(alphaRhoPhi.*,Yi\)"       Gauss upwind;//limitedLinear 1;
-    "div\(alphaRhoPhi.*,(h|e|f).*\)"  Gauss upwind;//limitedLinear 1;
-    "div\(alphaRhoPhi.*,K.*\)"      Gauss upwind;//limitedLinear 1;
-    "div\(alphaPhi.*,p\)"           Gauss upwind;//limitedLinear 1;
+    "div\(alphaRhoPhi.*,Yi\)"       Gauss upwind;
+    "div\(alphaRhoPhi.*,(h|e|f).*\)"  Gauss upwind;
+    "div\(alphaRhoPhi.*,K.*\)"      Gauss upwind;
+    "div\(alphaPhi.*,p\)"           Gauss upwind;
 
     "div\(alphaRhoPhi.*,(k|epsilon).*\)"  Gauss upwind;
     "div\(phim,(k|epsilon)m\)"      Gauss upwind;
 
     "div\(\(\(\(alpha.*\*thermo:rho.*\)\*nuEff.*\)\*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
-
-/*
-    div(phi,U)      Gauss upwind;
-    div(phi,K)      Gauss linear;
-    div(phi,h)      Gauss upwind;
-    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) Gauss upwind;
-    div(phi,R)      Gauss upwind;
-    div(R)          Gauss linear;
-    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
-*/
 }
 
 laplacianSchemes