diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 48aecba36956..cb0c97cf3e13 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -920,7 +920,6 @@ class CConfig {
   Initial_BCThrust,                /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
   Pressure_FreeStream,             /*!< \brief Total pressure of the fluid. */
   Pressure_Thermodynamic,          /*!< \brief Thermodynamic pressure of the fluid. */
-  Standard_Ref_Temperature,        /*!< \brief Standard reference temperature for multicomponent flows. */
   Temperature_FreeStream,          /*!< \brief Total temperature of the fluid.  */
   Temperature_ve_FreeStream;       /*!< \brief Total vibrational-electronic temperature of the fluid.  */
   unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
@@ -1960,12 +1959,6 @@ class CConfig {
    */
   su2double GetPressure_Thermodynamic(void) const { return Pressure_Thermodynamic; }
 
-  /*!
-   * \brief Get the value of the standard reference temperature for multicomponent flows.
-   * \return Standard reference temperature, Non-dimensionalized if it is needed for Non-Dimensional problems.
-   */
-  su2double GetStandard_RefTemperatureND(void) const { return Standard_Ref_Temperature / Temperature_Ref; }
-
   /*!
    * \brief Get the value of the non-dimensionalized thermodynamic pressure.
    * \return Non-dimensionalized thermodynamic pressure.
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 88bc7d4f8743..6c66780d84c2 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -96,7 +96,7 @@ const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598;   /*!< \brief Universal gas
 const su2double BOLTZMANN_CONSTANT = 1.3806503E-23;   /*!< \brief Boltzmann's constant [J K^-1] */
 const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
 const su2double FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */
-
+const su2double STD_REF_TEMP = 298.15;  /*!< \brief Standard reference temperature for enthalpy in Kelvin. */
 const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */
 
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 4d1a262d4b1f..9e30278ee10c 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1229,9 +1229,6 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("GAMMA_VALUE", Gamma, 1.4);
   /*!\brief THERMODYNAMIC_PRESSURE  \n DESCRIPTION: Thermodynamics(operating) Pressure (101325 Pa), only for incompressible flows) \ingroup Config*/
   addDoubleOption("THERMODYNAMIC_PRESSURE", Pressure_Thermodynamic, 101325.0);
-  /*!\brief STANDARD_REFERENCE_TEMPERATURE  \n DESCRIPTION: Standard reference temperature (298.15K), only for
-   * multicomponent incompressible flows) \ingroup Config*/
-  addDoubleOption("STANDARD_REFERENCE_TEMPERATURE", Standard_Ref_Temperature, 298.15);
   /*!\brief CP_VALUE  \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/
   addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp);
   /*!\brief THERMAL_EXPANSION_COEFF  \n DESCRIPTION: Thermal expansion coefficient (0.00347 K^-1 (air), used for Boussinesq approximation for liquids/non-ideal gases) \ingroup Config*/
diff --git a/SU2_CFD/include/fluid/CConstantDensity.hpp b/SU2_CFD/include/fluid/CConstantDensity.hpp
index acfeed2e0cf0..a796762a2713 100644
--- a/SU2_CFD/include/fluid/CConstantDensity.hpp
+++ b/SU2_CFD/include/fluid/CConstantDensity.hpp
@@ -56,16 +56,16 @@ class CConstantDensity final : public CFluidModel {
        decoupled equation. Hence, we update the value.
        Note Cp = Cv, (gamma = 1).*/
     Temperature = t;
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - STD_REF_TEMP);  // Sensible enthalpy relative to STD_REF_TEMP
   }
 
   /*!
    * \brief Set the Dimensionless State using Enthalpy.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - not used here. 
+   * \param[in] val_scalars - not used here.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + STD_REF_TEMP;  // Temperature from sensible enthalpy
   }
 };
diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp
index 164043872d6a..ba97fa5451e4 100644
--- a/SU2_CFD/include/fluid/CFluidScalar.hpp
+++ b/SU2_CFD/include/fluid/CFluidScalar.hpp
@@ -42,7 +42,6 @@ class CFluidScalar final : public CFluidModel {
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
   su2double Gas_Constant;                 /*!< \brief Specific gas constant. */
   const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
-  const su2double Ref_Temperature;        /*!< \brief Standard Reference temperature, usually set to 298.15 K. */
   const su2double GasConstant_Ref;        /*!< \brief Gas constant reference needed for Nondimensional problems. */
   const su2double Prandtl_Turb_Number;    /*!< \brief Prandlt turbulent number.*/
   const su2double Schmidt_Turb_Number;    /*!< \brief Schmidt turbulent number.*/
@@ -177,7 +176,7 @@ class CFluidScalar final : public CFluidModel {
   /*!
    * \brief Virtual member.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - Scalar mass fractions. 
+   * \param[in] val_scalars - Scalar mass fractions.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
 };
diff --git a/SU2_CFD/include/fluid/CIncIdealGas.hpp b/SU2_CFD/include/fluid/CIncIdealGas.hpp
index 16266e569568..4b35735f096d 100644
--- a/SU2_CFD/include/fluid/CIncIdealGas.hpp
+++ b/SU2_CFD/include/fluid/CIncIdealGas.hpp
@@ -58,7 +58,7 @@ class CIncIdealGas final : public CFluidModel {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - STD_REF_TEMP);  // Sensible enthalpy relative to REF_TEMP
   }
 
   /*!
@@ -67,7 +67,7 @@ class CIncIdealGas final : public CFluidModel {
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + STD_REF_TEMP;  // Temperature from sensible enthalpy
     Density = Pressure / (Temperature * Gas_Constant);
   }
 
diff --git a/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp b/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
index e1f5f1d9105c..aff8085617d9 100644
--- a/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
+++ b/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
@@ -58,11 +58,10 @@ class CIncIdealGasPolynomial final : public CFluidModel {
    * \param[in] config - configuration container for the problem.
    */
   void SetCpModel(const CConfig* config) override {
-    const su2double t_ref = config->GetStandard_RefTemperatureND();
     Enthalpy_Ref = 0.0;
     su2double t_i = 1.0;
     for (int i = 0; i < N; ++i) {
-      t_i *= t_ref;
+      t_i *= STD_REF_TEMP;
       coeffs_[i] = config->GetCp_PolyCoeffND(i);
       Enthalpy_Ref += coeffs_[i] * t_i / (i + 1);
     }
diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp
index 5573949ea62b..dd9821f8610d 100644
--- a/SU2_CFD/include/variables/CIncEulerVariable.hpp
+++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp
@@ -72,7 +72,6 @@ class CIncEulerVariable : public CFlowVariable {
   VectorType Streamwise_Periodic_RecoveredPressure,    /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */
              Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */
   su2double TemperatureLimits[2];                      /*!< \brief Temperature limits [K]. */
-  su2double TemperatureInc = 0.0;                      /*!< \brief Temperature [K] imposed when energy equation is switch off. */
  public:
   /*!
    * \brief Constructor of the class.
diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp
index b42cf7b70433..4b384c521f9c 100644
--- a/SU2_CFD/src/fluid/CFluidScalar.cpp
+++ b/SU2_CFD/src/fluid/CFluidScalar.cpp
@@ -45,7 +45,6 @@ CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* co
     : CFluidModel(),
       n_species_mixture(config->GetnSpecies() + 1),
       Pressure_Thermodynamic(value_pressure_operating),
-      Ref_Temperature(config->GetStandard_RefTemperatureND()),
       GasConstant_Ref(config->GetGas_Constant_Ref()),
       Prandtl_Turb_Number(config->GetPrandtl_Turb()),
       Schmidt_Turb_Number(config->GetSchmidt_Number_Turbulent()),
@@ -219,14 +218,14 @@ su2double CFluidScalar::ComputeEnthalpyFromT(const su2double val_temperature, co
    * depend on temperature, but it does depend on mixture composition, enthalpy is directly computed from
    * the expression h_s = Cp(T - T_ref).
    */
-  su2double val_Enthalpy = Cp * (val_temperature - Ref_Temperature);
+  su2double val_Enthalpy = Cp * (val_temperature - STD_REF_TEMP);
   return val_Enthalpy;
 }
 
 void CFluidScalar::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const {
-  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - Ref_Temperature);
+  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - STD_REF_TEMP);
   for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) {
-    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - Ref_Temperature);
+    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - STD_REF_TEMP);
     enthalpy_diffusions[iVar] = Density * (enthalpy_species_i * massDiffusivity[iVar] -
                                            enthalpy_species_N * massDiffusivity[n_species_mixture - 1]);
     enthalpy_diffusions[iVar] += Mu_Turb * (enthalpy_species_i - enthalpy_species_N) / Schmidt_Turb_Number;
@@ -271,7 +270,7 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
    * depend on temperature, but it does depend on mixture composition, temperature is directly solved from the
    * expression h_s = Cp(T - T_ref).
    */
-  Temperature = val_enthalpy / Cp + Ref_Temperature;
+  Temperature = val_enthalpy / Cp + STD_REF_TEMP;
   Density = Pressure_Thermodynamic / (Temperature * Gas_Constant);
   Cv = Cp - Gas_Constant;
 
@@ -282,4 +281,5 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
   }
 
   Kt = WilkeConductivity(val_scalars);
+  ComputeMassDiffusivity();
 }
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
index 9d408892cf6e..f27c77cf5569 100644
--- a/SU2_CFD/src/output/CFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -330,7 +330,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("LAMINAR_VISCOSITY", "Laminar_Viscosity", "PRIMITIVE", "Laminar viscosity");
     AddVolumeOutput("HEAT_CAPACITY", "Heat_Capacity", "PRIMITIVE", "Heat capacity");
     AddVolumeOutput("THERMAL_CONDUCTIVITY", "Thermal_Conductivity", "PRIMITIVE", "Thermal conductivity");
-    if (heat || flamelet) AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
+    if (!weakly_coupled_heat)
+      AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
 
     AddVolumeOutput("SKIN_FRICTION-X", "Skin_Friction_Coefficient_x", "PRIMITIVE", "x-component of the skin friction vector");
     AddVolumeOutput("SKIN_FRICTION-Y", "Skin_Friction_Coefficient_y", "PRIMITIVE", "y-component of the skin friction vector");
@@ -354,7 +355,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component");
   if (nDim == 3)
     AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component");
-  if (config->GetEnergy_Equation()){
+  if (heat){
     AddVolumeOutput("RES_ENTHALPY", "Residual_Enthalpy", "RESIDUAL", "Residual of the enthalpy");
   }
   SetVolumeOutputFieldsScalarResidual(config);
@@ -440,7 +441,8 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
     SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint));
     SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint));
     SetVolumeOutputValue("THERMAL_CONDUCTIVITY", iPoint, Node_Flow->GetThermalConductivity(iPoint));
-    if (heat || flamelet) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
+    if (!weakly_coupled_heat)
+      SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
   }
 
   SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0));
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 956af94ab58d..b9d8c013bb65 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -2145,9 +2145,9 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo
      Therefore, we build inv(Precon) here and multiply by the residual
      later in the R-K and Euler Explicit time integration schemes. ---*/
 
-    
+
     Preconditioner[0][0] = Enthalpy * BetaInc2 * dRhodh / Density + BetaInc2;
-    
+
     for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][0] = -1.0 * Velocity[iDim] / Density;
 
     if (energy) {
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 28747c7e1d71..383c55945ec8 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -892,14 +892,15 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
       SU2_MPI::Allreduce(&globalCounter1, &notConvergedCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       SU2_MPI::Allreduce(&globalCounter2, &smallYPlusCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
 
-      if (rank == MASTER_NODE) {
+      const bool write_iter = (config->GetTimeIter() % config->GetVolumeOutputFrequency(0) == 0);
+      if (rank == MASTER_NODE && write_iter) {
         if (notConvergedCounter)
           cout << "Warning: Computation of wall coefficients (y+) did not converge in "
                << notConvergedCounter << " points." << endl;
 
         if (smallYPlusCounter)
-          cout << "Warning: y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
-               << " points, for which the wall model is not active." << endl;
+          cout << "y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
+               << " points. No problem, but you can increase your near-wall mesh size." << endl;
       }
     }
     END_SU2_OMP_SAFE_GLOBAL_ACCESS
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
index ab90660ef22e..be79d762c7c8 100644
--- a/SU2_CFD/src/solvers/CNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -1028,14 +1028,16 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
       SU2_MPI::Allreduce(&globalCounter1, &notConvergedCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       SU2_MPI::Allreduce(&globalCounter2, &smallYPlusCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
 
-      if (rank == MASTER_NODE) {
+      const bool write_iter = (config->GetInnerIter() % config->GetVolumeOutputFrequency(0) == 0);
+
+      if (rank == MASTER_NODE && write_iter) {
         if (notConvergedCounter)
           cout << "Warning: Computation of wall coefficients (y+) did not converge in "
                << notConvergedCounter << " points." << endl;
 
         if (smallYPlusCounter)
-          cout << "Warning: y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
-               << " points, for which the wall model is not active." << endl;
+          cout << "y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
+               << " points. No problem, but you can increase your near-wallmesh size." << endl;
       }
     }
     END_SU2_OMP_SAFE_GLOBAL_ACCESS
diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp
index 3f25af7d1183..985b67426c46 100644
--- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp
+++ b/SU2_CFD/src/solvers/CSpeciesSolver.cpp
@@ -207,12 +207,11 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi
 
     const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE);
     const bool energy = config->GetEnergy_Equation();
-    const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET);
     const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat();
 
     /*--- for the flamelet model, the temperature is saved to file, but the energy equation is off ---*/
 
-   if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--;
+   if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--;
 
     /*--- Load data from the restart into correct containers. ---*/
 
@@ -296,22 +295,23 @@ void CSpeciesSolver::Preprocessing(CGeometry* geometry, CSolver** solver_contain
                                    bool Output) {
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
-  /*--- Set the laminar mass Diffusivity for the species solver. ---*/
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
     const su2double temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
     const su2double* scalar = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint);
     solver_container[FLOW_SOL]->GetFluidModel()->SetMassDiffusivityModel(config);
     solver_container[FLOW_SOL]->GetFluidModel()->SetTDState_T(temperature, scalar);
+    /*--- Recompute viscosity, important  to get diffusivity correct across MPI ranks. ---*/
+    nodes->SetLaminarViscosity(iPoint, solver_container[FLOW_SOL]->GetFluidModel()->GetLaminarViscosity());
+    /*--- Set the laminar mass Diffusivity for the species solver. ---*/
     for (auto iVar = 0u; iVar <= nVar; iVar++) {
       const su2double mass_diffusivity = solver_container[FLOW_SOL]->GetFluidModel()->GetMassDiffusivity(iVar);
       nodes->SetDiffusivity(iPoint, mass_diffusivity, iVar);
     }
-
   }  // iPoint
   END_SU2_OMP_FOR
 
-  /*--- Clear Residual and Jacobian. Upwind second order reconstruction and gradients ---*/
+  /*--- Clear Residual and Jacobian. Upwind second order reconstruction and gradients. ---*/
   CommonPreprocessing(geometry, config, Output);
 }
 
diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp
index 7d7e6fd26730..05f456883fb2 100644
--- a/SU2_CFD/src/solvers/CTurbSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSolver.cpp
@@ -131,9 +131,8 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig*
     const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE);
     const bool energy = config->GetEnergy_Equation();
     const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat();
-    const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET);
 
-    if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--;
+    if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--;
 
     /*--- Load data from the restart into correct containers. ---*/
 
@@ -252,7 +251,7 @@ void CTurbSolver::ComputeUnderRelaxationFactorHelper(su2double allowableRatio) {
 
   /* Loop over the solution update given by relaxing the linear
    system for this nonlinear iteration. */
-  
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
     su2double localUnderRelaxation = 1.0;
@@ -265,16 +264,16 @@ void CTurbSolver::ComputeUnderRelaxationFactorHelper(su2double allowableRatio) {
       turbulence variables can change over a nonlinear iteration. */
       if (ratio > allowableRatio) {
         localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
-      
+
       }
     }
 
     /* Threshold the relaxation factor in the event that there is
      a very small value. This helps avoid catastrophic crashes due
      to non-realizable states by canceling the update. */
-    
+
     if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0;
-    
+
     /* Store the under-relaxation factor for this point. */
 
     nodes->SetUnderRelaxation(iPoint, localUnderRelaxation);
diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp
index 881de86562e4..0c117127c8b6 100644
--- a/SU2_CFD/src/variables/CIncEulerVariable.cpp
+++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp
@@ -75,11 +75,10 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
   SetPressure(iPoint);
 
-  /*--- Set the value of the enthalpy directly ---*/
-
   su2double Enthalpy = Solution(iPoint, nDim +1);
   FluidModel->SetTDState_h(Enthalpy);
   su2double Temperature = FluidModel->GetTemperature();
+
   const auto check_temp = SetTemperature(iPoint, Temperature, TemperatureLimits);
 
   /*--- Use the fluid model to compute the new value of density.
@@ -103,8 +102,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
     Enthalpy = Solution(iPoint, nDim+1);
     FluidModel->SetTDState_h(Enthalpy);
-    Temperature = FluidModel->GetTemperature();
-    SetTemperature(iPoint, Temperature, TemperatureLimits);
+    SetTemperature(iPoint, FluidModel->GetTemperature(), TemperatureLimits);
     SetDensity(iPoint, FluidModel->GetDensity());
 
     /*--- Flag this point as non-physical. ---*/
@@ -117,10 +115,13 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
   SetVelocity(iPoint);
 
-  /*--- Set enthalpy and specific heats (only necessary for consistency with preconditioning). ---*/
+  /*--- Set specific heats ---*/
 
   SetSpecificHeatCp(iPoint, FluidModel->GetCp());
   SetSpecificHeatCv(iPoint, FluidModel->GetCv());
+
+  /*--- Set enthalpy ---*/
+
   SetEnthalpy(iPoint, FluidModel->GetEnthalpy());
 
   return physical;
diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp
index e31109fcd662..6a335e103df4 100644
--- a/SU2_CFD/src/variables/CIncNSVariable.cpp
+++ b/SU2_CFD/src/variables/CIncNSVariable.cpp
@@ -47,34 +47,26 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su
     AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);
     Grad_AuxVar.resize(nPoint,nAuxVar,nDim);
   }
-  if(!Energy) TemperatureInc = config->GetInc_Temperature_Init();
 }
 
 bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel, const su2double *scalar) {
 
   bool physical = true;
-  su2double Temperature;
-  su2double Enthalpy;
 
   /*--- Set the value of the pressure ---*/
 
   SetPressure(iPoint);
 
-  if (Energy) {
-    Enthalpy = Solution(iPoint, nDim + 1);
-    FluidModel->SetTDState_h(Enthalpy, scalar);
-    Temperature = FluidModel->GetTemperature();
-  } else {
-    /*--- When energy equation is switch off, a constant temperature is imposed, and enthalpy is recomputed based on
-     * this temperature. As in the fluid flamelet model, the temperature is retrieved from a look-up table, then the
-     * temperature is obtained directly from the fluidmodel. For the other fluid models, GetTemperature provides the
-     * same value as TemperatureInc ---*/
-    FluidModel->SetTDState_T(TemperatureInc, scalar);
-    Enthalpy = Solution(iPoint, nDim + 1) = FluidModel->GetEnthalpy();
-    Temperature = FluidModel->GetTemperature();
-  }
+  su2double Enthalpy = Solution(iPoint, nDim + 1);
+  FluidModel->SetTDState_h(Enthalpy, scalar);
+  su2double Temperature = FluidModel->GetTemperature();
+
   auto check_temp = SetTemperature(iPoint, Temperature, TemperatureLimits);
 
+  /*--- Use the fluid model to compute the new value of density.
+  Note that the thermodynamic pressure is constant and decoupled
+  from the dynamic pressure being iterated. ---*/
+
   /*--- Set the value of the density ---*/
 
   const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity());
@@ -93,7 +85,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do
     Enthalpy = Solution(iPoint, nDim + 1);
     FluidModel->SetTDState_h(Enthalpy, scalar);
     SetTemperature(iPoint, FluidModel->GetTemperature(), TemperatureLimits);
-
     SetDensity(iPoint, FluidModel->GetDensity());
 
     /*--- Flag this point as non-physical. ---*/
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index fc6ee002821f..5b1642e92f35 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -375,7 +375,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-4.432484, 0.507522, 0.000000, 0.000000]
+    inc_buoyancy.test_vals = [-3.860456, -3.122396, 4.074516, -23.415000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -665,7 +665,7 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals         = [19.000000, -1.048446, -1.324274]
+    slinc_steady.test_vals         = [19.000000, -1.154873, -1.378119]
     slinc_steady.test_vals_aarch64 = [19.000000, -1.048446, -1.324274]
     slinc_steady.multizone = True
     test_list.append(slinc_steady)
diff --git a/TestCases/incomp_navierstokes/buoyancy_cavity/lam_buoyancy_cavity.cfg b/TestCases/incomp_navierstokes/buoyancy_cavity/lam_buoyancy_cavity.cfg
index 47474c4a0a74..f7d53887c85f 100644
--- a/TestCases/incomp_navierstokes/buoyancy_cavity/lam_buoyancy_cavity.cfg
+++ b/TestCases/incomp_navierstokes/buoyancy_cavity/lam_buoyancy_cavity.cfg
@@ -62,14 +62,14 @@ REF_AREA= 1.0
 MARKER_HEATFLUX= ( upper, 0.0, lower, 0.0 )
 MARKER_ISOTHERMAL= ( left, 461.04, right, 115.26 )
 MARKER_PLOTTING= ( upper, left, right, lower )
-MARKER_MONITORING= ( NONE )
+MARKER_MONITORING= ( left, right,upper )
 
 % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
 %
 NUM_METHOD_GRAD= GREEN_GAUSS
-CFL_NUMBER= 1e2
-CFL_ADAPT= NO
-CFL_ADAPT_PARAM= ( 1.5, 0.5, 15.0, 1e10)
+CFL_NUMBER= 20
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.95, 1.01, 10, 1000, 0.001, 0)
 MAX_DELTA_TIME= 1E6
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
 ITER= 99999
@@ -79,8 +79,8 @@ ITER= 99999
 LINEAR_SOLVER= FGMRES
 LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ILU_FILL_IN= 0
-LINEAR_SOLVER_ERROR= 1E-15
-LINEAR_SOLVER_ITER= 20
+LINEAR_SOLVER_ERROR= 1E-06
+LINEAR_SOLVER_ITER= 25
 
 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
 %
@@ -114,7 +114,6 @@ GRAD_OBJFUNC_FILENAME= of_grad
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
 OUTPUT_WRT_FREQ= 100
-SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_ENTHALPY, LIFT, DRAG)
-
+SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_ENTHALPY, TOTAL_HEATFLUX)
 % ----------------------- GEOMETRY EVALUATION PARAMETERS ----------------------%
 GEO_BOUNDS= ( 0.499, 0.501)
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 9227d447d597..643649b7af55 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -603,7 +603,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-4.435827, 0.508037, 0.000000, 0.000000]
+    inc_buoyancy.test_vals = [-3.860462, -3.122401, 4.074510, -23.414000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -635,7 +635,7 @@ def main():
     inc_heatTransfer_BC.cfg_dir   = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d"
     inc_heatTransfer_BC.cfg_file  = "BC_HeatTransfer.cfg"
     inc_heatTransfer_BC.test_iter = 50
-    inc_heatTransfer_BC.test_vals = [-8.747445, -7.653823, -8.079512, -0.654920, -1670.100000]
+    inc_heatTransfer_BC.test_vals = [-8.686723, -7.601245, -8.072766, 1.864213, -1669.800000]
     test_list.append(inc_heatTransfer_BC)
 
     ############################
@@ -1209,7 +1209,7 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals = [19.000000, -1.012869, -1.222087]
+    slinc_steady.test_vals = [19.000000, -1.131557, -1.370471]
     slinc_steady.timeout   = 100
     slinc_steady.tol       = 0.00002
     slinc_steady.multizone = True
@@ -1318,7 +1318,7 @@ def main():
     p1rad           = TestCase('p1rad')
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
-    p1rad.test_iter = 100
+    p1rad.test_iter = 50
     p1rad.test_vals = [-7.743666, -7.921411, -2.111848, 0.098302, -47.897000]
     test_list.append(p1rad)
 
@@ -1363,7 +1363,7 @@ def main():
     sp_pinArray_cht_2d_dp_hf.cfg_dir   = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d"
     sp_pinArray_cht_2d_dp_hf.cfg_file  = "configMaster.cfg"
     sp_pinArray_cht_2d_dp_hf.test_iter = 100
-    sp_pinArray_cht_2d_dp_hf.test_vals = [0.379544, 2.491817, -1.255558, -0.613405, 208.023676, 349.870000, -0.000000, -0.613410, 0.613410]
+    sp_pinArray_cht_2d_dp_hf.test_vals = [0.353740, 2.436726, -1.270003, -0.615079, 208.023676, 350.100000, -0.000000, -0.615080, 0.615080]
     sp_pinArray_cht_2d_dp_hf.multizone = True
     test_list.append(sp_pinArray_cht_2d_dp_hf)
 
@@ -1372,7 +1372,7 @@ def main():
     sp_pinArray_3d_cht_mf_hf_tp.cfg_dir   = "incomp_navierstokes/streamwise_periodic/chtPinArray_3d"
     sp_pinArray_3d_cht_mf_hf_tp.cfg_file  = "configMaster.cfg"
     sp_pinArray_3d_cht_mf_hf_tp.test_iter = 30
-    sp_pinArray_3d_cht_mf_hf_tp.test_vals = [0.261726, 3.878700, 0.298987, -0.009537, 104.749181, 354.800000, 0.000000]
+    sp_pinArray_3d_cht_mf_hf_tp.test_vals = [-1.279608, 3.011923, -0.641797, -0.009777, 104.609505, 410.330000, 0.000000]
     sp_pinArray_3d_cht_mf_hf_tp.multizone = True
     test_list.append(sp_pinArray_3d_cht_mf_hf_tp)
 
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg b/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg
index 537eee85ce8b..c01c2de84cb1 100644
--- a/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg
+++ b/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg
@@ -25,13 +25,12 @@ FREESTREAM_TEMPERATURE = 673
 FREESTREAM_DENSITY = 2.55
 % ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
 %
-%INC_DENSITY_MODEL= CONSTANT
 INC_DENSITY_MODEL= VARIABLE
 INC_DENSITY_INIT= 2.55
 %
 INC_VELOCITY_INIT= (40.00, 0.0, 0.0 )
 %
-INC_ENERGY_EQUATION= YES
+INC_ENERGY_EQUATION= NO
 INC_TEMPERATURE_INIT= 673.0
 %
 INC_NONDIM= DIMENSIONAL
@@ -41,20 +40,26 @@ INC_NONDIM= DIMENSIONAL
 %FLUID_MODEL= CONSTANT_DENSITY
 FLUID_MODEL= INC_IDEAL_GAS
 %
-CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY
-THERMAL_CONDUCTIVITY_CONSTANT= 0.0357
+CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
+%THERMAL_CONDUCTIVITY_CONSTANT= 0.0357
 %
 PRANDTL_LAM= 0.72
-TURBULENT_CONDUCTIVITY_MODEL= NONE
+TURBULENT_CONDUCTIVITY_MODEL= CONSTANT_PRANDTL_TURB
+%TURBULENT_CONDUCTIVITY_MODEL= NONE 
 PRANDTL_TURB= 0.90
 %
-VISCOSITY_MODEL= SUTHERLAND
+%VISCOSITY_MODEL= SUTHERLAND
+VISCOSITY_MODEL= POLYNOMIAL_VISCOSITY 
+% polynomial viscosity for methane-air (phi=0.50) 
+MU_POLYCOEFFS= (1.097e-5, 3.437e-8, 3.332e-12, 0, 0)
+
+
 MU_CONSTANT=  1.716E-5
 MU_REF = 1.716e-5
 MU_T_REF= 273.15
 SUTHERLAND_CONSTANT = 110.4
 
-SPECIFIC_HEAT_CP = 1150
+SPECIFIC_HEAT_CP = 1350
 %
 % -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
 %
@@ -108,9 +113,12 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT
 % -------------------- SCALAR TRANSPORT ---------------------------------------%
 %
 KIND_SCALAR_MODEL= SPECIES_TRANSPORT
-DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY
+%DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY
+DIFFUSIVITY_MODEL= CONSTANT_SCHMIDT
+%DIFFUSIVITY_MODEL= CONSTANT_LEWIS
 SCHMIDT_NUMBER_LAMINAR= 1.0
-DIFFUSIVITY_CONSTANT= 7.56e-5
+CONSTANT_LEWIS_NUMBER= 1.0 1.0
+DIFFUSIVITY_CONSTANT= 5.0e-5
 
 % according to the paper
 SCHMIDT_NUMBER_TURBULENT= 0.7
@@ -152,7 +160,7 @@ CONV_FILENAME= history
 MARKER_ANALYZE= gas_inlet, air_axial_inlet, outlet
 MARKER_ANALYZE_AVERAGE= AREA
 %
-OUTPUT_FILES= RESTART, PARAVIEW_MULTIBLOCK
+OUTPUT_FILES= RESTART, RESTART_ASCII, PARAVIEW_MULTIBLOCK
 VOLUME_OUTPUT= RESIDUAL, PRIMITIVE, SPECIES_UDS_0
 OUTPUT_WRT_FREQ= 100
 %
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/run.py b/TestCases/py_wrapper/turbulent_premixed_psi/run.py
index dc1133169a24..1f355f135e6f 100644
--- a/TestCases/py_wrapper/turbulent_premixed_psi/run.py
+++ b/TestCases/py_wrapper/turbulent_premixed_psi/run.py
@@ -36,6 +36,7 @@
 # without mpi:
 #  comm = 0
 
+Tref = 298.15
 # flame temperature of the methane-air mixture (phi=0.5, P=5)
 Tf = 1777
 
@@ -48,7 +49,7 @@
 # unburnt thermal conductivity of methane-air (phi=0.5, P=5)
 k_u = 0.0523
 # unburnt heat capacity of methane-air (phi=0.5, P=5)
-cp_u = 1311.0
+cp_u = 1350.0
 
 # P = rho*R*T
 # 5 = 2.55 * R * 673
@@ -92,17 +93,23 @@ def update_temperature(SU2Driver, iPoint):
     # Note: returns a list
     C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
     T = Tu*(1-C) + Tf*C
+
+    #iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
+    #solindex = getsolvar(SU2Driver)
+    #iTEMP = solindex.get("TEMPERATURE")
+    #T = SU2Driver.Solution(iFLOWSOLVER)(iPoint,iTEMP)
+    #SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,cp_u*T)
+
+    prim_indices = SU2Driver.GetPrimitiveIndices()
+    iTemp = prim_indices['TEMPERATURE']
+    SU2Driver.Primitives().Set(iPoint,iTemp, T)
+    #ih = prim_indices['ENTHALPY']
+    #SU2Driver.Primitives().Set(iPoint,ih, cp_u*(T-Tref))
+        
     iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
-    # the list with names
-    solindex = getsolvar(SU2Driver)
-    primindex = SU2Driver.GetPrimitiveIndices()
-    # get Cp
-    iCp = primindex.get("CP_TOTAL")
-    Cp = SU2Driver.Primitives()(iPoint,iCp)
-    # get solution index for energy equation
-    iTEMP = solindex.get("TEMPERATURE")
-    # set enthalpy, for Ideal Gas model, enthalpy is Cp*T
-    SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,Cp*T)
+    iENTH = 3
+    #h = 
+    SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iENTH, cp_u*(T-Tref))
 
 
 # ################################################################## #
@@ -216,11 +223,9 @@ def main():
   sys.stdout.flush()
 
   # run N iterations
-  for inner_iter in range(2):
+  for inner_iter in range(10):
     if (rank==0):
       print("python iteration ", inner_iter)
-    driver.Preprocess(inner_iter)
-    driver.Run()
 
     Source = driver.UserDefinedSource(iSPECIESSOLVER)
 
@@ -236,6 +241,10 @@ def main():
       # set the temperature to T = c*Tf + (1-c)*Tu
       update_temperature(driver, i_node)
 
+
+    driver.Preprocess(inner_iter)
+    driver.Run()
+
     driver.Postprocess()
     driver.Update()
     # Monitor the solver and output solution to file if required.
diff --git a/TestCases/radiation/p1model/configp1.cfg b/TestCases/radiation/p1model/configp1.cfg
index 6c779b16d0fb..be5787656d31 100644
--- a/TestCases/radiation/p1model/configp1.cfg
+++ b/TestCases/radiation/p1model/configp1.cfg
@@ -67,7 +67,7 @@ MARKER_MONITORING= ( upper, left )
 LINEAR_SOLVER= FGMRES
 LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ILU_FILL_IN= 0
-LINEAR_SOLVER_ERROR= 1E-8
+LINEAR_SOLVER_ERROR= 1E-6
 LINEAR_SOLVER_ITER= 10
 
 %%%%%%%%%%%%%%%%%%%%%%%
@@ -80,12 +80,14 @@ MUSCL_FLOW= YES
 SLOPE_LIMITER_FLOW= NONE
 TIME_DISCRE_FLOW= EULER_IMPLICIT
 CFL_NUMBER= 100
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.9, 1.10, 100.0, 1000, 0.1, 0)
 
 %%%%%%%%%%%%%%%%%%%%%%%
 % CONVERGENCE CRITERIA
 %%%%%%%%%%%%%%%%%%%%%%%
 
-INNER_ITER= 101
+INNER_ITER= 51
 
 CONV_RESIDUAL_MINVAL= -10
 CONV_STARTITER= 10
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index cfd25784fa03..8dd870bfd7ba 100755
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -407,7 +407,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-4.436657, 0.507847, 0.000000, 0.000000]
+    inc_buoyancy.test_vals = [-3.860463, -3.122402, 4.074509, -23.414000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -458,7 +458,7 @@ def main():
     inc_turb_wallfunction_flatplate_sst.cfg_dir   = "wallfunctions/flatplate/incompressible_SST"
     inc_turb_wallfunction_flatplate_sst.cfg_file  = "turb_SST_flatplate.cfg"
     inc_turb_wallfunction_flatplate_sst.test_iter = 10
-    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.881656, -5.732047, -6.724784, -4.242639, -7.189741, -2.050900, 10.000000, -2.877778, 0.001149, 0.003172, 0.000000]
+    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.881887, -5.731983, -6.725181, -2.618754, -7.189701, -2.050891, 10.000000, -2.027265, 0.001155, 0.003172, 0.000000]
     test_list.append(inc_turb_wallfunction_flatplate_sst)
 
     # FLAT PLATE, WALL FUNCTIONS, INCOMPRESSIBLE SA
@@ -466,7 +466,7 @@ def main():
     inc_turb_wallfunction_flatplate_sa.cfg_dir   = "wallfunctions/flatplate/incompressible_SA"
     inc_turb_wallfunction_flatplate_sa.cfg_file  = "turb_SA_flatplate.cfg"
     inc_turb_wallfunction_flatplate_sa.test_iter = 10
-    inc_turb_wallfunction_flatplate_sa.test_vals = [-6.894226, -5.716031, -6.743814, -4.242551, -9.550079, 10.000000, -2.879412, 0.001021, 0.003759]
+    inc_turb_wallfunction_flatplate_sa.test_vals = [-6.894473, -5.715971, -6.744286, -2.631311, -9.550134, 10.000000, -2.017104, 0.001027, 0.003759]
     test_list.append(inc_turb_wallfunction_flatplate_sa)
 
     ####################
@@ -986,7 +986,7 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals = [19.000000, -1.049468, -1.303013]
+    slinc_steady.test_vals = [19.000000, -1.148249, -1.398402]
     slinc_steady.timeout   = 100
     slinc_steady.multizone = True
     test_list.append(slinc_steady)
@@ -1097,8 +1097,8 @@ def main():
     p1rad           = TestCase('p1rad')
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
-    p1rad.test_iter = 100
-    p1rad.test_vals = [-7.751309, -7.923059, -2.119084, 0.091733, -47.387000]
+    p1rad.test_iter = 50
+    p1rad.test_vals = [-8.284673, -8.008659, -2.424204, 0.389030, -56.560000]
     test_list.append(p1rad)
 
     # ###############################
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 0cb1b5715891..4bffc9ec0aa7 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -194,7 +194,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals = [-2.685957, 0.665986, 0.000000, -9.155800]
+    discadj_heat.test_vals = [-2.677870, 0.674825, 0.000000, -9.215500]
     test_list.append(discadj_heat)
 
     ###################################
@@ -219,7 +219,7 @@ def main():
     discadj_cht.cfg_dir   = "coupled_cht/disc_adj_incomp_2d"
     discadj_cht.cfg_file  = "cht_2d_3cylinders.cfg"
     discadj_cht.test_iter = 10
-    discadj_cht.test_vals = [-6.131112, -2.564634, -2.565220, -2.565304]
+    discadj_cht.test_vals = [-6.308199, -3.086023, -3.086024, -3.086024]
     test_list.append(discadj_cht)
 
     ######################################
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
index 0768bc29370d..eccb5f3e630c 100644
--- a/TestCases/tutorials.py
+++ b/TestCases/tutorials.py
@@ -49,7 +49,7 @@ def main():
     cht_incompressible_unsteady.cfg_dir   = "../Tutorials/multiphysics/unsteady_cht/"
     cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
     cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-2.661440, -2.534489, -0.080399, -0.080399, -0.080399, -12.421979, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000] #last columns
+    cht_incompressible_unsteady.test_vals = [-1.742894, -3.075372, -0.080399, -0.080399, -0.080399, -11.163219, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000] #last columns
     cht_incompressible_unsteady.multizone = True
     cht_incompressible_unsteady.unsteady  = True
     test_list.append(cht_incompressible_unsteady)
@@ -69,7 +69,7 @@ def main():
     cht_CR.cfg_dir   = "../Tutorials/multiphysics/contact_resistance_cht"
     cht_CR.cfg_file  = "master.cfg"
     cht_CR.test_iter = 80
-    cht_CR.test_vals = [-8.606916, -9.227614, -10.411674, -2.116659]
+    cht_CR.test_vals = [-8.606938, -9.227901, -10.410107, -2.129377]
     cht_CR.multizone = True
     test_list.append(cht_CR)
 
@@ -80,7 +80,7 @@ def main():
     sp_pinArray_2d_mf_hf.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Streamwise_Periodic"
     sp_pinArray_2d_mf_hf.cfg_file  = "sp_pinArray_2d_mf_hf.cfg"
     sp_pinArray_2d_mf_hf.test_iter = 25
-    sp_pinArray_2d_mf_hf.test_vals = [-4.685131, 1.388616, -0.755475, 241.878084]
+    sp_pinArray_2d_mf_hf.test_vals = [-4.684477, 0.844369, -0.755565, 241.872175]
     sp_pinArray_2d_mf_hf.test_vals_aarch64 = [-4.686092, 1.387918, -0.755447, 241.878841]
     test_list.append(sp_pinArray_2d_mf_hf)
 
@@ -89,7 +89,7 @@ def main():
     sp_pinArray_2d_dp_hf_tp.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Streamwise_Periodic"
     sp_pinArray_2d_dp_hf_tp.cfg_file  = "sp_pinArray_2d_dp_hf_tp.cfg"
     sp_pinArray_2d_dp_hf_tp.test_iter = 25
-    sp_pinArray_2d_dp_hf_tp.test_vals = [-4.737196, 1.322378, -0.713373, 208.023676]
+    sp_pinArray_2d_dp_hf_tp.test_vals = [-4.739886, 0.520841, -0.713547, 208.023676]
     sp_pinArray_2d_dp_hf_tp.test_vals_aarch64 = [-4.733643, 1.325195, -0.713411, 208.023676]
     test_list.append(sp_pinArray_2d_dp_hf_tp)
 
@@ -136,7 +136,7 @@ def main():
     species3_primitiveVenturi.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Species_Transport"
     species3_primitiveVenturi.cfg_file  = "species3_primitiveVenturi.cfg"
     species3_primitiveVenturi.test_iter = 50
-    species3_primitiveVenturi.test_vals = [-5.673008, -4.673131, -4.692722, -5.445021, -1.069687, -5.960203, -6.065022, 5.000000, -0.570021, 5.000000, -2.512831, 5.000000, -0.525792, 1.660602, 0.502278, 0.603347, 0.554978]
+    species3_primitiveVenturi.test_vals = [-5.673011, -4.673114, -4.692686, -5.445043, -1.069650, -5.960206, -6.065010, 5.000000, -0.569992, 5.000000, -2.512966, 5.000000, -0.525813, 1.660601, 0.502277, 0.603346, 0.554977]
     test_list.append(species3_primitiveVenturi)
 
     # 3 species (2 eq) primitive venturi mixing
@@ -144,7 +144,7 @@ def main():
     DAspecies3_primitiveVenturi.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Species_Transport"
     DAspecies3_primitiveVenturi.cfg_file  = "DAspecies3_primitiveVenturi.cfg"
     DAspecies3_primitiveVenturi.test_iter = 50
-    DAspecies3_primitiveVenturi.test_vals = [-9.819097, -8.643457, -8.676919, -8.347340, -12.926241, -9.739487, -8.947991]
+    DAspecies3_primitiveVenturi.test_vals = [-9.821887, -8.644118, -8.677534, -8.347228, -12.926295, -9.739487, -8.947991]
     DAspecies3_primitiveVenturi.command   = TestCase.Command("mpirun -n 2", "SU2_CFD_AD")
     test_list.append(DAspecies3_primitiveVenturi)
 
@@ -153,7 +153,7 @@ def main():
     kenics_mixer_tutorial.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Species_Transport_Composition_Dependent_Model"
     kenics_mixer_tutorial.cfg_file  = "kenics_mixer_tutorial.cfg"
     kenics_mixer_tutorial.test_iter = 10
-    kenics_mixer_tutorial.test_vals = [-7.490303, -6.823796, -6.838369, -6.383332, -7.907780, -3.062276, -7.451176, 5.000000, -1.858354, 4.000000, -5.318193, 3.000000, -6.371744, 0.025671, 0.000000, 0.025671, 0.000000, 62.846000, 8.470600, 46.847000, 7.527900]
+    kenics_mixer_tutorial.test_vals = [-7.490472, -6.823803, -6.838373, -6.383764, -7.907928, -3.062404, -7.452184, 5.000000, -1.858094, 4.000000, -5.318066, 3.000000, -6.371967, 0.025661, 0.000000, 0.025661, 0.000000, 62.736000, 8.470600, 46.738000, 7.527900]
     kenics_mixer_tutorial.command   = TestCase.Command("mpirun -n 2", "SU2_CFD")
     test_list.append(kenics_mixer_tutorial)
 
diff --git a/TestCases/vandv.py b/TestCases/vandv.py
index f7b197632ca4..a379e775c0f0 100644
--- a/TestCases/vandv.py
+++ b/TestCases/vandv.py
@@ -111,7 +111,7 @@ def main():
     sandiajet_sst.cfg_dir   = "vandv/species_transport/sandia_jet"
     sandiajet_sst.cfg_file  = "validation.cfg"
     sandiajet_sst.test_iter = 5
-    sandiajet_sst.test_vals = [-8.271400, -6.305945, -6.123143, -5.634565, -1.687556, -7.869775, 5.000000, -1.527449, 5.000000, -4.053788, 5.000000, -2.303502, 0.000257, 0.000000, 0.000000, 0.000257, 4020.500000, 3919.900000, 49.151000, 51.436000]
+    sandiajet_sst.test_vals = [-4.930056, -3.072757, -3.531917, -3.632319, 1.186060, -6.175518, 5.000000, -1.413277, 5.000000, -4.735257, 5.000000, -1.924560, 0.000256, 0.000000, 0.000000, 0.000256, 4126.200000, 4008.100000, 62.399000, 55.713000]
     sandiajet_sst.test_vals_aarch64 = [-8.271400, -6.305945, -6.123143, -5.634565, -1.687556, -7.869775, 5.000000, -1.527449, 5.000000, -4.053788, 5.000000, -2.303502, 0.000257, 0.000000, 0.000000, 0.000257, 4020.500000, 3919.900000, 49.151000, 51.436000]
     test_list.append(sandiajet_sst)