diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index c203c00285ba..f2accc9f637f 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -254,6 +254,7 @@ class CConfig {
   *Marker_Outlet,                 /*!< \brief Outlet flow markers. */
   *Marker_Isothermal,             /*!< \brief Isothermal wall markers. */
   *Marker_HeatFlux,               /*!< \brief Constant heat flux wall markers. */
+  *Marker_RoughWall,              /*!< \brief Constant heat flux wall markers. */
   *Marker_EngineInflow,           /*!< \brief Engine Inflow flow markers. */
   *Marker_EngineExhaust,          /*!< \brief Engine Exhaust flow markers. */
   *Marker_Clamped,                /*!< \brief Clamped markers. */
@@ -308,6 +309,7 @@ class CConfig {
   su2double *Outlet_Pressure;                /*!< \brief Specified back pressures (static) for outlet boundaries. */
   su2double *Isothermal_Temperature;         /*!< \brief Specified isothermal wall temperatures (static). */
   su2double *Heat_Flux;                      /*!< \brief Specified wall heat fluxes. */
+  su2double *Roughness_Height;               /*!< \brief Equivalent sand grain roughness for the marker according to config file. */
   su2double *Displ_Value;                    /*!< \brief Specified displacement for displacement boundaries. */
   su2double *Load_Value;                     /*!< \brief Specified force for load boundaries. */
   su2double *Damper_Constant;                /*!< \brief Specified constant for damper boundaries. */
@@ -555,6 +557,8 @@ class CConfig {
   *Kind_Inc_Outlet,
   *Kind_Data_Riemann,
   *Kind_Data_Giles;                /*!< \brief Kind of inlet boundary treatment. */
+  unsigned short *Kind_Wall;       /*!< \brief Type of wall treatment. */
+  unsigned short nWall_Types;      /*!< \brief Number of wall treatment types listed. */
   unsigned short nInc_Inlet;       /*!< \brief Number of inlet boundary treatment types listed. */
   unsigned short nInc_Outlet;      /*!< \brief Number of inlet boundary treatment types listed. */
   su2double Inc_Inlet_Damping;     /*!< \brief Damping factor applied to the iterative updates to the velocity at a pressure inlet in incompressible flow. */
@@ -918,7 +922,8 @@ class CConfig {
   nMarkerPitching_Ampl,           /*!< \brief Number of values provided for pitching amplitude of marker. */
   nMarkerPitching_Phase,          /*!< \brief Number of values provided for pitching phase offset of marker. */
   nMarkerPlunging_Omega,          /*!< \brief Number of values provided for angular frequency of marker. */
-  nMarkerPlunging_Ampl;           /*!< \brief Number of values provided for plunging amplitude of marker. */
+  nMarkerPlunging_Ampl,           /*!< \brief Number of values provided for plunging amplitude of marker. */
+  nRough_Wall;                    /*!< \brief Number of rough walls. */
   su2double  *Omega_HB;           /*!< \brief Frequency for Harmonic Balance Operator (in rad/s). */
   unsigned short
   nOmega_HB,                      /*!< \brief Number of frequencies in Harmonic Balance Operator. */
@@ -1077,7 +1082,8 @@ class CConfig {
   default_body_force[3],         /*!< \brief Default body force vector for the COption class. */
   default_nacelle_location[5],   /*!< \brief Location of the nacelle. */
   default_hs_axes[3],            /*!< \brief Default principal axes (x, y, z) of the ellipsoid containing the heat source. */
-  default_hs_center[3];          /*!< \brief Default position of the center of the heat source. */
+  default_hs_center[3],          /*!< \brief Default position of the center of the heat source. */
+  default_roughness[1];
 
   unsigned short Riemann_Solver_FEM;         /*!< \brief Riemann solver chosen for the DG method. */
   su2double Quadrature_Factor_Straight;      /*!< \brief Factor applied during quadrature of elements with a constant Jacobian. */
@@ -2924,6 +2930,12 @@ class CConfig {
    */
   unsigned short GetnMarker_HeatFlux(void) const { return nMarker_HeatFlux; }
 
+  /*!
+   * \brief Get the total number of rough markers.
+   * \return Total number of heat flux markers.
+   */
+  unsigned short GetnRoughWall(void) const { return nRough_Wall; }
+
   /*!
    * \brief Get the total number of objectives in kind_objective list
    * \return Total number of objectives in kind_objective list
@@ -5004,6 +5016,11 @@ class CConfig {
    */
   unsigned short GetKind_ActDisk(void) const { return Kind_ActDisk; }
 
+  /*!
+   * \brief Set the kind of wall - rough or smooth.
+   */
+  void SetKindWall(string val_marker, unsigned short val_kindwall);
+
   /*!
    * \brief Get the number of sections.
    * \return Number of sections
@@ -6754,6 +6771,13 @@ class CConfig {
    */
   su2double* GetWallFunction_DoubleInfo(string val_marker);
 
+  /*!
+   * \brief Get the type of wall and roughness height on a wall boundary (Heatflux or Isothermal).
+   * \param[in] val_index - Index corresponding to the boundary.
+   * \return The wall type and roughness height.
+   */
+  pair<unsigned short, su2double> GetWallRoughnessProperties(string val_marker) const;
+
   /*!
    * \brief Get the target (pressure, massflow, etc) at an engine inflow boundary.
    * \param[in] val_index - Index corresponding to the engine inflow boundary.
diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp
index e85f6cb384bf..f66d37ea8fbb 100644
--- a/Common/include/geometry/CPhysicalGeometry.hpp
+++ b/Common/include/geometry/CPhysicalGeometry.hpp
@@ -31,6 +31,7 @@
 #include "meshreader/CMeshReaderFVM.hpp"
 #include "../toolboxes/C2DContainer.hpp"
 
+
 /*!
  * \class CPhysicalGeometry
  * \brief Class for reading a defining the primal grid which is read from the grid file in .su2 or .cgns format.
@@ -105,6 +106,9 @@ class CPhysicalGeometry final : public CGeometry {
   unsigned long *Elem_ID_BoundTria_Linear{nullptr};
   unsigned long *Elem_ID_BoundQuad_Linear{nullptr};
 
+  vector<int> GlobalMarkerStorageDispl;
+  vector<su2double> GlobalRoughness_Height;
+
 public:
   /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/
   using CGeometry::SetVertex;
@@ -804,4 +808,9 @@ class CPhysicalGeometry final : public CGeometry {
     }
   }
 
+  /*!
+   * \brief Set roughness values for markers in a global array.
+   */
+  void SetGlobalMarkerRoughness(const CConfig* config);
+
 };
diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp
index 2e0d14530e40..b211baaa401d 100644
--- a/Common/include/geometry/dual_grid/CPoint.hpp
+++ b/Common/include/geometry/dual_grid/CPoint.hpp
@@ -85,6 +85,7 @@ class CPoint {
   su2activevector SharpEdge_Distance;     /*!< \brief Distance to a sharp edge. */
   su2activevector Curvature;              /*!< \brief Value of the surface curvature (SU2_GEO). */
   su2activevector MaxLength;              /*!< \brief The maximum cell-center to cell-center length. */
+  su2activevector RoughnessHeight;          /*!< \brief Roughness of the nearest wall. */
 
   su2matrix<int> AD_InputIndex;           /*!< \brief Indices of Coord variables in the adjoint vector. */
   su2matrix<int> AD_OutputIndex;          /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */
@@ -400,6 +401,20 @@ class CPoint {
    */
   inline su2double GetWall_Distance(unsigned long iPoint) const { return Wall_Distance(iPoint); }
 
+  /*!
+   * \brief Set the value of the distance to the nearest wall.
+   * \param[in] iPoint - Index of the point.
+   * \param[in] distance - Value of the distance.
+   */
+  inline void SetRoughnessHeight(unsigned long iPoint, su2double roughheight) { RoughnessHeight(iPoint) = roughheight; }
+
+  /*!
+   * \brief Get the value of the distance to the nearest wall.
+   * \param[in] iPoint - Index of the point.
+   * \return Value of the distance to the nearest wall.
+   */
+  inline su2double GetRoughnessHeight(unsigned long iPoint) const { return RoughnessHeight(iPoint); }
+
   /*!
    * \brief Set the value of the distance to a sharp edge.
    * \param[in] iPoint - Index of the point.
diff --git a/Common/include/geometry/dual_grid/CTurboVertex.hpp b/Common/include/geometry/dual_grid/CTurboVertex.hpp
index 498cfcf1ea7d..8b36513ec9b8 100644
--- a/Common/include/geometry/dual_grid/CTurboVertex.hpp
+++ b/Common/include/geometry/dual_grid/CTurboVertex.hpp
@@ -117,7 +117,7 @@ class CTurboVertex final : public CVertex {
    * \brief get global index for ordered span-wise turbovertex.
    */
   inline int GetGlobalVertexIndex(void) const {return GlobalIndex;}
- 
+
   /*!
    * \brief set angular coord.
    */
diff --git a/Common/include/mpi_structure.hpp b/Common/include/mpi_structure.hpp
index 58d4d58f9ae3..fc41129cd542 100644
--- a/Common/include/mpi_structure.hpp
+++ b/Common/include/mpi_structure.hpp
@@ -412,7 +412,7 @@ class CBaseMPIWrapper {
                         void *recvbuf, int recvcnt, Datatype recvtype, Comm comm);
 
   static void Allgatherv(void *sendbuf, int sendcnt, Datatype sendtype,
-                         void *recvbuf, int recvcnt, int *displs, Datatype recvtype, Comm comm);
+                         void *recvbuf, int *recvcnt, int *displs, Datatype recvtype, Comm comm);
 
   static void Sendrecv(void *sendbuf, int sendcnt, Datatype sendtype,
                        int dest, int sendtag, void *recvbuf, int recvcnt,
diff --git a/Common/include/mpi_structure.inl b/Common/include/mpi_structure.inl
index ac3c681dc8c7..6ecc766e0bdc 100644
--- a/Common/include/mpi_structure.inl
+++ b/Common/include/mpi_structure.inl
@@ -614,7 +614,7 @@ inline void CBaseMPIWrapper::Scatter(void *sendbuf, int sendcnt, Datatype sendty
 }
 
 inline void CBaseMPIWrapper::Allgatherv(void *sendbuf, int sendcnt, Datatype sendtype,
-                                   void *recvbuf, int recvcnt, int *displs, Datatype recvtype, Comm comm){
+                                   void *recvbuf, int *recvcnt, int *displs, Datatype recvtype, Comm comm){
   CopyData(sendbuf, recvbuf, sendcnt, sendtype);
 }
 
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 597a0886a012..a0ba49dda6df 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1394,6 +1394,18 @@ static const MapType<string, ACTDISK_TYPE> ActDisk_Map = {
   MakePair("VARIABLE_LOAD", VARIABLE_LOAD)
 };
 
+/*!
+ * \brief types of wall boundary condition - smooth or rough
+ */
+enum WALL_TYPE {
+  SMOOTH = 1,    /*!< \brief Smooth wall */
+  ROUGH = 2,   /*!< \brief Rough wall */
+};
+static const MapType<string, WALL_TYPE> WallType_Map = {
+  MakePair("SMOOTH", SMOOTH)
+  MakePair("ROUGH", ROUGH)
+};
+
 /*!
  * \brief Types of geometric entities based on VTK nomenclature
  */
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index c5db0cd8367b..920f36173fb7 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -806,6 +806,7 @@ void CConfig::SetPointersNull(void) {
   Kind_WallFunctions       = nullptr;
   IntInfo_WallFunctions    = nullptr;
   DoubleInfo_WallFunctions = nullptr;
+  Kind_Wall                = nullptr;
 
   Config_Filenames = nullptr;
 
@@ -818,7 +819,7 @@ void CConfig::SetPointersNull(void) {
   Marker_Inlet                = nullptr;    Marker_Outlet           = nullptr;
   Marker_Supersonic_Inlet     = nullptr;    Marker_Supersonic_Outlet= nullptr;
   Marker_Isothermal           = nullptr;    Marker_HeatFlux         = nullptr;    Marker_EngineInflow   = nullptr;
-  Marker_Load                 = nullptr;    Marker_Disp_Dir         = nullptr;
+  Marker_Load                 = nullptr;    Marker_Disp_Dir         = nullptr;    Marker_RoughWall      = nullptr;
   Marker_EngineExhaust        = nullptr;    Marker_Displacement     = nullptr;    Marker_Load           = nullptr;
   Marker_Load_Dir             = nullptr;    Marker_Load_Sine        = nullptr;    Marker_Clamped        = nullptr;
   Marker_FlowLoad             = nullptr;    Marker_Internal         = nullptr;
@@ -831,8 +832,9 @@ void CConfig::SetPointersNull(void) {
     /*--- Boundary Condition settings ---*/
 
   Isothermal_Temperature = nullptr;
-  Heat_Flux              = nullptr;    Displ_Value            = nullptr;    Load_Value      = nullptr;
-  FlowLoad_Value         = nullptr;    Damper_Constant        = nullptr;    Wall_Emissivity = nullptr;
+  Heat_Flux              = nullptr;    Displ_Value            = nullptr;    Load_Value             = nullptr;
+  FlowLoad_Value         = nullptr;    Damper_Constant        = nullptr;    Wall_Emissivity        = nullptr;
+  Roughness_Height       = nullptr;
 
   /*--- Inlet Outlet Boundary Condition settings ---*/
 
@@ -863,8 +865,6 @@ void CConfig::SetPointersNull(void) {
   Inlet_FlowDir             = nullptr;     Inlet_Temperature           = nullptr;     Inlet_Pressure        = nullptr;
   Inlet_Velocity            = nullptr;     Inflow_Mach                 = nullptr;     Inflow_Pressure       = nullptr;
   Exhaust_Pressure          = nullptr;     Outlet_Pressure             = nullptr;     Isothermal_Temperature= nullptr;
-  Heat_Flux                 = nullptr;     Displ_Value                 = nullptr;     Load_Value            = nullptr;
-  FlowLoad_Value            = nullptr;
 
   ElasticityMod             = nullptr;     PoissonRatio                = nullptr;     MaterialDensity       = nullptr;
 
@@ -1499,6 +1499,9 @@ void CConfig::SetConfig_Options() {
   /*!\brief MARKER_HEATFLUX  \n DESCRIPTION: Specified heat flux wall boundary marker(s)
    Format: ( Heat flux marker, wall heat flux (static), ... ) \ingroup Config*/
   addStringDoubleListOption("MARKER_HEATFLUX", nMarker_HeatFlux, Marker_HeatFlux, Heat_Flux);
+  /*!\brief WALL_ROUGHNESS  \n DESCRIPTION: Specified roughness heights at wall boundary marker(s)
+   Format: ( Wall marker, roughness_height (static), ... ) \ingroup Config*/
+  addStringDoubleListOption("WALL_ROUGHNESS", nRough_Wall, Marker_RoughWall, Roughness_Height);
   /*!\brief MARKER_ENGINE_INFLOW  \n DESCRIPTION: Engine inflow boundary marker(s)
    Format: ( nacelle inflow marker, fan face Mach, ... ) \ingroup Config*/
   addStringDoubleListOption("MARKER_ENGINE_INFLOW", nMarker_EngineInflow, Marker_EngineInflow, EngineInflow_Target);
@@ -4602,6 +4605,86 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
     }
   }
 
+  /*--- Check that if the wall roughness array are compatible and set deafult values if needed. ---*/
+   if ((nMarker_HeatFlux > 0) || (nMarker_Isothermal > 0) || (nMarker_CHTInterface > 0)) {
+
+     /*--- The total number of wall markers. ---*/
+     unsigned short nWall = nMarker_HeatFlux + nMarker_Isothermal + nMarker_CHTInterface;
+
+     /*--- If no roughness is specified all walls are assumed to be smooth. ---*/
+     if (nRough_Wall == 0) {
+
+       nRough_Wall = nWall;
+       Roughness_Height = new su2double [nWall];
+       Kind_Wall = new unsigned short [nWall];
+       for (unsigned short iMarker = 0; iMarker < nMarker_HeatFlux; iMarker++) {
+         Roughness_Height[iMarker] = 0.0;
+         Kind_Wall[iMarker] = SMOOTH;
+       }
+       for (iMarker = 0; iMarker < nMarker_Isothermal; iMarker++) {
+         Roughness_Height[nMarker_HeatFlux + iMarker] = 0.0;
+         Kind_Wall[nMarker_HeatFlux + iMarker] = SMOOTH;
+       }
+       for (iMarker = 0; iMarker < nMarker_CHTInterface; iMarker++) {
+         Roughness_Height[nMarker_HeatFlux + nMarker_Isothermal + iMarker] = 0.0;
+         Kind_Wall[nMarker_HeatFlux + nMarker_Isothermal + iMarker] = SMOOTH;
+       }
+
+       /*--- Check for mismatch in number of rough walls and solid walls. ---*/
+     } else if (nRough_Wall > nWall) {
+        SU2_MPI::Error("Mismatch in number of rough walls and solid walls. Number of rough walls cannot be more than solid walls.", CURRENT_FUNCTION);
+       /*--- Check name of the marker and assign the corresponding roughness. ---*/
+     } else {
+       /*--- Store roughness heights in a temp array. ---*/
+       vector<su2double> temp_rough;
+       for (iMarker = 0; iMarker < nRough_Wall; iMarker++)
+         temp_rough.push_back(Roughness_Height[iMarker]);
+
+       /*--- Reallocate the roughness arrays in case not all walls are rough. ---*/
+       delete Roughness_Height;
+       delete Kind_Wall;
+       Roughness_Height = new su2double [nWall];
+       Kind_Wall = new unsigned short [nWall];
+       unsigned short jMarker, chkRough = 0;
+
+       /*--- Initialize everything to smooth. ---*/
+       for (iMarker = 0; iMarker < nWall; iMarker++) {
+         Roughness_Height[iMarker] = 0.0;
+         Kind_Wall[iMarker] = SMOOTH;
+       }
+
+       /*--- Look through heat flux, isothermal and cht_interface markers and assign proper values. ---*/
+       for (iMarker = 0; iMarker < nRough_Wall; iMarker++) {
+         for (jMarker = 0; jMarker < nMarker_HeatFlux; jMarker++)
+           if (Marker_HeatFlux[jMarker].compare(Marker_RoughWall[iMarker]) == 0) {
+             Roughness_Height[jMarker] = temp_rough[iMarker];
+             chkRough++;
+           }
+
+         for (jMarker = 0; jMarker < nMarker_Isothermal; jMarker++)
+           if (Marker_Isothermal[jMarker].compare(Marker_RoughWall[iMarker]) == 0) {
+             Roughness_Height[nMarker_HeatFlux + jMarker] = temp_rough[iMarker];
+             chkRough++;
+           }
+
+         for (jMarker = 0; jMarker < nMarker_CHTInterface; jMarker++)
+           if (Marker_CHTInterface[jMarker].compare(Marker_RoughWall[iMarker]) == 0) {
+             Roughness_Height[nMarker_HeatFlux + nMarker_Isothermal + jMarker] = temp_rough[iMarker];
+             chkRough++;
+           }
+       }
+
+       /*--- Update kind_wall when a non zero roughness value is specified. ---*/
+       for (iMarker = 0; iMarker < nWall; iMarker++)
+         if (Roughness_Height[iMarker] != 0.0)
+           Kind_Wall[iMarker] = ROUGH;
+
+       /*--- Check if a non solid wall marker was specified as rough. ---*/
+       if (chkRough != nRough_Wall)
+         SU2_MPI::Error("Only solid walls can be rough.", CURRENT_FUNCTION);
+     }
+   }
+
   /*--- Handle default options for topology optimization ---*/
 
   if (topology_optimization && top_optim_nKernel==0) {
@@ -7352,23 +7435,25 @@ CConfig::~CConfig(void) {
   delete [] Marker_CfgFile_MixingPlaneInterface;
   delete [] Marker_All_MixingPlaneInterface;
 
-                delete[] Marker_DV;
-            delete[] Marker_Moving;
-       delete[] Marker_Monitoring;
-        delete[] Marker_Designing;
-          delete[] Marker_GeoEval;
-         delete[] Marker_Plotting;
-         delete[] Marker_Analyze;
-   delete[] Marker_WallFunctions;
-         delete[] Marker_ZoneInterface;
-              delete [] Marker_PyCustom;
-     delete[] Marker_All_SendRecv;
-
-       delete[] Kind_Inc_Inlet;
-       delete[] Kind_Inc_Outlet;
+  delete[] Marker_DV;
+  delete[] Marker_Moving;
+  delete[] Marker_Monitoring;
+  delete[] Marker_Designing;
+  delete[] Marker_GeoEval;
+  delete[] Marker_Plotting;
+  delete[] Marker_Analyze;
+  delete[] Marker_WallFunctions;
+  delete[] Marker_ZoneInterface;
+  delete [] Marker_PyCustom;
+  delete[] Marker_All_SendRecv;
+
+  delete[] Kind_Inc_Inlet;
+  delete[] Kind_Inc_Outlet;
 
   delete[] Kind_WallFunctions;
 
+  delete[] Kind_Wall;
+
   delete[] Config_Filenames;
 
   if (IntInfo_WallFunctions != nullptr) {
@@ -7554,8 +7639,8 @@ CConfig::~CConfig(void) {
      delete[] Load_Sine_Amplitude;
      delete[] Load_Sine_Frequency;
      delete[] FlowLoad_Value;
+     delete[] Roughness_Height;
      delete[] Wall_Emissivity;
-
   /*--- related to periodic boundary conditions ---*/
 
   for (iMarker = 0; iMarker < nMarker_PerBound; iMarker++) {
@@ -8666,7 +8751,38 @@ su2double CConfig::GetWall_HeatFlux(string val_marker) const {
   return Heat_Flux[iMarker_HeatFlux];
 }
 
+pair<unsigned short, su2double> CConfig::GetWallRoughnessProperties(string val_marker) const {
+  unsigned short iMarker = 0;
+  short          flag = -1;
+  pair<unsigned short, su2double> WallProp;
+
+  if (nMarker_HeatFlux > 0 || nMarker_Isothermal > 0 || nMarker_CHTInterface > 0) {
+    for (iMarker = 0; iMarker < nMarker_HeatFlux; iMarker++)
+      if (Marker_HeatFlux[iMarker] == val_marker) {
+        flag = iMarker;
+        WallProp = make_pair(Kind_Wall[flag], Roughness_Height[flag]);
+        return WallProp;
+      }
+    for (iMarker = 0; iMarker < nMarker_Isothermal; iMarker++)
+      if (Marker_Isothermal[iMarker] == val_marker) {
+        flag = nMarker_HeatFlux + iMarker;
+        WallProp = make_pair(Kind_Wall[flag], Roughness_Height[flag]);
+        return WallProp;
+      }
+    for (iMarker = 0; iMarker < nMarker_CHTInterface; iMarker++)
+      if (Marker_CHTInterface[iMarker] == val_marker) {
+        flag = nMarker_HeatFlux + nMarker_Isothermal + iMarker;
+        WallProp = make_pair(Kind_Wall[flag], Roughness_Height[flag]);
+        return WallProp;
+      }
+  }
+
+  WallProp = make_pair(SMOOTH, 0.0);
+  return WallProp;
+}
+
 unsigned short CConfig::GetWallFunction_Treatment(string val_marker) const {
+
   unsigned short WallFunction = NO_WALL_FUNCTION;
 
   for(unsigned short iMarker=0; iMarker<nMarker_WallFunctions; iMarker++) {
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
index 8cf5a7709042..ba277a6971fc 100644
--- a/Common/src/geometry/CPhysicalGeometry.cpp
+++ b/Common/src/geometry/CPhysicalGeometry.cpp
@@ -156,7 +156,7 @@ CPhysicalGeometry::CPhysicalGeometry(CConfig *config, unsigned short val_iZone,
           if (bound[iMarker][iElem_Bound]->GetVTK_Type() == VERTEX) {
             boundary_file << bound[iMarker][iElem_Bound]->GetRotation_Type() << "\t";
           }
-          boundary_file	<< iElem_Bound << endl;
+          boundary_file << iElem_Bound << endl;
         }
       }
 
@@ -169,7 +169,7 @@ CPhysicalGeometry::CPhysicalGeometry(CConfig *config, unsigned short val_iZone,
           if (bound[iMarker][iElem_Bound]->GetVTK_Type() == VERTEX) {
             boundary_file << bound[iMarker][iElem_Bound]->GetRotation_Type() << "\t";
           }
-          boundary_file	<< iElem_Bound << endl;
+          boundary_file << iElem_Bound << endl;
         }
       }
 
@@ -3331,7 +3331,7 @@ void CPhysicalGeometry::SetSendReceive(CConfig *config) {
   }
 
   /*--- First compute the Send/Receive boundaries ---*/
-  Counter_Send = 0; 	Counter_Receive = 0;
+  Counter_Send = 0;   Counter_Receive = 0;
   for (iDomain = 0; iDomain < nDomain; iDomain++)
     if (SendDomainLocal[iDomain].size() != 0) Counter_Send++;
 
@@ -4785,7 +4785,7 @@ void CPhysicalGeometry::SetPositive_ZArea(CConfig *config) {
     if (config->GetSystemMeasurements() == SI) cout <<" m"; else cout <<" ft";
 
     if (nDim == 3) {
-    	cout << ", z-direction = "<< TotalMaxCoordZ;
+      cout << ", z-direction = "<< TotalMaxCoordZ;
       if (config->GetSystemMeasurements() == SI) cout <<" m." << endl; else cout <<" ft."<< endl;
     }
     else cout << "." << endl;
@@ -4797,7 +4797,7 @@ void CPhysicalGeometry::SetPositive_ZArea(CConfig *config) {
     if (config->GetSystemMeasurements() == SI) cout <<" m"; else cout <<" ft";
 
     if (nDim == 3) {
-    	cout << ", z-direction = "<< TotalMinCoordZ;
+      cout << ", z-direction = "<< TotalMinCoordZ;
       if (config->GetSystemMeasurements() == SI) cout <<" m." << endl; else cout <<" ft."<< endl;
     }
     else cout << "." << endl;
@@ -5166,7 +5166,7 @@ void CPhysicalGeometry::SetVertex(CConfig *config) {
 }
 
 void CPhysicalGeometry::ComputeNSpan(CConfig *config, unsigned short val_iZone, unsigned short marker_flag, bool allocate) {
-  unsigned short iMarker, jMarker, iMarkerTP, iSpan, jSpan, kSpan = 0;
+unsigned short iMarker, jMarker, iMarkerTP, iSpan, jSpan, kSpan = 0;
   unsigned long iPoint, iVertex;
   long jVertex;
   int nSpan, nSpan_loc;
@@ -6121,7 +6121,7 @@ void CPhysicalGeometry::SetTurboVertex(CConfig *config, unsigned short val_iZone
         }
       }
 
-      delete [] x_gb;	delete [] y_gb; delete [] z_gb;	 delete [] angCoord_gb; delete [] deltaAngCoord_gb; delete[] checkAssign_gb;
+      delete [] x_gb; delete [] y_gb; delete [] z_gb;  delete [] angCoord_gb; delete [] deltaAngCoord_gb; delete[] checkAssign_gb;
 
     }
   }
@@ -6407,7 +6407,7 @@ void CPhysicalGeometry::SetAvgTurboValue(CConfig *config, unsigned short val_iZo
 
     /*--- Forces initialization for contenitors to zero ---*/
     for (iDim=0; iDim<nDim; iDim++) {
-      TotalTurboNormal[iDim]	=0.0;
+      TotalTurboNormal[iDim]  =0.0;
       TotalNormal[iDim]         =0.0;
       TotalGridVel[iDim]        =0.0;
     }
@@ -6458,7 +6458,7 @@ void CPhysicalGeometry::SetAvgTurboValue(CConfig *config, unsigned short val_iZo
     MyTotalGridVel         = new su2double[nDim];
 
     for (iDim = 0; iDim < nDim; iDim++) {
-      MyTotalTurboNormal[iDim]	= TotalTurboNormal[iDim];
+      MyTotalTurboNormal[iDim]  = TotalTurboNormal[iDim];
       TotalTurboNormal[iDim]    = 0.0;
       MyTotalNormal[iDim]       = TotalNormal[iDim];
       TotalNormal[iDim]         = 0.0;
@@ -6559,7 +6559,7 @@ void CPhysicalGeometry::SetAvgTurboValue(CConfig *config, unsigned short val_iZo
       if (config->GetMarker_All_Turbomachinery(iMarker) == iMarkerTP){
         if (config->GetMarker_All_TurbomachineryFlag(iMarker) == marker_flag){
           turboNormal2 = 0.0;
-          Normal2 		= 0.0;
+          Normal2     = 0.0;
 
           for (iDim = 0; iDim < nDim; iDim++){
             turboNormal2 += AverageTurboNormal[iMarker][nSpanWiseSections[marker_flag-1]][iDim]*AverageTurboNormal[iMarker][nSpanWiseSections[marker_flag-1]][iDim];
@@ -8933,7 +8933,7 @@ void CPhysicalGeometry::SetBoundSensitivity(CConfig *config) {
 
     /*--- Write file name with extension if unsteady or steady ---*/
     if (config->GetTime_Marching() == HARMONIC_BALANCE)
-    	SPRINTF (buffer, "_%d.csv", SU2_TYPE::Int(iTimeIter));
+      SPRINTF (buffer, "_%d.csv", SU2_TYPE::Int(iTimeIter));
 
     if ((config->GetTime_Marching() && config->GetTime_Domain()) ||
         (config->GetTime_Marching() == HARMONIC_BALANCE)) {
@@ -9755,7 +9755,7 @@ su2double CPhysicalGeometry::Compute_MaxThickness(su2double *Plane_P0, su2double
     unsigned short index = 2;
 
     /*--- Removing the trailing edge from list of points that we are going to use in the interpolation,
-			to be sure that a blunt trailing edge do not affect the interpolation ---*/
+      to be sure that a blunt trailing edge do not affect the interpolation ---*/
 
     if ((Normal[index] >= 0.0) && (fabs(Xcoord_Airfoil_[iVertex]) > MaxDistance*0.01)) {
       Xcoord.push_back(Xcoord_Airfoil_[iVertex]);
@@ -9981,8 +9981,8 @@ su2double CPhysicalGeometry::Compute_WaterLineWidth(su2double *Plane_P0, su2doub
   for (iVertex = 0; iVertex < Xcoord_Airfoil.size(); iVertex++) {
     Distance = fabs(Zcoord_Airfoil[iVertex] - WaterLine);
     if (Distance < MinDistance) {
-    	MinDistance = Distance;
-    	WaterLineWidth = fabs(Xcoord_Airfoil[iVertex] - Xcoord_Airfoil[Trailing_Point]);
+      MinDistance = Distance;
+      WaterLineWidth = fabs(Xcoord_Airfoil[iVertex] - Xcoord_Airfoil[Trailing_Point]);
     }
   }
 
@@ -10441,21 +10441,21 @@ void CPhysicalGeometry::Compute_Wing(CConfig *config, bool original_surface,
     /*--- Plot the geometrical quatities ---*/
 
     for (iPlane = 0; iPlane < nPlane; iPlane++) {
-    	if (Xcoord_Airfoil[iPlane].size() > 1) {
-    		if (config->GetTabular_FileFormat() == TAB_CSV) {
-    			Wing_File  << Ycoord_Airfoil[iPlane][0]/SemiSpan <<", "<< Area[iPlane] <<", "<< MaxThickness[iPlane] <<", "<< Chord[iPlane] <<", "<< LERadius[iPlane] <<", "<< ToC[iPlane]
-    			           <<", "<< Twist[iPlane] <<", "<< Curvature[iPlane] <<", "<< Dihedral[iPlane]
-    			           <<", "<< LeadingEdge[iPlane][0]/SemiSpan <<", "<< LeadingEdge[iPlane][2]/SemiSpan
-    			           <<", "<< TrailingEdge[iPlane][0]/SemiSpan <<", "<< TrailingEdge[iPlane][2]/SemiSpan << endl;
-    		}
-    		else  {
-    			Wing_File  << Ycoord_Airfoil[iPlane][0]/SemiSpan <<" "<< Area[iPlane] <<" "<< MaxThickness[iPlane] <<" "<< Chord[iPlane] <<" "<< LERadius[iPlane] <<" "<< ToC[iPlane]
-    			           <<" "<< Twist[iPlane] <<" "<< Curvature[iPlane]  <<" "<< Dihedral[iPlane]
-    			           <<" "<< LeadingEdge[iPlane][0]/SemiSpan <<" "<< LeadingEdge[iPlane][2]/SemiSpan
-    			           <<" "<< TrailingEdge[iPlane][0]/SemiSpan <<" "<< TrailingEdge[iPlane][2]/SemiSpan << endl;
-
-    		}
-    	}
+      if (Xcoord_Airfoil[iPlane].size() > 1) {
+        if (config->GetTabular_FileFormat() == TAB_CSV) {
+          Wing_File  << Ycoord_Airfoil[iPlane][0]/SemiSpan <<", "<< Area[iPlane] <<", "<< MaxThickness[iPlane] <<", "<< Chord[iPlane] <<", "<< LERadius[iPlane] <<", "<< ToC[iPlane]
+                     <<", "<< Twist[iPlane] <<", "<< Curvature[iPlane] <<", "<< Dihedral[iPlane]
+                     <<", "<< LeadingEdge[iPlane][0]/SemiSpan <<", "<< LeadingEdge[iPlane][2]/SemiSpan
+                     <<", "<< TrailingEdge[iPlane][0]/SemiSpan <<", "<< TrailingEdge[iPlane][2]/SemiSpan << endl;
+        }
+        else  {
+          Wing_File  << Ycoord_Airfoil[iPlane][0]/SemiSpan <<" "<< Area[iPlane] <<" "<< MaxThickness[iPlane] <<" "<< Chord[iPlane] <<" "<< LERadius[iPlane] <<" "<< ToC[iPlane]
+                     <<" "<< Twist[iPlane] <<" "<< Curvature[iPlane]  <<" "<< Dihedral[iPlane]
+                     <<" "<< LeadingEdge[iPlane][0]/SemiSpan <<" "<< LeadingEdge[iPlane][2]/SemiSpan
+                     <<" "<< TrailingEdge[iPlane][0]/SemiSpan <<" "<< TrailingEdge[iPlane][2]/SemiSpan << endl;
+
+        }
+      }
     }
 
     Wing_File.close();
@@ -10727,18 +10727,18 @@ void CPhysicalGeometry::Compute_Fuselage(CConfig *config, bool original_surface,
     /*--- Plot the geometrical quatities ---*/
 
     for (iPlane = 0; iPlane < nPlane; iPlane++) {
-    	if (Xcoord_Airfoil[iPlane].size() > 1) {
-    		if (config->GetTabular_FileFormat() == TAB_CSV) {
-    			Fuselage_File  << -Ycoord_Airfoil[iPlane][0] <<", "<< Area[iPlane] <<", "<< Length[iPlane] <<", "<< Width[iPlane] <<", "<< WaterLineWidth[iPlane] <<", "<< Height[iPlane] <<", "<< Curvature[iPlane]
-    			           <<", "<< -LeadingEdge[iPlane][1] <<", "<< LeadingEdge[iPlane][0]  <<", "<< LeadingEdge[iPlane][2]
-    			           <<", "<< -TrailingEdge[iPlane][1] <<", "<< TrailingEdge[iPlane][0]  <<", "<< TrailingEdge[iPlane][2]  << endl;
-    		}
-    		else  {
-    			Fuselage_File  << -Ycoord_Airfoil[iPlane][0] <<" "<< Area[iPlane] <<" "<< Length[iPlane] <<" "<< Width[iPlane] <<" "<< WaterLineWidth[iPlane] <<" "<< Height[iPlane] <<" "<< Curvature[iPlane]
-    			           <<" "<< -LeadingEdge[iPlane][1] <<" "<< LeadingEdge[iPlane][0]  <<" "<< LeadingEdge[iPlane][2]
-    			           <<" "<< -TrailingEdge[iPlane][1] <<" "<< TrailingEdge[iPlane][0]  <<" "<< TrailingEdge[iPlane][2] << endl;
-    		}
-    	}
+      if (Xcoord_Airfoil[iPlane].size() > 1) {
+        if (config->GetTabular_FileFormat() == TAB_CSV) {
+          Fuselage_File  << -Ycoord_Airfoil[iPlane][0] <<", "<< Area[iPlane] <<", "<< Length[iPlane] <<", "<< Width[iPlane] <<", "<< WaterLineWidth[iPlane] <<", "<< Height[iPlane] <<", "<< Curvature[iPlane]
+                     <<", "<< -LeadingEdge[iPlane][1] <<", "<< LeadingEdge[iPlane][0]  <<", "<< LeadingEdge[iPlane][2]
+                     <<", "<< -TrailingEdge[iPlane][1] <<", "<< TrailingEdge[iPlane][0]  <<", "<< TrailingEdge[iPlane][2]  << endl;
+        }
+        else  {
+          Fuselage_File  << -Ycoord_Airfoil[iPlane][0] <<" "<< Area[iPlane] <<" "<< Length[iPlane] <<" "<< Width[iPlane] <<" "<< WaterLineWidth[iPlane] <<" "<< Height[iPlane] <<" "<< Curvature[iPlane]
+                     <<" "<< -LeadingEdge[iPlane][1] <<" "<< LeadingEdge[iPlane][0]  <<" "<< LeadingEdge[iPlane][2]
+                     <<" "<< -TrailingEdge[iPlane][1] <<" "<< TrailingEdge[iPlane][0]  <<" "<< TrailingEdge[iPlane][2] << endl;
+        }
+      }
     }
 
     Fuselage_File.close();
@@ -10776,7 +10776,7 @@ void CPhysicalGeometry::Compute_Fuselage(CConfig *config, bool original_surface,
     Fuselage_Volume = 0.0;
     for (iPlane = 0; iPlane < nPlane-2; iPlane+=2) {
       if (Xcoord_Airfoil[iPlane].size() > 1) {
-      	Fuselage_Volume += (1.0/3.0)*dPlane*(Area[iPlane] + 4.0*Area[iPlane+1] + Area[iPlane+2]);
+        Fuselage_Volume += (1.0/3.0)*dPlane*(Area[iPlane] + 4.0*Area[iPlane+1] + Area[iPlane+2]);
       }
     }
 
@@ -10786,7 +10786,7 @@ void CPhysicalGeometry::Compute_Fuselage(CConfig *config, bool original_surface,
     if (Xcoord_Airfoil[0].size() > 1) Fuselage_WettedArea += (1.0/2.0)*dPlane*Length[0];
     for (iPlane = 1; iPlane < nPlane-1; iPlane++) {
       if (Xcoord_Airfoil[iPlane].size() > 1) {
-      	Fuselage_WettedArea += dPlane*Length[iPlane];
+        Fuselage_WettedArea += dPlane*Length[iPlane];
       }
     }
     if (Xcoord_Airfoil[nPlane-1].size() > 1) Fuselage_WettedArea += (1.0/2.0)*dPlane*Length[nPlane-1];
@@ -11193,7 +11193,6 @@ std::unique_ptr<CADTElemClass> CPhysicalGeometry::ComputeViscousWallADT(const CC
         markerIDs.push_back(iMarker);
         VTK_TypeElem.push_back(VTK_Type);
         elemIDs.push_back(iElem);
-
         for (unsigned short iNode = 0; iNode < nDOFsPerElem; iNode++)
           surfaceConn.push_back(bound[iMarker][iElem]->GetNode(iNode));
       }
@@ -11241,6 +11240,9 @@ void CPhysicalGeometry::SetWallDistance(const CConfig *config, CADTElemClass *Wa
   /*---        distance to a solid wall element                           ---*/
   /*--------------------------------------------------------------------------*/
 
+  /*--- Store marker list and roughness in a global array. ---*/
+  if (config->GetnRoughWall() > 0) SetGlobalMarkerRoughness(config);
+
   SU2_OMP_PARALLEL
   if (!WallADT->IsEmpty()) {
     /*--- Solid wall boundary nodes are present. Compute the wall
@@ -11252,9 +11254,51 @@ void CPhysicalGeometry::SetWallDistance(const CConfig *config, CADTElemClass *Wa
       unsigned long  elemID;
       int            rankID;
       su2double      dist;
+
       WallADT->DetermineNearestElement(nodes->GetCoord(iPoint), dist, markerID, elemID, rankID);
+
       nodes->SetWall_Distance(iPoint, min(dist,nodes->GetWall_Distance(iPoint)));
+
+      if (config->GetnRoughWall() > 0) {
+        auto index = GlobalMarkerStorageDispl[rankID] + markerID;
+        auto localRoughness = GlobalRoughness_Height[index];
+        nodes->SetRoughnessHeight(iPoint, localRoughness);
+      }
     }
+
   }
   // end SU2_OMP_PARALLEL
 }
+
+void CPhysicalGeometry::SetGlobalMarkerRoughness(const CConfig* config) {
+
+  const auto nMarker_All = config->GetnMarker_All();
+
+  vector<int> recvCounts(size);
+  auto sizeLocal = static_cast<int>(nMarker_All); // number of local markers
+
+  /*--- Communicate size of local marker array and make an array large enough to hold all data. ---*/
+  SU2_MPI::Allgather(&sizeLocal, 1, MPI_INT, recvCounts.data(), 1, MPI_INT, MPI_COMM_WORLD);
+
+  /*--- Set the global array of displacements, needed to access the correct roughness element. ---*/
+  GlobalMarkerStorageDispl.resize(size);
+  GlobalMarkerStorageDispl[0] = 0;
+  for (int iRank = 1; iRank < size; iRank++)
+    GlobalMarkerStorageDispl[iRank] = GlobalMarkerStorageDispl[iRank-1] + recvCounts[iRank-1];
+
+  /*--- Total size ---*/
+  const auto sizeGlobal = GlobalMarkerStorageDispl[size-1] + recvCounts[size-1];
+
+  /*--- Allocate local and global arrays to hold roughness. ---*/
+  vector<su2double> localRough(nMarker_All);   // local number of markers
+  GlobalRoughness_Height.resize(sizeGlobal);   // all markers including send recieve
+
+  for (auto iMarker = 0u; iMarker < nMarker_All; iMarker++) {
+    auto wallprop = config->GetWallRoughnessProperties(config->GetMarker_All_TagBound(iMarker));
+    localRough[iMarker] = wallprop.second;
+  }
+
+  /*--- Finally, gather the roughness of all markers. ---*/
+  SU2_MPI::Allgatherv(localRough.data(), sizeLocal, MPI_DOUBLE, GlobalRoughness_Height.data(),
+                      recvCounts.data(), GlobalMarkerStorageDispl.data(), MPI_DOUBLE, MPI_COMM_WORLD);
+}
diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp
index a7e5099bb496..fb8fb1bff814 100644
--- a/Common/src/geometry/dual_grid/CPoint.cpp
+++ b/Common/src/geometry/dual_grid/CPoint.cpp
@@ -126,6 +126,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig *config) {
   MaxLength.resize(npoint) = su2double(0.0);
   Curvature.resize(npoint) = su2double(0.0);
   Wall_Distance.resize(npoint) = su2double(0.0);
+  RoughnessHeight.resize(npoint) = su2double(0.0);
   SharpEdge_Distance.resize(npoint) = su2double(0.0);
 
 }
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index 696c202de82e..afd8310eb5ec 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -205,6 +205,8 @@ class CNumerics {
   su2double StrainMag_i, StrainMag_j;      /*!< \brief Strain rate magnitude. */
   su2double Dissipation_i, Dissipation_j;  /*!< \brief Dissipation. */
   su2double Dissipation_ij;
+  su2double roughness_i,                   /*!< \brief Roughness of the wall nearest to point i. */
+  roughness_j;                             /*!< \brief Roughness of the wall nearest to point j. */
 
   su2double *l, *m;
 
@@ -557,6 +559,16 @@ class CNumerics {
     dist_j = val_dist_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
+   * \param[in] val_dist_j - Value of of the roughness of the nearest wall from point j
+   */
+  void SetRoughness(su2double val_roughness_i, su2double val_roughness_j) {
+    roughness_i = val_roughness_i;
+    roughness_j = val_roughness_j;
+  }
+
   /*!
    * \brief Set coordinates of the points.
    * \param[in] val_coord_i - Coordinates of the point i.
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 26f1f6e40c33..9c562a7ccb9c 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -49,6 +49,7 @@ class CSourceBase_TurbSA : public CNumerics {
   su2double sigma;
   su2double cb2;
   su2double cw1;
+  su2double cr1;
 
   su2double gamma_BC;
   su2double intermittency;
@@ -60,6 +61,7 @@ class CSourceBase_TurbSA : public CNumerics {
 
 protected:
   const bool incompressible = false, rotating_frame = false;
+  bool roughwall = false;
 
 public:
   /*!
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index c6373fa0361e..8ffbec0648af 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -1974,7 +1974,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
             RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0},
-            delta[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+            delta[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}, TurbViscosity;
   su2double AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
 
@@ -1992,6 +1992,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   bool energy = config->GetEnergy_Equation();
   bool QCR = config->GetQCR();
   bool axisymmetric = config->GetAxisymmetric();
+  bool roughwall = (config->GetnRoughWall() > 0);
 
   /// TODO: Move these ifs to specialized functions.
 
@@ -2096,6 +2097,14 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
         }
 
         Viscosity = nodes->GetLaminarViscosity(iPoint);
+        if (roughwall) {
+          unsigned short WallType; su2double Roughness_Height;
+          tie(WallType, Roughness_Height) = config->GetWallRoughnessProperties(Marker_Tag);
+          TurbViscosity = 0.0;
+          if (WallType == ROUGH) TurbViscosity = nodes->GetEddyViscosity(iPoint);
+          Viscosity = nodes->GetLaminarViscosity(iPoint);
+          Viscosity += TurbViscosity;
+        }
         Density = nodes->GetDensity(iPoint);
 
         Area = 0.0;
diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp
index cc5de7f335e4..0003ed4bf9e7 100644
--- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
index 5e6f4fa6810f..69826888890a 100644
--- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
+++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
@@ -48,6 +48,7 @@ CSourceBase_TurbSA::CSourceBase_TurbSA(unsigned short val_nDim,
   cb2   = 0.622;
   cb2_sigma = cb2/sigma;
   cw1 = cb1/k2+(1.0+cb2)/sigma;
+  cr1 = 0.5;
 
   /*--- Setup the Jacobian pointer, we need to return su2double** but
    *    we know the Jacobian is 1x1 so we use this "trick" to avoid
@@ -79,6 +80,9 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig
   su2double vmag, rey, re_theta, re_theta_t, re_v;
   su2double tu , nu_cr, nu_t, nu_BC, chi_1, chi_2, term1, term2, term_exponential;
 
+  // Set the boolean here depending on whether the point is closest to a rough wall or not.
+  roughwall = (roughness_i > 0.0);
+
   if (incompressible) {
     Density_i = V_i[nDim+2];
     Laminar_Viscosity_i = V_i[nDim+4];
@@ -120,11 +124,20 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig
 
     dist_i_2 = dist_i*dist_i;
     nu = Laminar_Viscosity_i/Density_i;
-    Ji = TurbVar_i[0]/nu;
+
+    /*--- Modified values for roughness ---*/
+    /*--- Ref: Aupoix, B. and Spalart, P. R., "Extensions of the Spalart-Allmaras Turbulence Model to Account for Wall Roughness,"
+     * International Journal of Heat and Fluid Flow, Vol. 24, 2003, pp. 454-462. ---*/
+    /* --- See https://turbmodels.larc.nasa.gov/spalart.html#sarough for detailed explanation. ---*/
+
+    Ji = TurbVar_i[0]/nu  + cr1*(roughness_i/(dist_i+EPS)); //roughness_i = 0 for smooth walls and Ji remains the same, changes only if roughness is specified.
     Ji_2 = Ji*Ji;
     Ji_3 = Ji_2*Ji;
     fv1 = Ji_3/(Ji_3+cv1_3);
-    fv2 = 1.0 - Ji/(1.0+Ji*fv1);
+
+    /*--- Using a modified relation so as to not change the Shat that depends on fv2. ---*/
+    fv2 = 1.0 - TurbVar_i[0]/(nu+TurbVar_i[0]*fv1);   // From NASA turb modeling resource and 2003 paper
+
     ft2 = ct3*exp(-ct4*Ji_2);
     S = Omega;
     inv_k2_d2 = 1.0/(k2*dist_i_2);
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 050867370b9a..550af3f0fb4f 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -300,7 +300,7 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
 
 void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) {
 
-  const su2double cv1_3 = 7.1*7.1*7.1;
+  const su2double cv1_3 = 7.1*7.1*7.1, cR1 = 0.5, rough_const = 0.03;
 
   const bool neg_spalart_allmaras = (config->GetKind_Turb_Model() == SA_NEG);
 
@@ -314,8 +314,15 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain
 
     su2double nu  = mu/rho;
     su2double nu_hat = nodes->GetSolution(iPoint,0);
+    su2double roughness = geometry->nodes->GetRoughnessHeight(iPoint);
+    su2double dist = geometry->nodes->GetWall_Distance(iPoint);
+
+    dist += rough_const*roughness;
+
+    su2double Ji   = nu_hat/nu ;
+    if (roughness > 1.0e-10)
+      Ji+= cR1*roughness/(dist+EPS);
 
-    su2double Ji   = nu_hat/nu;
     su2double Ji_3 = Ji*Ji*Ji;
     su2double fv1  = Ji_3/(Ji_3+cv1_3);
 
@@ -329,6 +336,7 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain
 
 }
 
+
 void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,
                                     CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
 
@@ -338,6 +346,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
 
   CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
 
+
   /*--- Pick one numerics object per thread. ---*/
   CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS];
 
@@ -379,9 +388,22 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
 
     if (config->GetKind_HybridRANSLES() == NO_HYBRIDRANSLES) {
 
+    /*--- For the SA model, wall roughness is accounted by modifying the computed wall distance
+       *                              d_new = d + 0.03 k_s
+       *    where k_s is the equivalent sand grain roughness height that is specified in cfg file.
+       *    For smooth walls, wall roughness is zero and computed wall distance remains the same. */
+
+      su2double modifiedWallDistance = geometry->nodes->GetWall_Distance(iPoint);
+
+      modifiedWallDistance += 0.03*geometry->nodes->GetRoughnessHeight(iPoint);
+
       /*--- Set distance to the surface ---*/
 
-      numerics->SetDistance(geometry->nodes->GetWall_Distance(iPoint), 0.0);
+      numerics->SetDistance(modifiedWallDistance, 0.0);
+
+      /*--- Set the roughness of the closest wall. ---*/
+
+      numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), 0.0 );
 
     } else {
 
@@ -442,10 +464,16 @@ void CTurbSASolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_conta
     SU2_OMP_BARRIER
     return;
   }
+  
+  bool rough_wall = false;
+  string Marker_Tag = config->GetMarker_All_TagBound(val_marker); 
+  unsigned short WallType; su2double Roughness_Height;
+  tie(WallType, Roughness_Height) = config->GetWallRoughnessProperties(Marker_Tag);
+  if (WallType == ROUGH ) rough_wall = true;
 
   /*--- The dirichlet condition is used only without wall function, otherwise the
    convergence is compromised as we are providing nu tilde values for the
-   first point of the wall. ---*/
+   first point of the wall  ---*/
 
   SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
   for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
@@ -455,18 +483,49 @@ void CTurbSASolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_conta
     /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
 
     if (geometry->nodes->GetDomain(iPoint)) {
-
-      for (auto iVar = 0u; iVar < nVar; iVar++)
-        nodes->SetSolution_Old(iPoint,iVar,0.0);
-
-      LinSysRes.SetBlock_Zero(iPoint);
-
-      /*--- Includes 1 in the diagonal ---*/
-
-      Jacobian.DeleteValsRowi(iPoint);
+      if (!rough_wall) {
+        for (auto iVar = 0u; iVar < nVar; iVar++)
+          nodes->SetSolution_Old(iPoint,iVar,0.0);
+
+        LinSysRes.SetBlock_Zero(iPoint);
+
+        /*--- Includes 1 in the diagonal ---*/
+
+        Jacobian.DeleteValsRowi(iPoint);
+       } else {
+         /*--- For rough walls, the boundary condition is given by
+          * (\frac{\partial \nu}{\partial n})_wall = \frac{\nu}{0.03*k_s}
+          * where \nu is the solution variable, $n$ is the wall normal direction
+          * and k_s is the equivalent sand grain roughness specified. ---*/
+
+         /*--- Compute dual-grid area and boundary normal ---*/
+         su2double Normal[MAXNDIM] = {0.0};
+         for (auto iDim = 0u; iDim < nDim; iDim++)
+           Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim);
+          
+         su2double Area = GeometryToolbox::Norm(nDim, Normal);
+
+         /*--- Get laminar_viscosity and density ---*/
+         su2double sigma = 2.0/3.0;
+         su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
+         su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
+
+         su2double nu_total = (laminar_viscosity/density + nodes->GetSolution(iPoint,0));
+
+         su2double coeff = (nu_total/sigma);
+         su2double RoughWallBC = nodes->GetSolution(iPoint,0)/(0.03*Roughness_Height);
+
+         su2double Res_Wall;// = new su2double [nVar];
+         Res_Wall = coeff*RoughWallBC*Area;
+         LinSysRes.SubtractBlock(iPoint, &Res_Wall);
+
+         su2double Jacobian_i = (laminar_viscosity*Area)/(0.03*Roughness_Height*sigma);
+         Jacobian_i += 2.0*RoughWallBC*Area/sigma;
+         Jacobian_i = -Jacobian_i;
+         Jacobian.AddVal2Diag(iPoint, Jacobian_i);
+      }
     }
   }
-
 }
 
 void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
index b17bdef9b1ef..f83637f20b2a 100644
--- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp
@@ -390,43 +390,90 @@ void CTurbSSTSolver::Source_Template(CGeometry *geometry, CSolver **solver_conta
 
 void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
                                       CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
+  bool rough_wall = false;
+  string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+  unsigned short WallType; su2double Roughness_Height;
+  tie(WallType, Roughness_Height) = config->GetWallRoughnessProperties(Marker_Tag);
+  if (WallType == ROUGH ) rough_wall = true;
 
   SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
   for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
 
+
     const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
     /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
     if (geometry->nodes->GetDomain(iPoint)) {
 
-      /*--- distance to closest neighbor ---*/
-      const auto jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
 
-      su2double distance2 = GeometryToolbox::SquaredDistance(nDim,
+      if (rough_wall) {
+
+        /*--- Set wall values ---*/
+        su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
+        su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
+
+        su2double WallShearStress = 0.0;
+        for (auto iDim = 0; iDim < nDim; iDim++)
+          WallShearStress += pow(solver_container[FLOW_SOL]->GetCSkinFriction(val_marker, iVertex, iDim),2.0);
+
+        WallShearStress = sqrt(WallShearStress);
+        /*--- Compute non-dimensional velocity ---*/
+        su2double FrictionVel = sqrt(fabs(WallShearStress)/density);
+
+        /*--- Compute roughness in wall units. ---*/
+        //su2double Roughness_Height = config->GetWall_RoughnessHeight(Marker_Tag);
+        su2double kPlus = FrictionVel*Roughness_Height*density/laminar_viscosity;
+
+        su2double S_R= 0.0;
+        /*--- Reference 1 original Wilcox (1998) ---*/
+        /*if (kPlus <= 25)
+            S_R = (50/(kPlus+EPS))*(50/(kPlus+EPS));
+          else
+            S_R = 100/(kPlus+EPS);*/
+
+        /*--- Reference 2 from D.C. Wilcox Turbulence Modeling for CFD (2006) ---*/
+        if (kPlus <= 5)
+          S_R = (200/(kPlus+EPS))*(200/(kPlus+EPS));
+        else
+          S_R = 100/(kPlus+EPS) + ((200/(kPlus+EPS))*(200/(kPlus+EPS)) - 100/(kPlus+EPS))*exp(5-kPlus);
+
+        /*--- Modify the omega to account for a rough wall. ---*/
+        su2double solution[2];
+        solution[0] = 0.0;
+        solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);
+
+        /*--- Set the solution values and zero the residual ---*/
+        nodes->SetSolution_Old(iPoint,solution);
+        nodes->SetSolution(iPoint,solution);
+        LinSysRes.SetBlock_Zero(iPoint);
+      } else {
+        /*--- distance to closest neighbor ---*/
+        const auto jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+
+        su2double distance2 = GeometryToolbox::SquaredDistance(nDim,
                                                              geometry->nodes->GetCoord(iPoint),
                                                              geometry->nodes->GetCoord(jPoint));
-      /*--- Set wall values ---*/
+        /*--- Set wall values ---*/
 
-      su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint);
-      su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint);
+        su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint);
+        su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint);
 
-      su2double beta_1 = constants[4];
+        su2double beta_1 = constants[4];
+        su2double solution[2];
+        solution[0] = 0.0;
+        solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance2);
 
-      su2double solution[2];
-      solution[0] = 0.0;
-      solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance2);
-
-      /*--- Set the solution values and zero the residual ---*/
-      nodes->SetSolution_Old(iPoint,solution);
-      nodes->SetSolution(iPoint,solution);
-      LinSysRes.SetBlock_Zero(iPoint);
+        /*--- Set the solution values and zero the residual ---*/
+        nodes->SetSolution_Old(iPoint,solution);
+        nodes->SetSolution(iPoint,solution);
+        LinSysRes.SetBlock_Zero(iPoint);
+      }
 
       /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/
       Jacobian.DeleteValsRowi(iPoint*nVar);
       Jacobian.DeleteValsRowi(iPoint*nVar+1);
     }
   }
-
 }
 
 void CTurbSSTSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp
index 4c418f48e33a..179d659d1cfc 100644
--- a/SU2_CFD/src/solvers/CTurbSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSolver.cpp
@@ -275,6 +275,11 @@ void CTurbSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSo
     numerics->SetF1blending(nodes->GetF1blending(iPoint),
                             nodes->GetF1blending(jPoint));
 
+  /*--- Roughness heights. ---*/
+    if (config->GetKind_Turb_Model() == SA)
+      numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint),
+                             geometry->nodes->GetRoughnessHeight(jPoint));
+
   /*--- Compute residual, and Jacobians ---*/
 
   auto residual = numerics->ComputeResidual(config);
diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp
index da1e836d282a..7ea2c0c5b92a 100644
--- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp
index 0889748e5210..2f6ce8a3e3bd 100644
--- a/SU2_CFD/src/variables/CTurbVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbVariable.cpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -45,7 +45,7 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
   if (config->GetReconstructionGradientRequired()) {
     Gradient_Aux.resize(nPoint,nVar,nDim,0.0);
   }
-  
+
   if (config->GetLeastSquaresRequired()) {
     Rmatrix.resize(nPoint,nDim,nDim,0.0);
   }
@@ -62,5 +62,5 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
   /* Under-relaxation parameter. */
   UnderRelaxation.resize(nPoint) = su2double(1.0);
   LocalCFL.resize(nPoint) = su2double(0.0);
-  
+
 }
diff --git a/TestCases/incomp_rans/rough_flatplate/rough_flatplate_incomp.cfg b/TestCases/incomp_rans/rough_flatplate/rough_flatplate_incomp.cfg
new file mode 100644
index 000000000000..0142f508d274
--- /dev/null
+++ b/TestCases/incomp_rans/rough_flatplate/rough_flatplate_incomp.cfg
@@ -0,0 +1,291 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Turbulent flow over rough flat plate with zero             %
+%                   pressure gradient                                          %
+% Author: Akshay Koodly                                                        %
+% Date: 2020.07.07                                                             %
+% File Version 7.0.x "Blackbird"                                               %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+% Physical governing equations (EULER, NAVIER_STOKES,
+%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
+%                               POISSON_EQUATION)
+SOLVER= INC_RANS
+%
+% If Navier-Stokes, kind of turbulent model (NONE, SA)
+KIND_TURB_MODEL= SA
+
+% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
+MATH_PROBLEM= DIRECT
+%
+% Restart solution (NO, YES)
+RESTART_SOL= NO
+%
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+% Density model within the incompressible flow solver.
+% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE,
+% an appropriate fluid model must be selected.
+INC_DENSITY_MODEL= CONSTANT
+%
+% Solve the energy equation in the incompressible flow solver
+INC_ENERGY_EQUATION = NO
+%
+% Initial density for incompressible flows (1.2886 kg/m^3 by default)
+INC_DENSITY_INIT= 1.32905
+%
+% Initial velocity for incompressible flows (1.0,0,0 m/s by default)
+INC_VELOCITY_INIT= ( 69.4448, 0.0, 0.0 )
+%
+% Initial temperature for incompressible flows that include the
+% energy equation (288.15 K by default). Value is ignored if
+% INC_ENERGY_EQUATION is false.
+INC_TEMPERATURE_INIT= 300.0
+%
+% Non-dimensionalization scheme for incompressible flows. Options are
+% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL.
+% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen.
+INC_NONDIM= INITIAL_VALUES
+%
+% Reference density for incompressible flows (1.0 kg/m^3 by default)
+INC_DENSITY_REF= 1.0
+%
+% Reference velocity for incompressible flows (1.0 m/s by default)
+INC_VELOCITY_REF= 1.0
+%
+% Reference temperature for incompressible flows that include the
+% energy equation (1.0 K by default)
+INC_TEMPERATURE_REF = 1.0
+%
+% List of inlet types for incompressible flows. List length must
+% match number of inlet markers. Options: VELOCITY_INLET, PRESSURE_INLET.
+INC_INLET_TYPE= VELOCITY_INLET
+
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY).
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+%
+% Molecular Viscosity that would be constant (1.716E-5 by default)
+MU_CONSTANT= 1.84592e-05
+%
+% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI)
+MU_REF= 1.716E-5
+%
+% Sutherland Temperature Ref (273.15 K default value for AIR SI)
+MU_T_REF= 273.15
+%
+% Sutherland constant (110.4 default value for AIR SI)
+SUTHERLAND_CONSTANT= 110.4
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+% Reference origin for moment computation
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+%
+% Reference length for pitching, rolling, and yawing non-dimensional moment
+REF_LENGTH= 1.0
+%
+% Reference area for force coefficients (0 implies automatic calculation)
+REF_AREA= 2.0
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Navier-Stokes wall boundary marker(s) (NONE = no marker)
+MARKER_HEATFLUX= ( wall, 0.0 )
+%WALL_ROUGHNESS = (wall, 0.000061)
+%WALL_ROUGHNESS = (wall, 0.000123)
+WALL_ROUGHNESS = (wall, 0.000246)
+%WALL_ROUGHNESS = (wall, 0.000984)
+%
+%
+% Inlet boundary marker(s) (NONE = no marker) 
+% Format: ( inlet marker, total temperature, total pressure, flow_direction_x, 
+%           flow_direction_y, flow_direction_z, ... )
+MARKER_INLET= ( inlet, 300.0, 69.4448, 1.0, 0.0, 0.0 )
+%
+% Outlet boundary marker(s) (NONE = no marker)
+% Format: ( outlet marker, back pressure, ... )
+MARKER_OUTLET= ( outlet, 0.0, farfield, 0.0 )
+%
+INC_OUTLET_TYPE= PRESSURE_OUTLET,PRESSURE_OUTLET
+%
+% Symmetry boundary marker(s) (NONE = no marker)
+MARKER_SYM= ( symmetry )
+%
+% Marker(s) of the surface to be plotted or designed
+MARKER_PLOTTING= ( wall )
+%
+% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated
+MARKER_MONITORING= ( wall )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+% Numerical method for spatial gradients (GREEN_GAUSS, LEAST_SQUARES, 
+%                                         WEIGHTED_LEAST_SQUARES)
+NUM_METHOD_GRAD= GREEN_GAUSS
+%
+% Courant-Friedrichs-Lewy condition of the finest grid
+CFL_NUMBER= 100.0
+%
+% Adaptive CFL number (NO, YES)
+CFL_ADAPT= NO
+%
+% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
+%                                        CFL max value )
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.1, 100.0 )
+%
+% Runge-Kutta alpha coefficients
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+%
+% Number of total iterations
+ITER= 1000
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+% Linear solver for implicit formulations (BCGSTAB, FGMRES)
+LINEAR_SOLVER= FGMRES
+%
+% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS)
+LINEAR_SOLVER_PREC= ILU
+%
+% Linael solver ILU preconditioner fill-in level (0 by default)
+LINEAR_SOLVER_ILU_FILL_IN= 0
+%
+% Minimum error of the linear solver for implicit formulations
+LINEAR_SOLVER_ERROR= 1E-12
+%
+% Max number of iterations of the linear solver for the implicit formulation
+LINEAR_SOLVER_ITER= 20
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+%
+% Coefficient for the limiter
+VENKAT_LIMITER_COEFF= 0.1
+%
+% Coefficient for the sharp edges limiter
+ADJ_SHARP_LIMITER_COEFF= 3.0
+%
+% Reference coefficient (sensitivity) for detecting sharp edges.
+REF_SHARP_EDGES= 3.0
+%
+% Remove sharp edges from the sensitivity evaluation (NO, YES)
+SENS_REMOVE_SHARP= NO
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+% Multi-Grid Levels (0 = no multi-grid)
+MGLEVEL= 0
+%
+% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
+MGCYCLE= V_CYCLE
+%
+% Multi-grid pre-smoothing level
+MG_PRE_SMOOTH= ( 1, 1, 1, 1 )
+%
+% Multi-grid post-smoothing level
+MG_POST_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Jacobi implicit smoothing of the correction
+MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Damping factor for the residual restriction
+MG_DAMP_RESTRICTION= 0.8
+%
+% Damping factor for the correction prolongation
+MG_DAMP_PROLONGATION= 0.8
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
+%                              TURKEL_PREC, MSW)
+CONV_NUM_METHOD_FLOW= FDS
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_FLOW= YES
+%
+% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG,
+%                BARTH_JESPERSEN, VAN_ALBADA_EDGE)
+SLOPE_LIMITER_FLOW= NONE
+%
+% 2nd and 4th order artificial dissipation coefficients
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+%
+% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+% Convective numerical method (SCALAR_UPWIND)
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_TURB= NO
+%
+% Slope limiter (VENKATAKRISHNAN, MINMOD)
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+%
+% Time discretization (EULER_IMPLICIT)
+TIME_DISCRE_TURB= EULER_IMPLICIT
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence criteria (CAUCHY, RESIDUAL)
+%
+CONV_FIELD= RMS_VELOCITY-X
+%
+% Min value of the residual (log10 of the residual)
+CONV_RESIDUAL_MINVAL= -14
+%
+% Start convergence criteria at iteration number
+CONV_STARTITER= 10
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+% Mesh input file
+%MESH_FILENAME= mesh_flatplate_turb_69x49.su2
+MESH_FILENAME= mesh_flatplate_turb_137x97.su2
+%MESH_FILENAME= mesh_flatplate_turb_273x193.su2
+%
+% Mesh input file format (SU2, CGNS, NETCDF_ASCII)
+MESH_FORMAT= SU2
+%
+% Mesh output file
+MESH_OUT_FILENAME= mesh_out.su2
+%
+% Restart flow input file
+SOLUTION_FILENAME= restart_flow
+%
+% Output file format (PARAVIEW, TECPLOT, SLT)
+TABULAR_FORMAT= CSV
+
+OUTPUT_FILES= PARAVIEW, RESTART, SURFACE_PARAVIEW
+%
+% Output file convergence history (w/o extension) 
+CONV_FILENAME= history
+%
+% Output file restart flow
+RESTART_FILENAME= restart_flow
+%
+% Output file flow (w/o extension) variables
+VOLUME_FILENAME= flow
+%
+% Output file surface flow coefficient (w/o extension)
+SURFACE_FILENAME= surface_flow
+%
+% Writing solution file frequency
+OUTPUT_WRT_FREQ= 100
+%
+% Writing convergence history frequency
+WRT_CON_FREQ= 1
+%
+SCREEN_OUTPUT= (WALL_TIME,INNER_ITER, RMS_VELOCITY-X, RMS_NU_TILDE, LIFT,DRAG)
+%
+WRT_FORCES_BREAKDOWN= YES
diff --git a/config_template.cfg b/config_template.cfg
index 14f8749e6322..47825d0caede 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -812,7 +812,11 @@ GILES_EXTRA_RELAXFACTOR= ( 0.05, 0.05)
 %
 % YES Non reflectivity activated, NO the Giles BC behaves as a normal 1D characteristic-based BC
 SPATIAL_FOURIER= NO
-
+% ------------------------ WALL ROUGHNESS DEFINITION --------------------------%
+% The equivalent sand grain roughness height (k_s) on each of the wall. This must be in m. 
+% This is a list of (string, double) each element corresponding to the MARKER defined in WALL_TYPE.
+%WALL_ROUGHNESS = (wall1, ks1, wall2, ks2)
+%WALL_ROUGHNESS = (wall1, ks1, wall2, 0.0) %is also allowed
 % ------------------------ SURFACES IDENTIFICATION ----------------------------%
 %
 % Marker(s) of the surface in the surface flow solution file
@@ -843,6 +847,11 @@ MARKER_ANALYZE_AVERAGE = MASSFLUX
 %
 % Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
 NUM_METHOD_GRAD= GREEN_GAUSS
+
+% Numerical method for spatial gradients to be used for MUSCL reconstruction
+% Options are (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES, LEAST_SQUARES). Default value is
+% NONE and the method specified in NUM_METHOD_GRAD is used. 
+NUM_METHOD_GRAD_RECON = LEAST_SQUARES
 %
 % CFL number (initial value for the adaptive CFL number)
 CFL_NUMBER= 15.0