diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index bf99257b1387..9124405da691 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1062,7 +1062,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   const bool sst_compWilcox = IsPresent(SST_OPTIONS::COMP_Wilcox);
   const bool sst_compSarkar = IsPresent(SST_OPTIONS::COMP_Sarkar);
   const bool sst_dll = IsPresent(SST_OPTIONS::DLL);
-
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
   } else if (sst_2003) {
@@ -1473,10 +1472,16 @@ static const MapType<std::string, WINDOW_FUNCTION> Window_Map = {
  */
 enum ENUM_HYBRIDRANSLES {
   NO_HYBRIDRANSLES = 0,  /*!< \brief No turbulence model. */
-  SA_DES   = 1,          /*!< \brief Kind of Hybrid RANS/LES (SA - Detached Eddy Simulation (DES)). */
-  SA_DDES  = 2,          /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Delta_max SGS ). */
-  SA_ZDES  = 3,          /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Vorticity based SGS like Zonal DES). */
-  SA_EDDES = 4           /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */
+  SA_DES           = 1,  /*!< \brief Kind of Hybrid RANS/LES (SA - Detached Eddy Simulation (DES)). */
+  SA_DDES          = 2,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Delta_max SGS ). */
+  SA_ZDES          = 3,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Vorticity based SGS like Zonal DES). */
+  SA_EDDES         = 4,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */
+  SA_EDDES_UNSTR   = 5,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES and for Unstructured grids). */
+  SST_DDES         = 6,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): DDES). */
+  SST_IDDES        = 7,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */
+  SST_SIDDES       = 8,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */
+  SST_EDDES        = 9,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
+  SST_EDDES_UNSTR  = 10  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
 };
 static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("NONE", NO_HYBRIDRANSLES)
@@ -1484,6 +1489,12 @@ static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("SA_DDES", SA_DDES)
   MakePair("SA_ZDES", SA_ZDES)
   MakePair("SA_EDDES", SA_EDDES)
+  MakePair("SA_EDDES_UNSTR", SA_EDDES_UNSTR)
+  MakePair("SST_DDES", SST_DDES)
+  MakePair("SST_IDDES", SST_IDDES)
+  MakePair("SST_SIDDES", SST_SIDDES)
+  MakePair("SST_EDDES", SST_EDDES)
+  MakePair("SST_EDDES_UNSTR", SST_EDDES_UNSTR)
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 0a9cca0f63f7..91bfec3e81cf 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -6469,10 +6469,14 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         cout << "Hybrid RANS/LES: ";
         switch (Kind_HybridRANSLES) {
           case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
-          case SA_DES:   cout << "Detached Eddy Simulation (DES97) " << endl; break;
-          case SA_DDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
-          case SA_ZDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
-          case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
+          case SA_DES:         cout << "Detached Eddy Simulation (DES97) " << endl; break;
+          case SA_DDES:        cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
+          case SA_ZDES:        cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
+          case SA_EDDES:       cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
+          case SA_EDDES_UNSTR: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS (modified for Unstructured grids)" << endl; break;
+          case SST_DDES:       cout << "Delayed Detached Eddy Simulation (DDES)" << endl; break;
+          case SST_IDDES:      cout << "Improved Delayed Detached Eddy Simulation (IDDES)" << endl; break;
+          case SST_SIDDES:     cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES)" << endl; break;
         }
         break;
       case MAIN_SOLVER::NEMO_EULER:
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index d5fb6a0ec11e..ba2ad7ec751f 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -135,6 +135,7 @@ class CNumerics {
   const su2double
   *TurbPsi_i,  /*!< \brief Vector of adjoint turbulent variables at point i. */
   *TurbPsi_j;  /*!< \brief Vector of adjoint turbulent variables at point j. */
+  su2double lengthScale_i, lengthScale_j; /*!< \brief length scale for SST */
   CMatrixView<const su2double>
   ConsVar_Grad_i,  /*!< \brief Gradient of conservative variables at point i. */
   ConsVar_Grad_j,  /*!< \brief Gradient of conservative variables at point j. */
@@ -827,6 +828,16 @@ class CNumerics {
     dist_j = val_dist_j;
   }
 
+  /*!
+   * \brief Set the value of the length scale for SST.
+   * \param[in] val_lengthScale_i - Value of of the length scale for SST from point i.
+   * \param[in] val_lengthScale_j - Value of of the length scale for SST from point j.
+   */
+  void SetLengthScale(su2double val_lengthScale_i, su2double val_lengthScale_j) {
+    lengthScale_i = val_lengthScale_i;
+    lengthScale_j = val_lengthScale_j;
+  }
+
   /*!
    * \brief Set the value of the roughness from the nearest wall.
    * \param[in] val_dist_i - Value of of the roughness of the nearest wall from point i
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 2195f94bd21f..afadc4bd1282 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -885,14 +885,17 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       }
 
       if (sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) {
-        const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt;
+        const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * pow(Mt, 2);
         pk += Dilatation_Sarkar;
       }
 
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
-      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
+      if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES)
+        dk = Density_i * sqrt(pow(ScalarVar_i[0], 3)) / lengthScale_i;
+        
+      su2double dw = beta_blended * Density_i * pow(ScalarVar_i[1], 2) * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
       if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {
diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
index 707316f06094..9abb3d0115f4 100644
--- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
+++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
@@ -61,6 +61,17 @@ class CTurbSSTSolver final : public CTurbSolver {
    */
   void ComputeUnderRelaxationFactor(const CConfig *config);
 
+
+   /*!
+   * \brief A virtual member.
+   * \param[in] solver - Solver container
+   * \param[in] geometry - Geometrical definition.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetDES_LengthScale(CSolver** solver,
+                          CGeometry *geometry,
+                          CConfig *config);
+
 public:
   /*!
    * \brief Constructor.
diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp
index adf3cb5cf17b..6a533d4fb28b 100644
--- a/SU2_CFD/include/variables/CTurbSAVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp
@@ -39,8 +39,6 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
-  VectorType Vortex_Tilting;
 
 public:
   /*!
@@ -60,31 +58,4 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   ~CTurbSAVariable() override = default;
 
-  /*!
-   * \brief Get the DES length scale
-   * \param[in] iPoint - Point index.
-   * \return Value of the DES length Scale.
-   */
-  inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); }
-
-  /*!
-   * \brief Set the DES Length Scale.
-   * \param[in] iPoint - Point index.
-   */
-  inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; }
-
-  /*!
-   * \brief Set the vortex tilting measure for computation of the EDDES length scale
-   * \param[in] iPoint - Point index.
-   */
-  void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double>,
-                         const su2double* Vorticity, su2double LaminarViscosity) override;
-
-  /*!
-   * \brief Get the vortex tilting measure for computation of the EDDES length scale
-   * \param[in] iPoint - Point index.
-   * \return Value of the DES length Scale
-   */
-  inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
-
 };
diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp
index c0c33efd8c90..2bb514786223 100644
--- a/SU2_CFD/include/variables/CTurbVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbVariable.hpp
@@ -37,12 +37,14 @@
  */
 class CTurbVariable : public CScalarVariable {
 protected:
-  VectorType muT; /*!< \brief Eddy viscosity. */
+  VectorType muT;             /*!< \brief Eddy viscosity. */
+  VectorType DES_LengthScale; /*!< \brief DES Length scale. */
 
 public:
   static constexpr size_t MAXNVAR = 2;
-  VectorType turb_index;
-  VectorType intermittency;         /*!< \brief Value of the intermittency for the trans. model. */
+  VectorType turb_index;            /*!< \brief Value of the turbulence index for transition simulations. */
+  VectorType intermittency;         /*!< \brief Value of the intermittency for the transition model. */
+  VectorType Vortex_Tilting;        /*!< \brief Value of the vortex tilting measure for EDDES models. */
 
   /*!
    * \brief Constructor of the class.
@@ -100,6 +102,32 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \brief Get the DES length scale
+   * \param[in] iPoint - Point index.
+   * \return Value of the DES length Scale.
+   */
+  inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); }
+
+  /*!
+   * \brief Set the DES Length Scale.
+   * \param[in] iPoint - Point index.
+   */
+  inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; }
+
+  /*!
+   * \brief Set the vortex tilting measure for computation of the EDDES length scale
+   * \param[in] iPoint - Point index.
+   */
+  void SetVortex_Tilting(unsigned long iPoint, su2double **Strain,
+                         const su2double* Vorticity, su2double LaminarViscosity) override;
+
+  /*!
+   * \brief Get the vortex tilting measure for computation of the EDDES length scale
+   * \param[in] iPoint - Point index.
+   * \return Value of the DES length Scale
+   */
+  inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
   /*!
    * \brief Set the Diffusion Coefficients of TKE and omega equations.
    * \param[in] iPoint - Point index.
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index dc905d79fcb4..a19453754bf7 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -2235,7 +2235,7 @@ class CVariable {
 
   inline virtual su2double GetTau_Wall(unsigned long iPoint) const { return 0.0; }
 
-  inline virtual void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
+  inline virtual void SetVortex_Tilting(unsigned long iPoint, su2double **Strain,
                                         const su2double* Vorticity, su2double LaminarViscosity) {}
 
   inline virtual su2double GetVortex_Tilting(unsigned long iPoint) const { return 0.0; }
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index b355c0698afb..e087a1fd61a1 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1483,6 +1483,9 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value");
     AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value");
+    AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations");
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
+      AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "DDES", "desired grid size for Scale Resolving Simulations");
   }
 
   if (config->GetViscous()) {
@@ -1494,6 +1497,27 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
       AddVolumeOutput("VORTICITY", "Vorticity", "VORTEX_IDENTIFICATION", "Value of the vorticity");
     }
     AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion");
+    if (config->GetKind_HybridRANSLES() == NO_HYBRIDRANSLES) {
+      AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DEBUG", "Wall distance value");
+    }
+    AddVolumeOutput("GRAD_VEL_XX", "Grad_Vel_xx", "VELOCITY_GRADIENT", "X-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_XY", "Grad_Vel_xy", "VELOCITY_GRADIENT", "Y-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_YX", "Grad_Vel_yx", "VELOCITY_GRADIENT", "X-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_YY", "Grad_Vel_yy", "VELOCITY_GRADIENT", "Y-Gradient of V");
+    if (nDim == 3) {
+      AddVolumeOutput("GRAD_VEL_XZ", "Grad_Vel_xz", "VELOCITY_GRADIENT", "Z-Gradient of U");
+      AddVolumeOutput("GRAD_VEL_YZ", "Grad_Vel_yz", "VELOCITY_GRADIENT", "Z-Gradient of V");
+      AddVolumeOutput("GRAD_VEL_ZX", "Grad_Vel_zx", "VELOCITY_GRADIENT", "X-Gradient of W");
+      AddVolumeOutput("GRAD_VEL_ZY", "Grad_Vel_zy", "VELOCITY_GRADIENT", "Y-Gradient of W");
+      AddVolumeOutput("GRAD_VEL_ZZ", "Grad_Vel_zz", "VELOCITY_GRADIENT", "Z-Gradient of W");
+    }
+
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST){
+      AddVolumeOutput("CDkw", "CDkw", "SST_QUANTITIES", "Cross-Diffusion term");
+      AddVolumeOutput("F1", "F1", "SST_QUANTITIES", "F1 blending function");
+      AddVolumeOutput("F2", "F2", "SST_QUANTITIES", "F2 blending function");
+    }
+
   }
 
   // Timestep info
@@ -1519,6 +1543,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
   SetVolumeOutputValue("CFL", iPoint, Node_Flow->GetLocalCFL(iPoint));
 
   if (config->GetViscous()) {
+    const auto VelGrad = Node_Flow->GetVelocityGradient(iPoint);
     if (nDim == 3){
       SetVolumeOutputValue("VORTICITY_X", iPoint, Node_Flow->GetVorticity(iPoint)[0]);
       SetVolumeOutputValue("VORTICITY_Y", iPoint, Node_Flow->GetVorticity(iPoint)[1]);
@@ -1526,7 +1551,27 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     } else {
       SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     }
-    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint)));
+    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(VelGrad));
+    SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+
+    SetVolumeOutputValue("GRAD_VEL_XX", iPoint, VelGrad(0,0));
+    SetVolumeOutputValue("GRAD_VEL_XY", iPoint, VelGrad(0,1));
+    SetVolumeOutputValue("GRAD_VEL_YX", iPoint, VelGrad(1,0));
+    SetVolumeOutputValue("GRAD_VEL_YY", iPoint, VelGrad(1,1));
+    if (nDim == 3) {
+      SetVolumeOutputValue("GRAD_VEL_XZ", iPoint, VelGrad(0,2));
+      SetVolumeOutputValue("GRAD_VEL_YZ", iPoint, VelGrad(1,2));
+      SetVolumeOutputValue("GRAD_VEL_ZX", iPoint, VelGrad(2,0));
+      SetVolumeOutputValue("GRAD_VEL_ZY", iPoint, VelGrad(2,1));
+      SetVolumeOutputValue("GRAD_VEL_ZZ", iPoint, VelGrad(2,2));
+    }
+    
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST){
+      SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint));
+      SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint));
+      SetVolumeOutputValue("F2", iPoint, Node_Turb->GetF2blending(iPoint));
+    }
+
   }
 
   const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE);
@@ -1582,6 +1627,17 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint));
     SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+    const su2double mut = Node_Flow->GetEddyViscosity(iPoint);
+    const su2double mu = Node_Flow->GetLaminarViscosity(iPoint); 
+    const su2double LESIQ = 1.0/(1.0+0.05*pow((mut+mu)/mu, 0.53));
+    SetVolumeOutputValue("LESIQ", iPoint, LESIQ);
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST) {
+      const su2double betaStar = 0.09; // constants[6]
+      const su2double RANSLength = sqrt(Node_Flow->GetSolution(iPoint, 0)) / max(1e-20, (betaStar * Node_Flow->GetSolution(iPoint, 1)));
+      const su2double RatioL = 0.1;  // TODO:: it should be less or equal than 0.2 - 0.1. Should be taken as input from config?
+      const su2double SRSGridSize = RANSLength * RatioL;
+      SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, SRSGridSize);
+    }
   }
 
   switch (config->GetKind_Species_Model()) {
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index a9c3875406bf..d630213df6f6 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -968,7 +968,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             }
 
             break;
-
+            
           default:
             SU2_MPI::Error("Unrecognized quantity for periodic communication.",
                            CURRENT_FUNCTION);
@@ -1290,7 +1290,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
                 limiter(iPoint, iVar) = min(limiter(iPoint, iVar), bufDRecv[buf_offset+iVar]);
 
               break;
-
+              
             default:
 
               SU2_MPI::Error("Unrecognized quantity for periodic communication.",
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 7649165d0be3..689675422beb 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -186,15 +186,24 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
 
     /*--- Set the vortex tilting coefficient at every node if required ---*/
 
-    if (kind_hybridRANSLES == SA_EDDES){
+    if (kind_hybridRANSLES == SA_EDDES || kind_hybridRANSLES == SA_EDDES_UNSTR){
       auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+        su2double **Grad_Vel = new su2double* [nDim];
+        su2double **StrainMat = new su2double* [nDim];
         auto Vorticity = flowNodes->GetVorticity(iPoint);
-        auto PrimGrad_Flow = flowNodes->GetGradient_Primitive(iPoint);
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+          Grad_Vel[iDim] = new su2double [nDim];
+          StrainMat[iDim] = new su2double [nDim];
+          for (unsigned short jDim = 0; jDim < nDim; jDim++) {
+            Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
+          }
+        }
         auto Laminar_Viscosity = flowNodes->GetLaminarViscosity(iPoint);
-        nodes->SetVortex_Tilting(iPoint, PrimGrad_Flow, Vorticity, Laminar_Viscosity);
+        CNumerics::ComputeMeanRateOfStrainMatrix(3, StrainMat, Grad_Vel);
+        nodes->SetVortex_Tilting(iPoint, StrainMat, Vorticity, Laminar_Viscosity);
       }
       END_SU2_OMP_FOR
     }
@@ -1419,8 +1428,8 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
 
         const su2double maxDelta = geometry->nodes->GetMaxLength(iPoint);
 
-        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
-        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
+        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
+        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
 
         const su2double distDES = constDES * maxDelta;
         lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
@@ -1455,8 +1464,8 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
                                   pow(ratioOmega[1], 2)*delta[0]*delta[2] +
                                   pow(ratioOmega[2], 2)*delta[0]*delta[1]);
 
-        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
-        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
+        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
+        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
 
         if (f_d < 0.99){
           maxDelta = deltaDDES;
@@ -1479,7 +1488,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
 
         su2double ratioOmega[MAXNDIM] = {};
 
-        for (auto iDim = 0; iDim < 3; iDim++){
+        for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
           ratioOmega[iDim] = vorticity[iDim]/omega;
         }
 
@@ -1506,8 +1515,64 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
                                    min(f_max,
                                        f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
 
-        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
-        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
+        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
+        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
+
+        su2double maxDelta = (ln_max/sqrt(3.0)) * f_kh;
+        if (f_d < 0.999){
+          maxDelta = deltaDDES;
+        }
+
+        const su2double distDES = constDES * maxDelta;
+        lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
+
+        break;
+      }
+      case SA_EDDES_UNSTR: {
+        /*--- An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows.
+         Shur et al.
+         Flow Turbulence Combust - 2015
+         ---*/
+
+        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);
+
+        const su2double omega = GeometryToolbox::Norm(3, vorticity);
+
+        su2double ratioOmega[MAXNDIM] = {};
+        for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
+          ratioOmega[iDim] = vorticity[iDim]/omega;
+        }
+
+        const su2double deltaDDES = geometry->nodes->GetMaxLength(iPoint);
+        
+        // Introduced correction for unstructured grids
+        su2double ln_max = -1;
+        for (const auto jPoint : geometry->nodes->GetPoints(iPoint)){
+          const auto coord_j = geometry->nodes->GetCoord(jPoint);
+
+          for (const auto kPoint : geometry->nodes->GetPoints(iPoint)){
+            const auto coord_k = geometry->nodes->GetCoord(kPoint);
+
+            su2double delta[MAXNDIM] = {};
+            for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
+              // TODO: Should I divide by 2 as I am interested in the dual volume (the edge is split at midpoint)?
+              delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0; 
+            }
+            su2double l_n_minus_m[MAXNDIM];
+            GeometryToolbox::CrossProduct(delta, ratioOmega, l_n_minus_m);
+            ln_max = max(ln_max, GeometryToolbox::Norm(3, l_n_minus_m));
+          }
+          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
+        }
+        
+        vortexTiltingMeasure /= (nNeigh + 1);
+
+        const su2double f_kh = max(f_min,
+                                   min(f_max,
+                                       f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
+
+        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2));
+        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3));
 
         su2double maxDelta = (ln_max/sqrt(3.0)) * f_kh;
         if (f_d < 0.999){
diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
index b36f5c2b8765..03e2e012d3ca 100644
--- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
@@ -199,8 +199,43 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain
          unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) {
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
+  const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES();
+
   /*--- Upwind second order reconstruction and gradients ---*/
   CommonPreprocessing(geometry, config, Output);
+
+  if (kind_hybridRANSLES != NO_HYBRIDRANSLES) {
+
+    /*--- Set the vortex tilting coefficient at every node if required ---*/
+
+    if (kind_hybridRANSLES == SST_EDDES || kind_hybridRANSLES == SST_EDDES_UNSTR){
+      auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+
+      SU2_OMP_FOR_STAT(omp_chunk_size)
+      for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+        su2double **Grad_Vel = new su2double* [nDim];
+        su2double **StrainMat = new su2double* [nDim];
+        auto Vorticity = flowNodes->GetVorticity(iPoint);
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+          Grad_Vel[iDim] = new su2double [nDim];
+          StrainMat[iDim] = new su2double [nDim];
+          for (unsigned short jDim = 0; jDim < nDim; jDim++) {
+            Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
+          }
+        }
+        auto Laminar_Viscosity = flowNodes->GetLaminarViscosity(iPoint);
+        CNumerics::ComputeMeanRateOfStrainMatrix(3, StrainMat, Grad_Vel);
+        nodes->SetVortex_Tilting(iPoint, StrainMat, Vorticity, Laminar_Viscosity);
+      }
+      END_SU2_OMP_FOR
+    }
+
+    /*--- Compute the DES length scale ---*/
+
+    SetDES_LengthScale(solver_container, geometry, config);
+
+  }
+
 }
 
 void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container,
@@ -369,12 +404,16 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta
       numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
     }
 
+    if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
+      numerics->SetLengthScale(nodes->GetDES_LengthScale(iPoint), 0.0);
+    }
+
     /*--- Compute the source term ---*/
 
     auto residual = numerics->ComputeResidual(config);
 
     /*--- Store the intermittency ---*/
-
+    
     if (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) {
       nodes->SetIntermittency(iPoint, numerics->GetIntermittencyEff());
     }
@@ -1027,6 +1066,230 @@ void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMark
 
 }
 
+
+void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CConfig *config){
+
+  const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES();
+
+  auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver[FLOW_SOL]->GetNodes());
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+
+    const auto coord_i       = geometry->nodes->GetCoord(iPoint);
+    const auto nNeigh        = geometry->nodes->GetnPoint(iPoint);
+
+    const su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12);
+    const auto Vorticity      = flowNodes->GetVorticity(iPoint);
+    const su2double VortMag   = max(GeometryToolbox::Norm(3, Vorticity), 1e-12);
+
+    const su2double KolmConst2 = 0.41*0.41;
+    const su2double wallDist2 = geometry->nodes->GetWall_Distance(iPoint)*geometry->nodes->GetWall_Distance(iPoint); 
+    
+    const su2double eddyVisc = nodes->GetmuT(iPoint)/flowNodes->GetDensity(iPoint);
+    const su2double lamVisc = nodes->GetLaminarViscosity(iPoint)/flowNodes->GetDensity(iPoint);
+
+    const su2double C_DES1 = 0.78;
+    const su2double C_DES2 = 0.61;
+
+    const su2double h_max = geometry->nodes->GetMaxLength(iPoint);
+    const su2double C_DES = C_DES1 * nodes->GetF1blending(iPoint) + C_DES2 * (1-nodes->GetF1blending(iPoint));
+    const su2double l_RANS = sqrt(nodes->GetSolution(iPoint, 0)) / (constants[6] * nodes->GetSolution(iPoint, 1));
+
+    su2double DES_lengthScale = 0.0;
+
+    switch(kind_hybridRANSLES){
+      /*--- Every model is taken from "Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model"
+                                      Mikhail S. Gritskevich et al. (DOI:10.1007/s10494-011-9378-4)
+        ---*/
+      case SST_DDES: {
+
+        const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double C_d1 = 20.0;
+        const su2double C_d2 = 3.0;
+
+        const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2));
+
+        const su2double l_LES = C_DES * h_max;
+        
+        DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES);
+
+        break;
+      }
+      case SST_IDDES: {
+        
+        // Constants
+        const su2double C_w = 0.15;
+        const su2double C_dt1 = 20.0;
+        const su2double C_dt2 = 3.0;
+        const su2double C_l = 5.0;
+        const su2double C_t = 1.87;
+
+        const su2double alpha = 0.25 - sqrt(wallDist2) / h_max;
+        const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0);
+        const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2));
+        const su2double ftilda_d = max(1.0 - f_dt, f_b);
+
+        const su2double r_dl = lamVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double f_l = tanh(pow(C_l*C_l*r_dl, 10.0));
+        const su2double f_t = tanh(pow(C_t*C_t*r_dt, 3.0));
+        const su2double f_e2 = 1.0 - max(f_t, f_l);
+        const su2double f_e1 = alpha >= 0.0 ? 2.0 * exp(-11.09*alpha*alpha) : 2.0 * exp(-9.0*alpha*alpha);
+        const su2double f_e = f_e2 * max((f_e1 - 1.0), 0.0);
+          
+
+        const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max);
+        const su2double l_LES = C_DES * Delta;
+
+        DES_lengthScale = ftilda_d *(1.0+f_e)*l_RANS + (1.0 - ftilda_d) * l_LES;
+
+        break;
+      }
+      case SST_SIDDES: {
+        
+        // Constants
+        const su2double C_w = 0.15;
+        const su2double C_dt1 = 20.0;
+        const su2double C_dt2 = 3.0;
+
+        const su2double alpha = 0.25 - sqrt(wallDist2) / h_max;
+        const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0);
+        const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2));
+        const su2double ftilda_d = max(1.0 - f_dt, f_b);
+
+        const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max);
+        const su2double l_LES = C_DES * Delta;
+
+        DES_lengthScale = ftilda_d*l_RANS + (1.0 - ftilda_d) * l_LES;
+
+        break;
+      }
+      case SST_EDDES: {
+
+        // Improved DDES version with the Shear-Layer-Adapted augmentation 
+        // found in Detached Eddy Simulation: Recent Development and Application to Compressor Tip Leakage Flow, Xiao He, Fanzhou Zhao, Mehdi Vahdati
+        // originally from Application of SST-Based SLA-DDES Formulation to Turbomachinery Flows, Guoping Xia, Zifei Yin and Gorazd Medic
+        // I could be naming it either as SST_EDDES to follow the same notation as for the SA model or as SST_SLA_DDES to follow the paper notation
+
+        const su2double f_max = 1.0, f_min = 0.1, a1 = 0.15, a2 = 0.3;
+
+        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);
+
+        const su2double omega = GeometryToolbox::Norm(3, Vorticity);
+
+        su2double ratioOmega[MAXNDIM] = {};
+
+        for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
+          ratioOmega[iDim] = Vorticity[iDim]/omega;
+        }
+
+        const su2double deltaDDES = geometry->nodes->GetMaxLength(iPoint);
+
+        su2double ln_max = 0.0;
+        for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+          const auto coord_j = geometry->nodes->GetCoord(jPoint);
+          su2double delta[MAXNDIM] = {};
+          for (auto iDim = 0u; iDim < nDim; iDim++){
+            delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);
+          }
+          su2double ln[3];
+          ln[0] = delta[1]*ratioOmega[2] - delta[2]*ratioOmega[1];
+          ln[1] = delta[2]*ratioOmega[0] - delta[0]*ratioOmega[2];
+          ln[2] = delta[0]*ratioOmega[1] - delta[1]*ratioOmega[0];
+          const su2double aux_ln = sqrt(ln[0]*ln[0] + ln[1]*ln[1] + ln[2]*ln[2]);
+          ln_max = max(ln_max, aux_ln);
+          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
+        }
+        vortexTiltingMeasure /= (nNeigh + 1);
+
+        
+
+        const su2double f_kh = max(f_min,
+                                   min(f_max,
+                                       f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
+
+        const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double C_d1 = 20.0;
+        const su2double C_d2 = 3.0;
+
+        const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2));
+
+        su2double delta = (ln_max/sqrt(3.0)) * f_kh;
+        if (f_d < 0.99){
+          delta = h_max;
+        }
+
+        const su2double l_LES = C_DES * delta;
+        DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES);
+      }
+      case SST_EDDES_UNSTR: {
+
+        // Improved DDES version with the Shear-Layer-Adapted augmentation 
+        // found in Detached Eddy Simulation: Recent Development and Application to Compressor Tip Leakage Flow, Xiao He, Fanzhou Zhao, Mehdi Vahdati
+        // originally from Application of SST-Based SLA-DDES Formulation to Turbomachinery Flows, Guoping Xia, Zifei Yin and Gorazd Medic
+        // I could be naming it either as SST_EDDES to follow the same notation as for the SA model or as SST_SLA_DDES to follow the paper notation
+
+        const su2double f_max = 1.0, f_min = 0.1, a1 = 0.15, a2 = 0.3;
+        
+        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);
+
+        su2double deltaOmega = -1.0;
+        su2double vorticityDir[MAXNDIM] = {};
+
+        for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
+          vorticityDir[iDim] = Vorticity[iDim]/VortMag;
+        }
+        
+        for (const auto jPoint : geometry->nodes->GetPoints(iPoint)){
+          const auto coord_j = geometry->nodes->GetCoord(jPoint);
+
+          for (const auto kPoint : geometry->nodes->GetPoints(iPoint)){
+            const auto coord_k = geometry->nodes->GetCoord(kPoint);
+
+            su2double delta[MAXNDIM] = {};
+            for (auto iDim = 0u; iDim < MAXNDIM; iDim++){
+              // TODO: Should I divide by 2 as I am interested in the dual volume (the edge is split at midpoint)?
+              delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0; 
+            }
+            su2double l_n_minus_m[MAXNDIM];
+            GeometryToolbox::CrossProduct(delta, vorticityDir, l_n_minus_m);
+            deltaOmega = max(deltaOmega, GeometryToolbox::Norm(3, l_n_minus_m));
+          }
+
+          // Add to VTM(iPoint) to perform the average
+          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
+        }
+        deltaOmega /= sqrt(3.0);
+        vortexTiltingMeasure /= (nNeigh+1);
+
+        const su2double f_kh = max(f_min,
+                                   min(f_max,
+                                       f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
+
+        const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double C_d1 = 20.0;
+        const su2double C_d2 = 3.0;
+
+        const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2));
+
+        su2double delta = deltaOmega * f_kh;
+        if (f_d < 0.99){
+          delta = h_max;
+        }
+
+        const su2double l_LES = C_DES * delta;
+        DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES);
+      }
+    }
+
+    nodes->SetDES_LengthScale(iPoint, DES_lengthScale);
+
+  }
+  END_SU2_OMP_FOR
+}
+
 void CTurbSSTSolver::ComputeUnderRelaxationFactor(const CConfig *config) {
 
   const su2double allowableRatio = config->GetMaxUpdateFractionSST();
diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp
index 75119c89ed5f..a9227af91400 100644
--- a/SU2_CFD/src/variables/CTurbSAVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp
@@ -45,54 +45,4 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi
     Solution_time_n  = Solution;
     Solution_time_n1 = Solution;
   }
-
-  DES_LengthScale.resize(nPoint) = su2double(0.0);
-  Vortex_Tilting.resize(nPoint);
-}
-
-void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
-                                        const su2double* Vorticity, su2double LaminarViscosity) {
-
-  su2double Strain[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}, Omega, StrainDotVort[3], numVecVort[3];
-  su2double numerator, trace0, trace1, denominator;
-
-  AD::StartPreacc();
-  AD::SetPreaccIn(PrimGrad_Flow, nDim+1, nDim);
-  AD::SetPreaccIn(Vorticity, 3);
-  /*--- Eddy viscosity ---*/
-  AD::SetPreaccIn(muT(iPoint));
-  /*--- Laminar viscosity --- */
-  AD::SetPreaccIn(LaminarViscosity);
-
-  Strain[0][0] = PrimGrad_Flow[1][0];
-  Strain[1][0] = 0.5*(PrimGrad_Flow[2][0] + PrimGrad_Flow[1][1]);
-  Strain[0][1] = 0.5*(PrimGrad_Flow[1][1] + PrimGrad_Flow[2][0]);
-  Strain[1][1] = PrimGrad_Flow[2][1];
-  if (nDim == 3){
-    Strain[0][2] = 0.5*(PrimGrad_Flow[3][0] + PrimGrad_Flow[1][2]);
-    Strain[1][2] = 0.5*(PrimGrad_Flow[3][1] + PrimGrad_Flow[2][2]);
-    Strain[2][0] = 0.5*(PrimGrad_Flow[1][2] + PrimGrad_Flow[3][0]);
-    Strain[2][1] = 0.5*(PrimGrad_Flow[2][2] + PrimGrad_Flow[3][1]);
-    Strain[2][2] = PrimGrad_Flow[3][2];
-  }
-
-  Omega = sqrt(Vorticity[0]*Vorticity[0] + Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]);
-
-  StrainDotVort[0] = Strain[0][0]*Vorticity[0]+Strain[0][1]*Vorticity[1]+Strain[0][2]*Vorticity[2];
-  StrainDotVort[1] = Strain[1][0]*Vorticity[0]+Strain[1][1]*Vorticity[1]+Strain[1][2]*Vorticity[2];
-  StrainDotVort[2] = Strain[2][0]*Vorticity[0]+Strain[2][1]*Vorticity[1]+Strain[2][2]*Vorticity[2];
-
-  numVecVort[0] = StrainDotVort[1]*Vorticity[2] - StrainDotVort[2]*Vorticity[1];
-  numVecVort[1] = StrainDotVort[2]*Vorticity[0] - StrainDotVort[0]*Vorticity[2];
-  numVecVort[2] = StrainDotVort[0]*Vorticity[1] - StrainDotVort[1]*Vorticity[0];
-
-  numerator = sqrt(6.0) * sqrt(numVecVort[0]*numVecVort[0] + numVecVort[1]*numVecVort[1] + numVecVort[2]*numVecVort[2]);
-  trace0 = 3.0*(pow(Strain[0][0],2.0) + pow(Strain[1][1],2.0) + pow(Strain[2][2],2.0));
-  trace1 = pow(Strain[0][0] + Strain[1][1] + Strain[2][2],2.0);
-  denominator = pow(Omega, 2.0) * sqrt(trace0-trace1);
-
-  Vortex_Tilting(iPoint) = (numerator/denominator) * max(1.0,0.2*LaminarViscosity/muT(iPoint));
-
-  AD::SetPreaccOut(Vortex_Tilting(iPoint));
-  AD::EndPreacc();
-}
+}
\ No newline at end of file
diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp
index 47ca7d3bd01b..171ec8a7ec35 100644
--- a/SU2_CFD/src/variables/CTurbVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbVariable.cpp
@@ -27,6 +27,7 @@
 
 
 #include "../../include/variables/CTurbVariable.hpp"
+#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
 
 CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config)
@@ -35,8 +36,46 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
     turb_index.resize(nPoint) = su2double(1.0);
     intermittency.resize(nPoint) = su2double(1.0);
 
+    DES_LengthScale.resize(nPoint) = su2double(0.0);
+    
+    Vortex_Tilting.resize(nPoint);
+
    }
 
+
+void CTurbVariable::SetVortex_Tilting(unsigned long iPoint, su2double **Strain,
+                                      const su2double* Vorticity, su2double LaminarViscosity) {
+
+  su2double Omega, StrainDotVort[3], numVecVort[3];
+  su2double numerator, trace0, trace1, denominator;
+
+  AD::StartPreacc();
+  AD::SetPreaccIn(Strain, nDim, nDim);
+  AD::SetPreaccIn(Vorticity, 3);
+  /*--- Eddy viscosity ---*/
+  AD::SetPreaccIn(muT(iPoint));
+  /*--- Laminar viscosity --- */
+  AD::SetPreaccIn(LaminarViscosity);
+
+  Omega = GeometryToolbox::Norm(3, Vorticity);
+
+  StrainDotVort[0] = GeometryToolbox::DotProduct(3, Strain[0], Vorticity);
+  StrainDotVort[1] = GeometryToolbox::DotProduct(3, Strain[1], Vorticity);
+  StrainDotVort[2] = GeometryToolbox::DotProduct(3, Strain[2], Vorticity);
+
+  GeometryToolbox::CrossProduct(StrainDotVort, Vorticity, numVecVort);
+
+  numerator = sqrt(6.0) * GeometryToolbox::Norm(3, numVecVort);
+  trace0 = 3.0*(pow(Strain[0][0],2.0) + pow(Strain[1][1],2.0) + pow(Strain[2][2],2.0));
+  trace1 = pow(Strain[0][0] + Strain[1][1] + Strain[2][2],2.0);
+  denominator = pow(Omega, 2.0) * sqrt(trace0-trace1);
+
+  Vortex_Tilting(iPoint) = (numerator/denominator) * max(1.0,0.2*LaminarViscosity/muT(iPoint));
+
+  AD::SetPreaccOut(Vortex_Tilting(iPoint));
+  AD::EndPreacc();
+}
+
 void CTurbVariable::RegisterEddyViscosity(bool input) {
   RegisterContainer(input, muT);
-}
+}
\ No newline at end of file
diff --git a/externals/ninja b/externals/ninja
index b4d51f6ed5be..52649de2c56b 160000
--- a/externals/ninja
+++ b/externals/ninja
@@ -1 +1 @@
-Subproject commit b4d51f6ed5bed09dd2b70324df0d9cb4ecad2638
+Subproject commit 52649de2c56b63f42bc59513d51286531c595b44