diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index ddf570ab697..8fca4a72895 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -994,6 +994,8 @@ enum class SST_OPTIONS {
   COMP_Wilcox, /*!< \brief Menter k-w SST model with Compressibility correction of Wilcox. */
   COMP_Sarkar, /*!< \brief Menter k-w SST model with Compressibility correction of Sarkar. */
   DLL,         /*!< \brief Menter k-w SST model with dimensionless lower limit clipping of turbulence variables. */
+  SAS_TRAVIS,         /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */
+  SAS_BABU,         /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */
 };
 static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("NONE", SST_OPTIONS::NONE)
@@ -1009,6 +1011,8 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("COMPRESSIBILITY-WILCOX", SST_OPTIONS::COMP_Wilcox)
   MakePair("COMPRESSIBILITY-SARKAR", SST_OPTIONS::COMP_Sarkar)
   MakePair("DIMENSIONLESS_LIMIT", SST_OPTIONS::DLL)
+  MakePair("SAS_TRAVIS", SST_OPTIONS::SAS_TRAVIS)
+  MakePair("SAS_BABU", SST_OPTIONS::SAS_BABU)
 };
 
 /*!
@@ -1017,6 +1021,7 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
 struct SST_ParsedOptions {
   SST_OPTIONS version = SST_OPTIONS::V1994;   /*!< \brief Enum SST base model. */
   SST_OPTIONS production = SST_OPTIONS::NONE; /*!< \brief Enum for production corrections/modifiers for SST model. */
+  SST_OPTIONS sasModel = SST_OPTIONS::NONE;   /*!< \brief Enum SST base model. */
   bool sust = false;                          /*!< \brief Bool for SST model with sustaining terms. */
   bool uq = false;                            /*!< \brief Bool for using uncertainty quantification. */
   bool modified = false;                      /*!< \brief Bool for modified (m) SST model. */
@@ -1061,6 +1066,15 @@ 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);
+  const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_TRAVIS);
+  const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_BABU);
+  if (sst_sas_simple && sst_sas_comp) {
+    SU2_MPI::Error("Two versions (Simple and Complicated) selected for SAS under SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
+  } else if (sst_sas_simple) {
+    SSTParsedOptions.sasModel = SST_OPTIONS::SAS_TRAVIS;
+  } else {
+    SSTParsedOptions.sasModel = SST_OPTIONS::SAS_BABU;
+  } 
 
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
@@ -1433,7 +1447,10 @@ enum ENUM_HYBRIDRANSLES {
   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 = 4,           /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */
+  SST_DDES = 5,           /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): DDES). */
+  SST_IDDES = 6,           /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */
+  SST_SIDDES = 7           /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */
 };
 static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("NONE", NO_HYBRIDRANSLES)
@@ -1441,6 +1458,9 @@ 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("SST_DDES", SST_DDES)
+  MakePair("SST_IDDES", SST_IDDES)
+  MakePair("SST_SIDDES", SST_SIDDES)
 };
 
 /*!
@@ -2535,6 +2555,7 @@ enum PERIODIC_QUANTITIES {
   PERIODIC_LIM_PRIM_1 ,  /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
   PERIODIC_LIM_PRIM_2 ,  /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
   PERIODIC_IMPLICIT   ,  /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
+  PERIODIC_VEL_LAPLACIAN  ,  /*!< \brief Velocity Laplacian communication for SAS (periodic only). */
 };
 
 /*!
@@ -2566,6 +2587,7 @@ enum class MPI_QUANTITIES {
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
+  VELOCITY_LAPLACIAN   ,  /*!< \brief Velocity Laplacian communication. */
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 8a60579189c..aa6ea868d22 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -6275,6 +6275,17 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
             else cout << "\nusing default hard coded lower limit clipping";
 
             cout << "." << endl;
+            if(sstParsedOptions.sasModel != SST_OPTIONS::NONE) cout << "Scale Adaptive Simulation model: ";
+            switch (sstParsedOptions.sasModel) {
+              case SST_OPTIONS::SAS_TRAVIS:
+                cout << "Travis et al. (2004)";
+                break;
+              case SST_OPTIONS::SAS_BABU:
+                cout << "Babu et al. (2016)";
+                break;
+              default:
+                break;
+            }
             break;
         }
         switch (Kind_Trans_Model) {
@@ -6316,6 +6327,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           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 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 535335947ca..c6297f6427e 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -136,6 +136,9 @@ 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 */
+  su2double FTrans;   /*!< \brief SAS function */
+  su2double *VelLapl; /*!< \brief Laplacian of the velocity */
   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. */
@@ -707,6 +710,18 @@ class CNumerics {
    */
   virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */};
 
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline virtual su2double GetFTrans() const { return 0.0; }
+
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline void SetVelLapl(su2double* val_VelLapl) {
+    VelLapl = val_VelLapl;
+  }
+
   /*!
    * \brief Set the value of the effective intermittency for the LM model.
    * \param[in] intermittency_eff_i - Value of the effective intermittency at point i.
@@ -828,6 +843,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 c9464b1dfc3..8b2bca6bcbd 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -652,6 +652,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   /*--- Closure constants ---*/
   const su2double sigma_k_1, sigma_k_2, sigma_w_1, sigma_w_2, beta_1, beta_2, beta_star, a1, alfa_1, alfa_2;
   const su2double prod_lim_const;
+  const su2double cTrans;
 
   /*--- Ambient values for SST-SUST. ---*/
   const su2double kAmb, omegaAmb;
@@ -736,6 +737,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
         alfa_1(constants[8]),
         alfa_2(constants[9]),
         prod_lim_const(constants[10]),
+        cTrans(1.25),
         kAmb(val_kine_Inf),
         omegaAmb(val_omega_Inf) {
     /*--- "Allocate" the Jacobian using the static buffer. ---*/
@@ -903,6 +905,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
+      if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES)
+        dk = Density_i * sqrt(ScalarVar_i[0]*ScalarVar_i[0]*ScalarVar_i[0]) / lengthScale_i;
+        
       su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
@@ -916,9 +921,43 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       Residual[0] += pk * Volume;
       Residual[1] += pw * Volume;
 
-      /*--- Add the dissipation  terms to the residuals.---*/
+      /*--- Compute SAS source term. ---*/
+      su2double Q_SAS = 0.0;
+      
+      if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) {
+                
+        const su2double KolmConst = 0.41;
+        const su2double csi2 = 3.51;
+        const su2double sigma_phi = 2.0/3.0;
+        const su2double C = 2.0;
+        const su2double C_s = 0.5; // Honestly I do not know if it is the right one.
+        const su2double gridSize = pow(Volume, 1.0/3.0);
+        // Scale of the modeled turbulence
+        const su2double L = sqrt(ScalarVar_i[0]) / (pow(beta_star, 0.25) * ScalarVar_i[1]);
+        // Von Karman Length Scale
+        const su2double VelLaplMag = GeometryToolbox::SquaredNorm(nDim, VelLapl);
+        const su2double L_vK_1 = KolmConst * StrainMag_i / sqrt(VelLaplMag);
+        const su2double L_vK_2 = C_s * sqrt(KolmConst * csi2 / (beta_blended/beta_star - alfa_blended)) * gridSize;
+        const su2double L_vK = max(L_vK_1, L_vK_2);
+        
+        const su2double gradTKE = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0])/(ScalarVar_i[0]*ScalarVar_i[0]);
+        const su2double gradOmega = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[1])/(ScalarVar_i[1]*ScalarVar_i[1]);
+
+        const su2double Q_SAS_1 = csi2 * KolmConst * StrainMag_i * StrainMag_i * (L/L_vK) * (L/L_vK);
+        const su2double Q_SAS_2 = C * (2*ScalarVar_i[0] / sigma_phi) * max(gradOmega, gradTKE);
+        
+        Q_SAS = max(Q_SAS_1 - Q_SAS_2, 0.0);
+
+        Residual[1] += Density_i * Q_SAS * Volume;
+
+      }
 
-      Residual[0] -= dk * Volume;
+      /*--- Add the dissipation  terms to the residuals.---*/
+      FTrans = 1.0;
+      if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) {
+        FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0));
+      }
+      Residual[0] -= dk * Volume * FTrans;
       Residual[1] -= dw * Volume;
 
       /*--- Cross diffusion ---*/
@@ -931,10 +970,14 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       /*--- Implicit part ---*/
 
-      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt);
-      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt);
+      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt) * FTrans;
+      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt) * FTrans;
       Jacobian_i[1][0] = 0.0;
       Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt);
+
+      if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) {
+        Jacobian_i[0][0] += Q_SAS * Volume / ScalarVar_i[0];
+      }
     }
 
     AD::SetPreaccOut(Residual, nVar);
@@ -942,4 +985,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
     return ResidualType<>(Residual, Jacobian_i, nullptr);
   }
+
+  /*!
+   * \brief Get the value of the FTrans.
+   */
+  inline su2double GetFTrans() const override { return FTrans; }
+
+
 };
diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
index 78aa7720095..9df3e0b0582 100644
--- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
+++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp
@@ -54,6 +54,16 @@ class CTurbSSTSolver final : public CTurbSolver {
                      const CConfig *config,
                      unsigned short val_marker);
 
+   /*!
+   * \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 f0e5ac2957a..5919face056 100644
--- a/SU2_CFD/include/variables/CTurbSAVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp
@@ -39,9 +39,10 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
   VectorType Vortex_Tilting;
 
+  VectorType k, Omega;  /*!< \brief SST variables as computed through SA solution. */
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -60,19 +61,6 @@ 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.
diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp
index dcc4772100e..42f938566d5 100644
--- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp
@@ -45,6 +45,16 @@ class CTurbSSTVariable final : public CTurbVariable {
   VectorType F2;    /*!< \brief Menter blending function for blending of k-w and k-eps. */
   VectorType CDkw;  /*!< \brief Cross-diffusion. */
   SST_ParsedOptions sstParsedOptions;
+
+  VectorType ftilda_d;
+  VectorType l_RANS;
+  VectorType l_LES;
+  VectorType r_dl;
+  VectorType r_dt;
+  VectorType r_d;
+  VectorType FTrans;
+  MatrixType VelocityLaplacian;
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -87,4 +97,99 @@ class CTurbSSTVariable final : public CTurbVariable {
    * \brief Get the value of the cross diffusion of tke and omega.
    */
   inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); }
+
+  /*!
+   * \brief Set the DES Length Scale.
+   * \param[in] iPoint - Point index.
+   */
+  inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) override { 
+    ftilda_d(iPoint) = val_ftilda_d;
+    l_RANS(iPoint) = val_l_RANS;
+    l_LES(iPoint) = val_l_LES;
+
+    if ( config->GetKind_HybridRANSLES() == SST_DDES)
+      r_d(iPoint) = val_r_d;
+    else
+      r_dt(iPoint) = val_r_d;
+  }
+
+  inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) override { 
+    ftilda_d(iPoint) = val_ftilda_d;
+    l_RANS(iPoint) = val_l_RANS;
+    l_LES(iPoint) = val_l_LES;
+    r_dt(iPoint) = val_r_dt;
+    r_dl(iPoint) = val_r_dl;
+  }
+
+  inline su2double Get_L_RANS(unsigned long iPoint) const override { return l_RANS(iPoint); }
+  inline su2double Get_L_LES(unsigned long iPoint) const override { return l_LES(iPoint); }
+  inline su2double Get_ftilda_d(unsigned long iPoint) const override { return ftilda_d(iPoint); }
+  inline su2double Get_r_dt(unsigned long iPoint) const override { return r_dt(iPoint); }
+  inline su2double Get_r_dl(unsigned long iPoint) const override { return r_dl(iPoint); }
+  inline su2double Get_r_d(unsigned long iPoint) const override { return r_d(iPoint); }
+
+
+  /*!
+   * \brief Get the value of the FTrans.
+   * \param[in] iPoint - Point index.
+   */
+  inline su2double GetFTrans(unsigned long iPoint) const override { return FTrans(iPoint); }
+
+  /*!
+   * \brief Set the value of the FTrans.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_FTrans - Value of the FTrans variable.
+   */
+  inline void SetFTrans(unsigned long iPoint, su2double val_FTrans) override { 
+    FTrans(iPoint) = val_FTrans;
+  };
+
+  /*!
+   * \brief Get the value of the velocity laplacian.
+   * \param[in] iPoint - Point index.
+   */
+  inline su2double* GetVelLapl(unsigned long iPoint) final { return VelocityLaplacian[iPoint]; }
+
+  /*!
+   * \brief Get the value of the velocity laplacian.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   */
+  inline su2double GetVelLapl(unsigned long iPoint, unsigned short iDim) const final { return VelocityLaplacian[iPoint][iDim]; }
+
+  /*!
+   * \brief Incrementally add the velocity laplacian vector.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_VelLapl_X - X-Component of the velocity laplacian.
+   * \param[in] val_VelLapl_Y - Y-Component of the velocity laplacian.
+   * \param[in] val_VelLapl_Z - Z-Component of the velocity laplacian.
+   */
+  inline void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y, su2double val_VelLapl_Z) override { 
+    VelocityLaplacian(iPoint, 0) += val_VelLapl_X;
+    VelocityLaplacian(iPoint, 1) += val_VelLapl_Y;
+    if(nDim == 3) VelocityLaplacian(iPoint, 2) += val_VelLapl_Z;
+  };
+
+  /*!
+   * \brief Set the value of the velocity laplacian.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_VelLapl - Vector of velocity laplacian.
+   */
+  inline void SetVelLapl(unsigned long iPoint, su2double* val_VelLapl) override { 
+    VelocityLaplacian(iPoint, 0) = val_VelLapl[0];
+    VelocityLaplacian(iPoint, 1) = val_VelLapl[1];
+    if(nDim == 3) VelocityLaplacian(iPoint, 2) = val_VelLapl[2];
+  };
+
+  /*!
+   * \brief Set the value of the velocity laplacian.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_VelLapl - value of the velocity laplacian.
+   */
+  inline void SetVelLapl(unsigned long iPoint, unsigned short iDim, su2double val_VelLapl) override { 
+    VelocityLaplacian(iPoint, iDim) = val_VelLapl;
+  };
+
+
 };
diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp
index 48f2087b2bf..794c8b03065 100644
--- a/SU2_CFD/include/variables/CTurbVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbVariable.hpp
@@ -38,11 +38,13 @@
 class CTurbVariable : public CScalarVariable {
 protected:
   VectorType muT; /*!< \brief Eddy viscosity. */
+  VectorType DES_LengthScale;
 
 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 SRSGridSize;           /*!< \brief alue of the desired grid size for Scale Resolving Simulations. */
 
   /*!
    * \brief Constructor of the class.
@@ -99,6 +101,33 @@ class CTurbVariable : public CScalarVariable {
    * \param[in] val_intermittency - New value of the intermittency.
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
+  
+  /*!
+   * \brief Get the desired grid size for Scale Resolving Simulations.
+   * \param[in] iPoint - Point index.
+   * \return the value of the desired grid size for Scale Resolving Simulations.
+   */
+  inline su2double GetSRSGridSize(unsigned long iPoint) const final { return SRSGridSize(iPoint); }
+
+  /*!
+   * \brief Set the value of the desired grid size for Scale Resolving Simulations.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_muT - Value of the desired grid size for Scale Resolving Simulations.
+   */
+  inline void SetSRSGridSize(unsigned long iPoint, su2double val_gridSize) final { SRSGridSize(iPoint) = val_gridSize; }
+
+  /*!
+   * \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; }
 
 };
 
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index bc855bc5dca..d01c9363326 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -381,6 +381,25 @@ class CVariable {
    */
   inline virtual void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) {}
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) {}
+  
+  inline virtual su2double Get_L_RANS(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_L_LES(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_ftilda_d(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_dt(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_dl(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_d(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) {}
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.
@@ -1657,6 +1676,41 @@ class CVariable {
    */
   inline virtual su2double GetCrossDiff(unsigned long iPoint) const { return 0.0; }
 
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline virtual su2double GetFTrans(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief Set the value of the value of FTrans.
+   */
+  inline virtual void SetFTrans(unsigned long iPoint, su2double val_FTrans) {}
+
+  /*!
+   * \brief Get the value of the velocity laplacian.
+   */
+  inline virtual su2double* GetVelLapl(unsigned long iPoint) { return nullptr; }
+
+  /*!
+   * \brief Get the value of the velocity laplacian.
+   */
+  inline virtual su2double GetVelLapl(unsigned long iPoint, unsigned short iDim) const { return 0.0; }
+
+  /*!
+   * \brief Incrementally add the velocity laplacian vector.
+   */
+  inline virtual void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y, su2double val_VelLapl_Z) {}
+  
+  /*!
+   * \brief Set the value of the velocity laplacian.
+   */
+  inline virtual void SetVelLapl(unsigned long iPoint, su2double* val_VelLapl) {}
+
+  /*!
+   * \brief Set the value of the velocity laplacian.
+   */
+  inline virtual void SetVelLapl(unsigned long iPoint, unsigned short iDim, su2double val_VelLapl) {}
+
   /*!
    * \brief Get the value of the eddy viscosity.
    * \return the value of the eddy viscosity.
@@ -1705,6 +1759,18 @@ class CVariable {
    */
   inline virtual void SetmuT(unsigned long iPoint, su2double val_muT) {}
 
+  /*!
+   * \brief Get the value of the desired grid size for Scale Resolving Simulations.
+    * \return the value of the desired grid size for Scale Resolving Simulations.
+  */
+  inline virtual su2double GetSRSGridSize(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief Set the value of the desired grid size for Scale Resolving Simulations.
+   * \param[in] val_gridSize
+   */
+  inline virtual void SetSRSGridSize(unsigned long iPoint, su2double val_gridSize) {}
+
   /*!
    * \brief Set the value of the turbulence index.
    * \param[in] val_turb_index - turbulence index
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 3c651269b02..40ea7e0fd23 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1247,11 +1247,19 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){
   switch (TurbModelFamily(config->GetKind_Turb_Model())) {
     case TURB_FAMILY::SA:
       AddVolumeOutput("NU_TILDE", "Nu_Tilde", "SOLUTION", "Spalart-Allmaras variable");
+      AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations");
       break;
 
     case TURB_FAMILY::KW:
       AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy");
       AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation");
+      AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations");
+      if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_TRAVIS) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation");
+      if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_BABU){
+        AddVolumeOutput("VEL-LAPLACIAN_X", "Vel Laplacian x", "SOLUTION", "value of laplacian of x-velocity for SAS simulation");
+        AddVolumeOutput("VEL-LAPLACIAN_Y", "Vel Laplacian y", "SOLUTION", "value of laplacian of y-velocity for SAS simulation");
+        if (nDim == 3) AddVolumeOutput("VEL-LAPLACIAN_Z", "Vel Laplacian z", "SOLUTION", "value of laplacian of z-velocity for SAS simulation");
+      }
       break;
 
     case TURB_FAMILY::NONE:
@@ -1468,6 +1476,20 @@ 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");
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST) {
+      AddVolumeOutput("F_D", "f_d", "DDES", "Empiric blending function");
+      AddVolumeOutput("L_RANS", "l_RANS", "DDES", "RANS length scale value");
+      AddVolumeOutput("L_LES", "l_LES", "DDES", "LES length scale value");
+    }
+    if ( config->GetKind_HybridRANSLES() == SST_DDES){
+      AddVolumeOutput("R_D", "r_d", "DDES", "r_d");
+    } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){
+      AddVolumeOutput("R_DT", "r_dt", "DDES", "turbulent r_d");
+      AddVolumeOutput("R_DL", "r_dl", "DDES", "laminar r_d");
+    } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){
+      AddVolumeOutput("R_DT", "r_dt", "DDES", "turbulent r_d");
+    }
+    AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations");
   }
 
   if (config->GetViscous()) {
@@ -1520,6 +1542,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     case TURB_FAMILY::SA:
       SetVolumeOutputValue("NU_TILDE", iPoint, Node_Turb->GetSolution(iPoint, 0));
       SetVolumeOutputValue("RES_NU_TILDE", iPoint, turb_solver->LinSysRes(iPoint, 0));
+      SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint));
       if (limiter) {
         SetVolumeOutputValue("LIMITER_NU_TILDE", iPoint, Node_Turb->GetLimiter(iPoint, 0));
       }
@@ -1530,6 +1553,13 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
       SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1));
       SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0));
       SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1));
+      SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, Node_Turb->GetSRSGridSize(iPoint));
+      if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_TRAVIS) SetVolumeOutputValue("FTRANS", iPoint, Node_Turb->GetFTrans(iPoint));
+      if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_BABU){
+        SetVolumeOutputValue("VEL-LAPLACIAN_X", iPoint, Node_Turb->GetVelLapl(iPoint, 0));
+        SetVolumeOutputValue("VEL-LAPLACIAN_Y", iPoint, Node_Turb->GetVelLapl(iPoint, 1));
+        if (nDim == 3) SetVolumeOutputValue("VEL-LAPLACIAN_Z", iPoint, Node_Turb->GetVelLapl(iPoint, 2));
+      } 
       if (limiter) {
         SetVolumeOutputValue("LIMITER_TKE", iPoint, Node_Turb->GetLimiter(iPoint, 0));
         SetVolumeOutputValue("LIMITER_DISSIPATION", iPoint, Node_Turb->GetLimiter(iPoint, 1));
@@ -1567,6 +1597,27 @@ 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));
+    if ( config->GetKind_HybridRANSLES() == SST_DDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_D", iPoint, Node_Turb->Get_r_d(iPoint));
+    } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint));
+      SetVolumeOutputValue("R_DL", iPoint, Node_Turb->Get_r_dl(iPoint));
+    } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(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);
   }
 
   switch (config->GetKind_Species_Model()) {
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index e446c3f1133..d19df52ab70 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -275,6 +275,11 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config,
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       ICOUNT           = nVar;
       break;
+    case PERIODIC_VEL_LAPLACIAN:
+      COUNT_PER_POINT  = nDim;
+      MPI_TYPE         = COMM_TYPE_DOUBLE;
+      ICOUNT           = nDim;
+      break;
     default:
       SU2_MPI::Error("Unrecognized quantity for periodic communication.",
                      CURRENT_FUNCTION);
@@ -363,6 +368,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
   auto *Diff      = new su2double[nVar];
   auto *Und_Lapl  = new su2double[nVar];
+  auto *Vel_Lapl  = new su2double[nDim];
   auto *Sol_Min   = new su2double[nPrimVarGrad];
   auto *Sol_Max   = new su2double[nPrimVarGrad];
   auto *rotPrim_i = new su2double[nPrimVar];
@@ -968,6 +974,58 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             }
 
             break;
+          case PERIODIC_VEL_LAPLACIAN:
+
+            /*--- For JST, the undivided Laplacian must be computed
+             consistently by using the complete control volume info
+             from both sides of the periodic face. ---*/
+
+            for (iDim = 0; iDim < nDim; iDim++)
+              Vel_Lapl[iDim] = 0.0;
+
+            for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+
+              /*--- Avoid periodic boundary points so that we do not
+               duplicate edges on both sides of the periodic BC. ---*/
+
+              if (!geometry->nodes->GetPeriodicBoundary(jPoint)) {
+
+                /*--- Solution differences ---*/
+                const su2double distance = GeometryToolbox::Distance(nDim, geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint));
+
+                for (iDim = 0; iDim < nDim; iDim++)
+                Diff[iDim] = (base_nodes->GetVelocity(iPoint, iDim) -
+                              base_nodes->GetVelocity(jPoint, iDim))/distance;
+
+
+                boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint);
+                boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);
+
+                /*--- Both points inside the domain, or both in the boundary ---*/
+                /*--- iPoint inside the domain, jPoint on the boundary ---*/
+
+                if (!boundary_i || boundary_j) {
+                  if (geometry->nodes->GetDomain(iPoint)){
+                    for (iDim = 0; iDim< nDim; iDim++)
+                    Vel_Lapl[iDim] -= Diff[iDim];
+                  }
+                }
+              }
+            }
+
+            /*--- Store the components to be communicated in the buffer. ---*/
+
+            for (iDim = 0; iDim < nDim; iDim++)
+              bufDSend[buf_offset+iDim] = Vel_Lapl[iDim];
+
+            /*--- Rotate the momentum components of the Laplacian. ---*/
+            
+            if (rotate_periodic) {
+              Rotate(zeros, &Vel_Lapl[0], &bufDSend[buf_offset+0]);
+            }
+
+            break;
+
 
           default:
             SU2_MPI::Error("Unrecognized quantity for periodic communication.",
@@ -1290,7 +1348,16 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
                 limiter(iPoint, iVar) = min(limiter(iPoint, iVar), bufDRecv[buf_offset+iVar]);
 
               break;
+            case PERIODIC_VEL_LAPLACIAN:
+            {
+              /*--- Adjust the undivided Laplacian. The accumulation was
+               with a subtraction before communicating, so now just add. ---*/
+              su2double bufDRecv_Z = 0.0;
+              if(nDim == 3) bufDRecv_Z = bufDRecv[buf_offset+2];
+              base_nodes->AddVelLapl(iPoint, bufDRecv[buf_offset+0], bufDRecv[buf_offset+1], bufDRecv_Z);
 
+              break;
+            }
             default:
 
               SU2_MPI::Error("Unrecognized quantity for periodic communication.",
@@ -1379,6 +1446,10 @@ void CSolver::GetCommCountAndType(const CConfig* config,
       COUNT_PER_POINT  = nVar;
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       break;
+    case MPI_QUANTITIES::VELOCITY_LAPLACIAN:
+      COUNT_PER_POINT  = nDim;
+      MPI_TYPE         = COMM_TYPE_DOUBLE;
+      break;
     default:
       SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.",
                      CURRENT_FUNCTION);
@@ -1528,6 +1599,10 @@ void CSolver::InitiateComms(CGeometry *geometry,
             for (iVar = 0; iVar < nVar; iVar++)
               bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar);
             break;
+          case MPI_QUANTITIES::VELOCITY_LAPLACIAN:
+            for(iDim = 0; iDim < nDim; iDim++)
+              bufDSend[buf_offset+iDim] = base_nodes->GetVelLapl(iPoint, iDim);
+            break;
           default:
             SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.",
                            CURRENT_FUNCTION);
@@ -1676,6 +1751,10 @@ void CSolver::CompleteComms(CGeometry *geometry,
             for (iVar = 0; iVar < nVar; iVar++)
               base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]);
             break;
+          case MPI_QUANTITIES::VELOCITY_LAPLACIAN:
+            for (iDim = 0; iDim < nDim; iDim++)
+              base_nodes->SetVelLapl(iPoint, iDim, bufDRecv[buf_offset+iDim]);
+            break;
           default:
             SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.",
                            CURRENT_FUNCTION);
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 6b465e802fd..141299f8940 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -224,7 +224,8 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain
     const su2double nu = mu/rho;
     const su2double nu_hat = nodes->GetSolution(iPoint,0);
     const su2double roughness = geometry->nodes->GetRoughnessHeight(iPoint);
-    const su2double dist = geometry->nodes->GetWall_Distance(iPoint) + rough_const * roughness;
+    const su2double wallDistance = geometry->nodes->GetWall_Distance(iPoint);
+    const su2double dist = wallDistance + rough_const * roughness;
 
     su2double Ji = nu_hat/nu;
     if (roughness > 1.0e-10)
diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
index 02b91c8ab2c..20f5f5602e2 100644
--- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
@@ -28,6 +28,7 @@
 #include "../../include/solvers/CTurbSSTSolver.hpp"
 #include "../../include/variables/CTurbSSTVariable.hpp"
 #include "../../include/variables/CFlowVariable.hpp"
+#include "../../include/variables/CMeshVariable.hpp"
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
@@ -197,8 +198,71 @@ 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) {
+
+    /*--- Compute the DES length scale ---*/
+
+    SetDES_LengthScale(solver_container, geometry, config);
+
+  }
+
+  if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU){
+    auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+    
+    SU2_OMP_FOR_DYN(256)
+    for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
+
+      const bool boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint);
+
+      /*--- Initialize. ---*/
+      for (unsigned short iDim = 0; iDim < nDim; iDim++)
+        nodes->SetVelLapl(iPoint, iDim, 0.0);
+
+      /*--- Loop over the neighbors of point i. ---*/
+      for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
+
+        bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);
+
+        /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/
+        if (boundary_i && !boundary_j) continue;
+
+        /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/
+        const su2double distance = GeometryToolbox::Distance(nDim, geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint));
+
+        const su2double delta_x = (flowNodes->GetVelocity(jPoint,0)-flowNodes->GetVelocity(iPoint,0))/distance;
+        const su2double delta_y = (flowNodes->GetVelocity(jPoint,1)-flowNodes->GetVelocity(iPoint,1))/distance;
+        su2double delta_z = 0.0;
+
+        if (nDim == 3) {
+          delta_z = (flowNodes->GetVelocity(jPoint,2)-flowNodes->GetVelocity(iPoint,2))/distance;
+        }
+        nodes->AddVelLapl(iPoint, delta_x, delta_y, delta_z);
+      
+      }
+
+      
+    }
+    END_SU2_OMP_FOR
+
+    /*--- Correct the Laplacian across any periodic boundaries. ---*/
+
+    for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) {
+      InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VEL_LAPLACIAN);
+      CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VEL_LAPLACIAN);
+    }
+
+    /*--- MPI parallelization ---*/
+
+    InitiateComms(geometry, config, MPI_QUANTITIES::VELOCITY_LAPLACIAN);
+    CompleteComms(geometry, config, MPI_QUANTITIES::VELOCITY_LAPLACIAN);
+
+  }
+
 }
 
 void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container,
@@ -245,6 +309,12 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
 
     nodes->SetmuT(iPoint, muT);
 
+    // Now compute desired cell size for Scale Resolving Simulations
+    const su2double RANSLength = sqrt(nodes->GetSolution(iPoint, 0)) / max(1e-20, (constants[6] * nodes->GetSolution(iPoint, 1)));
+    const su2double RatioL = 0.1;  // it should be less or equal than 0.2 - 0.1. Should be taken as input from config?
+    const su2double SRSGridSize = RANSLength * RatioL;
+    nodes->SetSRSGridSize(iPoint, SRSGridSize);
+
   }
   END_SU2_OMP_FOR
 
@@ -367,10 +437,23 @@ 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);
+    }
+
+    if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU){
+      numerics->SetVelLapl(nodes->GetVelLapl(iPoint));
+    }
+
     /*--- Compute the source term ---*/
 
     auto residual = numerics->ComputeResidual(config);
 
+    /*--- Store the SAS function ---*/
+    if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) {
+      nodes->SetFTrans(iPoint, numerics->GetFTrans());
+    }
+
     /*--- Store the intermittency ---*/
 
     if (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) {
@@ -1024,3 +1107,112 @@ 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 su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12);
+    const su2double VortMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 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){
+      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);
+
+        nodes->SetDebug_Quantities(config, iPoint, f_d, l_RANS, l_LES, r_d);
+
+        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;
+
+        nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dl, r_dt);
+
+        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 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 Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max);
+        const su2double l_LES = C_DES * Delta;
+
+        nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dt);
+
+        DES_lengthScale = ftilda_d*l_RANS + (1.0 - ftilda_d) * l_LES;
+
+        break;
+      }
+    }
+
+    nodes->SetDES_LengthScale(iPoint, DES_lengthScale);
+
+  }
+  END_SU2_OMP_FOR
+}
diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp
index d7091f8be59..d0466c83d6b 100644
--- a/SU2_CFD/src/variables/CTurbSAVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp
@@ -46,8 +46,10 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi
     Solution_time_n1 = Solution;
   }
 
-  DES_LengthScale.resize(nPoint) = su2double(0.0);
   Vortex_Tilting.resize(nPoint);
+
+  k.resize(nPoint) = su2double(0.0);
+  Omega.resize(nPoint) = su2double(0.0);
 }
 
 void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp
index d1fd98c8a14..2d24b5c4132 100644
--- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp
@@ -50,7 +50,19 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu
   F2.resize(nPoint) = su2double(0.0);
   CDkw.resize(nPoint) = su2double(0.0);
 
+  ftilda_d.resize(nPoint) = su2double(0.0);
+  l_RANS.resize(nPoint) = su2double(0.0);
+  l_LES.resize(nPoint) = su2double(0.0);
+  r_dl.resize(nPoint) = su2double(0.0);
+  r_dt.resize(nPoint) = su2double(0.0);
+  r_d.resize(nPoint) = su2double(0.0);
+
   muT.resize(nPoint) = mut;
+
+  if(sstParsedOptions.sasModel == SST_OPTIONS::SAS_TRAVIS) FTrans.resize(nPoint) = su2double(1.0);
+  if(sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) {
+    VelocityLaplacian.resize(nPoint, nDim) = su2double(0.0);
+  }
 }
 
 void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_viscosity,
diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp
index 46b686766a7..ef2bafed1c9 100644
--- a/SU2_CFD/src/variables/CTurbVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbVariable.cpp
@@ -35,4 +35,7 @@ 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);
+    SRSGridSize.resize(nPoint) = su2double(0.0);
+
    }
diff --git a/externals/codi b/externals/codi
index 762ba7698e3..8ee822a9b0b 160000
--- a/externals/codi
+++ b/externals/codi
@@ -1 +1 @@
-Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823
+Subproject commit 8ee822a9b0bb8235a2494467b774e27fb64ff14f
diff --git a/externals/medi b/externals/medi
index 7d550831e0e..aafc2d1966b 160000
--- a/externals/medi
+++ b/externals/medi
@@ -1 +1 @@
-Subproject commit 7d550831e0e233a85b9d9af9c181d7ecb2929946
+Subproject commit aafc2d1966ba1233640af737e71c77c1a86183fd
diff --git a/externals/opdi b/externals/opdi
index a6b9655c240..c42cca71a3d 160000
--- a/externals/opdi
+++ b/externals/opdi
@@ -1 +1 @@
-Subproject commit a6b9655c240af2a35454a61727e5bbbbaa3a425f
+Subproject commit c42cca71a3d0b44fb482e268ecd40b623e71776b
diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp
index c19c53ea2b8..665c45b7d35 160000
--- a/subprojects/MLPCpp
+++ b/subprojects/MLPCpp
@@ -1 +1 @@
-Subproject commit c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0
+Subproject commit 665c45b7d3533c977eb1f637918d5b8b75c07d3b