diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp
index 9c24b9cfed5..069911e11a1 100644
--- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp
+++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp
@@ -29,6 +29,7 @@
 
 #include "CVariable.hpp"
 #include "../fluid/CNEMOGas.hpp"
+#include "../../Common/include/toolboxes/geometry_toolbox.hpp"
 
 /*!
  * \class CNEMOEulerVariable
diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
index 4f88dee8d14..ee000e94f5a 100644
--- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
@@ -46,10 +46,9 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
     description = "Euler";
   }
 
-  unsigned long iPoint, counter_local, counter_global = 0;
-  unsigned short iDim, iMarker, iSpecies, nLineLets;
+  unsigned short iMarker, nLineLets;
   unsigned short nZone = geometry->GetnZone();
-  su2double *Mvec_Inf, Alpha, Beta, Soundspeed_Inf, sqvel;
+  su2double *Mvec_Inf, Alpha, Beta;
   bool restart   = (config->GetRestart() || config->GetRestart_Flow());
   unsigned short direct_diff = config->GetDirectDiff();
   int Unst_RestartIter = 0;
@@ -59,8 +58,6 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
   bool adjoint = config->GetDiscrete_Adjoint();
   string filename_ = "flow";
 
-  bool nonPhys;
-
   /*--- Store the multigrid level. ---*/
   MGLevel = iMesh;
 
@@ -216,52 +213,6 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
 
   node_infty->SetPrimVar(0, FluidModel);
 
-  /*--- Check that the initial solution is physical, report any non-physical nodes ---*/
-
-  counter_local = 0;
-
-  for (iPoint = 0; iPoint < nPoint; iPoint++) {
-
-    nonPhys = nodes->SetPrimVar(iPoint, FluidModel);
-
-    /*--- Set mixture state ---*/
-    FluidModel->SetTDStatePTTv(Pressure_Inf, MassFrac_Inf, Temperature_Inf, Temperature_ve_Inf);
-
-    /*--- Compute other freestream quantities ---*/
-    Density_Inf    = FluidModel->GetDensity();
-    Soundspeed_Inf = FluidModel->GetSoundSpeed();
-
-    sqvel = 0.0;
-    for (iDim = 0; iDim < nDim; iDim++){
-      sqvel += Mvec_Inf[iDim]*Soundspeed_Inf * Mvec_Inf[iDim]*Soundspeed_Inf;
-    }
-    const auto& Energies_Inf = FluidModel->ComputeMixtureEnergies();
-
-    /*--- Initialize Solution & Solution_Old vectors ---*/
-    for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-      Solution[iSpecies]      = Density_Inf*MassFrac_Inf[iSpecies];
-    }
-    for (iDim = 0; iDim < nDim; iDim++) {
-      Solution[nSpecies+iDim] = Density_Inf*Mvec_Inf[iDim]*Soundspeed_Inf;
-    }
-    Solution[nSpecies+nDim]     = Density_Inf*(Energies_Inf[0] + 0.5*sqvel);
-    Solution[nSpecies+nDim+1]   = Density_Inf*Energies_Inf[1];
-    nodes->SetSolution(iPoint,Solution);
-    nodes->SetSolution_Old(iPoint,Solution);
-
-    if(nonPhys)
-      counter_local++;
-  }
-
-  /*--- Warning message about non-physical points ---*/
-  if (config->GetComm_Level() == COMM_FULL) {
-
-    SU2_MPI::Reduce(&counter_local, &counter_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());
-
-    if ((rank == MASTER_NODE) && (counter_global != 0))
-      cout << "Warning. The original solution contains "<< counter_global << " points that are not physical." << endl;
-  }
-
   /*--- Initial comms. ---*/
 
   CommunicateInitialState(geometry, config);
@@ -303,7 +254,9 @@ void CNEMOEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver
       unsigned long tmp = ErrorCounter;
       SU2_MPI::Allreduce(&tmp, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       config->SetNonphysical_Points(ErrorCounter);
-
+      
+      if ((rank == MASTER_NODE) && (ErrorCounter != 0))
+        cout << "Warning. The initial solution contains "<< ErrorCounter << " points that are not physical." << endl;
   }
 
   /*--- Artificial dissipation ---*/
diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp
index 310b1ea2b7e..3e95b1afca1 100644
--- a/SU2_CFD/src/solvers/CSolverFactory.cpp
+++ b/SU2_CFD/src/solvers/CSolverFactory.cpp
@@ -429,7 +429,6 @@ CSolver* CSolverFactory::CreateFlowSolver(SUB_SOLVER_TYPE kindFlowSolver, CSolve
       break;
     case SUB_SOLVER_TYPE::NEMO_EULER:
       flowSolver = new CNEMOEulerSolver(geometry, config, iMGLevel);
-      flowSolver->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false);
       break;
     case SUB_SOLVER_TYPE::NEMO_NAVIER_STOKES:
       flowSolver = new CNEMONSSolver(geometry, config, iMGLevel);
diff --git a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp
index f9b73aaf6fd..38a77d2c585 100644
--- a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp
+++ b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp
@@ -45,9 +45,7 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
                                                                          config ),
                                        Gradient_Reconstruction(config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient_Primitive) {
 
-  vector<su2double> energies;
   unsigned short iDim, iSpecies;
-  su2double soundspeed, sqvel, rho;
 
   /*--- Setting variable amounts ---*/
   nDim            = ndim;
@@ -143,23 +141,18 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
   UnderRelaxation.resize(nPoint) = su2double(1.0);
   LocalCFL.resize(nPoint) = su2double(0.0);
 
+  /*--- Set mixture state ---*/
+  fluidmodel->SetTDStatePTTv(val_pressure, val_massfrac, val_temperature, val_temperature_ve);
+
+  /*--- Compute necessary quantities ---*/
+  const su2double rho = fluidmodel->GetDensity();
+  const su2double soundspeed = fluidmodel->ComputeSoundSpeed();
+  const su2double sqvel = GeometryToolbox::SquaredNorm(nDim, val_mach) * pow(soundspeed,2);
+  const auto& energies = fluidmodel->ComputeMixtureEnergies();
+  
   /*--- Loop over all points --*/
   for(unsigned long iPoint = 0; iPoint < nPoint; ++iPoint){
 
-    /*--- Reset velocity^2 [m2/s2] to zero ---*/
-    sqvel = 0.0;
-
-    /*--- Set mixture state ---*/
-    fluidmodel->SetTDStatePTTv(val_pressure, val_massfrac, val_temperature, val_temperature_ve);
-
-    /*--- Compute necessary quantities ---*/
-    rho = fluidmodel->GetDensity();
-    soundspeed = fluidmodel->ComputeSoundSpeed();
-    for (iDim = 0; iDim < nDim; iDim++){
-      sqvel += val_mach[iDim]*soundspeed * val_mach[iDim]*soundspeed;
-    }
-    energies = fluidmodel->ComputeMixtureEnergies();
-
     /*--- Initialize Solution & Solution_Old vectors ---*/
     for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
       Solution(iPoint,iSpecies)     = rho*val_massfrac[iSpecies];
@@ -168,14 +161,9 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
 
     Solution(iPoint,nSpecies+nDim)       = rho*(energies[0]+0.5*sqvel);
     Solution(iPoint,nSpecies+nDim+1)     = rho*(energies[1]);
-
-    Solution_Old = Solution;
-
-    /*--- Assign primitive variables ---*/
-    Primitive(iPoint,T_INDEX)   = val_temperature;
-    Primitive(iPoint,TVE_INDEX) = val_temperature_ve;
-    Primitive(iPoint,P_INDEX)   = val_pressure;
   }
+
+  Solution_Old = Solution;
 }
 
 void CNEMOEulerVariable::SetVelocity2(unsigned long iPoint) {