diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 433f019efe3..5fcc39e7798 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -699,10 +699,12 @@ class CConfig {
   *Marker_WallFunctions,              /*!< \brief Markers for which wall functions must be applied. */
   *Marker_SobolevBC;                  /*!< \brief Markers in the gradient solver */
 
-  unsigned short  nConfig_Files;          /*!< \brief Number of config files for multiphysics problems. */
-  string *Config_Filenames;               /*!< \brief List of names for configuration files. */
-  SST_OPTIONS *SST_Options;               /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/
-  unsigned short nSST_Options;            /*!< \brief number of SST options specified. */
+  unsigned short nConfig_Files;       /*!< \brief Number of config files for multiphysics problems. */
+  string *Config_Filenames;           /*!< \brief List of names for configuration files. */
+  SST_OPTIONS *SST_Options;           /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/
+  SA_OPTIONS *SA_Options;             /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/
+  unsigned short nSST_Options;        /*!< \brief Number of SST options specified. */
+  unsigned short nSA_Options;         /*!< \brief Number of SA options specified. */
   WALL_FUNCTIONS  *Kind_WallFunctions;        /*!< \brief The kind of wall function to use for the corresponding markers. */
   unsigned short  **IntInfo_WallFunctions;    /*!< \brief Additional integer information for the wall function markers. */
   su2double       **DoubleInfo_WallFunctions; /*!< \brief Additional double information for the wall function markers. */
@@ -1015,7 +1017,6 @@ class CConfig {
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
-  bool QCR;                    /*!< \brief Spalart-Allmaras with Quadratic Constitutive Relation, 2000 version (SA-QCR2000) . */
 
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
   unsigned short nSpanMaxAllZones;  /*!< \brief number of maximum span-wise sections for all zones */
@@ -1137,7 +1138,8 @@ class CConfig {
   unsigned short nScreenOutput,   /*!< \brief Number of screen output variables (max: 6). */
   nHistoryOutput, nVolumeOutput;  /*!< \brief Number of variables printed to the history file. */
   bool Multizone_Residual;        /*!< \brief Determines if memory should be allocated for the multizone residual. */
-  SST_ParsedOptions sstParsedOptions; /*!< \brief additional parameters for the SST turbulence model */
+  SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */
+  SA_ParsedOptions saParsedOptions;   /*!< \brief Additional parameters for the SA turbulence model. */
   su2double uq_delta_b;         /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */
   unsigned short eig_val_comp;  /*!< \brief Parameter used to determine type of eigenvalue perturbation */
   su2double uq_urlx;            /*!< \brief Under-relaxation factor */
@@ -9058,11 +9060,6 @@ class CConfig {
    */
   su2double GetConst_DES(void) const { return Const_DES; }
 
-  /*!
-   * \brief Get QCR (SA-QCR2000).
-   */
-  bool GetQCR(void) const { return QCR;}
-
   /*!
    * \brief Get if AD preaccumulation should be performed.
    */
@@ -9629,4 +9626,10 @@ class CConfig {
    */
   SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; }
 
+  /*!
+   * \brief Get parsed SA option data structure.
+   * \return SA option data structure.
+   */
+  SA_ParsedOptions GetSAParsedOptions() const { return saParsedOptions; }
+
 };
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index afb1fe643d4..8b09a93c737 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -912,19 +912,11 @@ static const MapType<std::string, LIMITER> Limiter_Map = {
 enum class TURB_MODEL {
   NONE,      /*!< \brief No turbulence model. */
   SA,        /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
-  SA_NEG,    /*!< \brief Kind of Turbulent model (Negative Spalart-Allmaras). */
-  SA_E,      /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards). */
-  SA_COMP,   /*!< \brief Kind of Turbulent model (Spalart-Allmaras Compressibility Correction). */
-  SA_E_COMP, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards with Compressibility Correction). */
   SST,       /*!< \brief Kind of Turbulence model (Menter SST). */
 };
 static const MapType<std::string, TURB_MODEL> Turb_Model_Map = {
   MakePair("NONE", TURB_MODEL::NONE)
   MakePair("SA", TURB_MODEL::SA)
-  MakePair("SA_NEG", TURB_MODEL::SA_NEG)
-  MakePair("SA_E", TURB_MODEL::SA_E)
-  MakePair("SA_COMP", TURB_MODEL::SA_COMP)
-  MakePair("SA_E_COMP", TURB_MODEL::SA_E_COMP)
   MakePair("SST", TURB_MODEL::SST)
 };
 
@@ -944,10 +936,6 @@ inline TURB_FAMILY TurbModelFamily(TURB_MODEL model) {
     case TURB_MODEL::NONE:
       return TURB_FAMILY::NONE;
     case TURB_MODEL::SA:
-    case TURB_MODEL::SA_NEG:
-    case TURB_MODEL::SA_E:
-    case TURB_MODEL::SA_COMP:
-    case TURB_MODEL::SA_E_COMP:
       return TURB_FAMILY::SA;
     case TURB_MODEL::SST:
       return TURB_FAMILY::KW;
@@ -998,6 +986,7 @@ struct SST_ParsedOptions {
  * \param[in] SST_Options - Selected SST option from config.
  * \param[in] nSST_Options - Number of options selected.
  * \param[in] rank - MPI rank.
+ * \return Struct with SST options.
  */
 inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigned short nSST_Options, int rank) {
   SST_ParsedOptions SSTParsedOptions;
@@ -1056,18 +1045,115 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   return SSTParsedOptions;
 }
 
+/*!
+ * \brief SA Options
+ */
+enum class SA_OPTIONS {
+  NONE,     /*!< \brief No option / default. */
+  NEG,      /*!< \brief Negative SA. */
+  EDW,      /*!< \brief Edwards version. */
+  FT2,      /*!< \brief Use FT2 term. */
+  QCR2000,  /*!< \brief Quadratic constitutive relation. */
+  COMP,     /*!< \brief Compressibility correction. */
+  ROT,      /*!< \brief Rotation correction. */
+  BC,       /*!< \brief Bas-Cakmakcioclu transition. */
+  EXP,      /*!< \brief Allow experimental combinations of options (according to TMR). */
+};
+static const MapType<std::string, SA_OPTIONS> SA_Options_Map = {
+  MakePair("NONE", SA_OPTIONS::NONE)
+  MakePair("NEGATIVE", SA_OPTIONS::NEG)
+  MakePair("EDWARDS", SA_OPTIONS::EDW)
+  MakePair("WITHFT2", SA_OPTIONS::FT2)
+  MakePair("QCR2000", SA_OPTIONS::QCR2000)
+  MakePair("COMPRESSIBILITY", SA_OPTIONS::COMP)
+  MakePair("ROTATION", SA_OPTIONS::ROT)
+  MakePair("BCM", SA_OPTIONS::BC)
+  MakePair("EXPERIMENTAL", SA_OPTIONS::EXP)
+};
+
+/*!
+ * \brief Structure containing parsed SA options.
+ */
+struct SA_ParsedOptions {
+  SA_OPTIONS version = SA_OPTIONS::NONE;  /*!< \brief SA base model. */
+  bool ft2 = false;                       /*!< \brief Use ft2 term. */
+  bool qcr2000 = false;                   /*!< \brief Use QCR-2000. */
+  bool comp = false;                      /*!< \brief Use compressibility correction. */
+  bool rot = false;                       /*!< \brief Use rotation correction. */
+  bool bc = false;                        /*!< \brief BC transition. */
+};
+
+/*!
+ * \brief Function to parse SA options.
+ * \param[in] SA_Options - Selected SA option from config.
+ * \param[in] nSA_Options - Number of options selected.
+ * \param[in] rank - MPI rank.
+ * \return Struct with SA options.
+ */
+inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned short nSA_Options, int rank) {
+  SA_ParsedOptions SAParsedOptions;
+
+  auto IsPresent = [&](SA_OPTIONS option) {
+    const auto sa_options_end = SA_Options + nSA_Options;
+    return std::find(SA_Options, sa_options_end, option) != sa_options_end;
+  };
+
+  const bool found_neg = IsPresent(SA_OPTIONS::NEG);
+  const bool found_edw = IsPresent(SA_OPTIONS::EDW);
+  const bool found_bsl = !found_neg && !found_edw;
+
+  if (found_neg && found_edw) {
+    SU2_MPI::Error("Two versions (Negative and Edwards) selected for SA_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
+  }
+
+  if (found_bsl) {
+    SAParsedOptions.version = SA_OPTIONS::NONE;
+  } else if (found_neg) {
+    SAParsedOptions.version = SA_OPTIONS::NEG;
+  } else {
+    SAParsedOptions.version = SA_OPTIONS::EDW;
+  }
+  SAParsedOptions.ft2 = IsPresent(SA_OPTIONS::FT2);
+  SAParsedOptions.qcr2000 = IsPresent(SA_OPTIONS::QCR2000);
+  SAParsedOptions.comp = IsPresent(SA_OPTIONS::COMP);
+  SAParsedOptions.rot = IsPresent(SA_OPTIONS::ROT);
+  SAParsedOptions.bc = IsPresent(SA_OPTIONS::BC);
+
+  /*--- Validate user settings when not in experimental mode. ---*/
+  if (!IsPresent(SA_OPTIONS::EXP)) {
+    const bool any_but_bc = SAParsedOptions.ft2 || SAParsedOptions.qcr2000 || SAParsedOptions.comp || SAParsedOptions.rot;
+
+    switch (SAParsedOptions.version) {
+      case SA_OPTIONS::NEG:
+        if (!SAParsedOptions.ft2 || SAParsedOptions.bc)
+          SU2_MPI::Error("A non-standard version of SA-neg was requested (see https://turbmodels.larc.nasa.gov/spalart.html).\n"
+                         "If you want to continue, add EXPERIMENTAL to SA_OPTIONS.", CURRENT_FUNCTION);
+        break;
+      case SA_OPTIONS::EDW:
+        if (any_but_bc || SAParsedOptions.bc)
+          SU2_MPI::Error("A non-standard version of SA-noft2-Edwards was requested (see https://turbmodels.larc.nasa.gov/spalart.html).\n"
+                         "If you want to continue, add EXPERIMENTAL to SA_OPTIONS.", CURRENT_FUNCTION);
+        break;
+      default:
+        if (SAParsedOptions.bc && any_but_bc)
+          SU2_MPI::Error("A non-standard version of SA-BCM was requested (see https://turbmodels.larc.nasa.gov/spalart.html).\n"
+                         "If you want to continue, add EXPERIMENTAL to SA_OPTIONS.", CURRENT_FUNCTION);
+        break;
+    }
+  }
+  return SAParsedOptions;
+}
+
 /*!
  * \brief Types of transition models
  */
 enum class TURB_TRANS_MODEL {
   NONE,  /*!< \brief No transition model. */
   LM,    /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */
-  BC     /*!< \brief Kind of transition model (BAS-CAKMAKCIOGLU (BC) for Spalart-Allmaras). */
 };
 static const MapType<std::string, TURB_TRANS_MODEL> Trans_Model_Map = {
   MakePair("NONE", TURB_TRANS_MODEL::NONE)
   MakePair("LM", TURB_TRANS_MODEL::LM)
-  MakePair("BC", TURB_TRANS_MODEL::BC)
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index b03edf579d4..e6699868e5e 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1091,6 +1091,8 @@ void CConfig::SetConfig_Options() {
   addEnumOption("KIND_TURB_MODEL", Kind_Turb_Model, Turb_Model_Map, TURB_MODEL::NONE);
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SST turbulence model options/corrections. \n Options: see \link SST_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumListOption("SST_OPTIONS", nSST_Options, SST_Options, SST_Options_Map);
+  /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
+  addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map);
 
   /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE);
@@ -2767,9 +2769,6 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
   addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS);
 
-  /* DESCRIPTION: Activate SA Quadratic Constitutive Relation, 2000 version */
-  addBoolOption("SA_QCR", QCR, false);
-
   /* DESCRIPTION: Compute Average for unsteady simulations */
   addBoolOption("COMPUTE_AVERAGE", Compute_Average, false);
 
@@ -2996,6 +2995,8 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){
             newString.append("SOLID_DENSITY is deprecated. Use MATERIAL_DENSITY instead.\n\n");
           else if (!option_name.compare("SOLID_TEMPERATURE_INIT"))
             newString.append("SOLID_TEMPERATURE_INIT is deprecated. Use FREESTREAM_TEMPERATURE instead.\n\n");
+          else if (!option_name.compare("SA_QCR"))
+            newString.append("SA_QCR is deprecated. Use SA_OPTIONS=QCR2000 instead.\n\n");
           else {
             /*--- Find the most likely candidate for the unrecognized option, based on the length
              of start and end character sequences shared by candidates and the option. ---*/
@@ -3054,9 +3055,21 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){
       string out = option_map[option_name]->SetValue(option_value);
       if (out.compare("") != 0) {
         /*--- valid option, but deprecated value ---*/
-        if ((!option_name.compare("KIND_TURB_MODEL")) && (option_value[0]=="SST_SUST"))
-          errorString.append("Option KIND_TURB_MODEL=SST_SUST is deprecated. Use KIND_TURB_MODEL=SST, SST_OPTIONS=SUSTAINING instead.\n");
-
+        if (!option_name.compare("KIND_TURB_MODEL")) {
+          if (option_value[0] == "SST_SUST")
+            errorString.append("Option KIND_TURB_MODEL=SST_SUST is deprecated. Use KIND_TURB_MODEL=SST, SST_OPTIONS=SUSTAINING instead.\n");
+          else if (option_value[0] == "SA_NEG")
+            errorString.append("Option KIND_TURB_MODEL=SA_NEG is deprecated. Use KIND_TURB_MODEL=SA, SA_OPTIONS=NEGATIVE instead.\n");
+          else if (option_value[0] == "SA_E")
+            errorString.append("Option KIND_TURB_MODEL=SA_E is deprecated. Use KIND_TURB_MODEL=SA, SA_OPTIONS=EDWARDS instead.\n");
+          else if (option_value[0] == "SA_COMP")
+            errorString.append("Option KIND_TURB_MODEL=SA_COMP is deprecated. Use KIND_TURB_MODEL=SA, SA_OPTIONS=COMPRESSIBILITY instead.\n");
+          else if (option_value[0] == "SA_E_COMP")
+            errorString.append("Option KIND_TURB_MODEL=SA_E_COMP is deprecated. Use KIND_TURB_MODEL=SA, SA_OPTIONS=EDWARDS,COMPRESSIBILITY instead.\n");
+        } else if (!option_name.compare("KIND_TRANS_MODEL")) {
+          if (option_value[0] == "BC")
+            errorString.append("Option KIND_TRANS_MODEL=BC is deprecated. Use KIND_TURB_MODEL=SA, SA_OPTIONS=BCM instead.\n");
+        }
         errorString.append(out);
         errorString.append("\n");
         err_count++;
@@ -3399,8 +3412,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   }
 
   /*--- Postprocess SST_OPTIONS into structure. ---*/
-  if (Kind_Turb_Model==TURB_MODEL::SST) {
+  if (Kind_Turb_Model == TURB_MODEL::SST) {
     sstParsedOptions = ParseSSTOptions(SST_Options, nSST_Options, rank);
+  } else if (Kind_Turb_Model == TURB_MODEL::SA) {
+    saParsedOptions = ParseSAOptions(SA_Options, nSA_Options, rank);
   }
 
   /*--- Check if turbulence model can be used for AXISYMMETRIC case---*/
@@ -4588,10 +4603,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     for (int i=0; i<7; ++i) eng_cyl[i] /= 12.0;
   }
 
-  if ((Kind_Turb_Model != TURB_MODEL::SA) && (Kind_Trans_Model == TURB_TRANS_MODEL::BC)){
-    SU2_MPI::Error("BC transition model currently only available in combination with SA turbulence model!", CURRENT_FUNCTION);
-  }
-
   if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) {
     SU2_MPI::Error("The LM transition model is under maintenance.", CURRENT_FUNCTION);
   }
@@ -5861,11 +5872,26 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         cout << "Turbulence model: ";
         switch (Kind_Turb_Model) {
           case TURB_MODEL::NONE: break;
-          case TURB_MODEL::SA:        cout << "Spalart Allmaras" << endl; break;
-          case TURB_MODEL::SA_NEG:    cout << "Negative Spalart Allmaras" << endl; break;
-          case TURB_MODEL::SA_E:      cout << "Edwards Spalart Allmaras" << endl; break;
-          case TURB_MODEL::SA_COMP:   cout << "Compressibility Correction Spalart Allmaras" << endl; break;
-          case TURB_MODEL::SA_E_COMP: cout << "Compressibility Correction Edwards Spalart Allmaras" << endl; break;
+          case TURB_MODEL::SA:
+            switch (saParsedOptions.version) {
+              case SA_OPTIONS::NEG:
+                cout << "Negative-";
+                break;
+              case SA_OPTIONS::EDW:
+                cout << "Edwards-";
+                break;
+              default:
+                break;
+            }
+            cout << "Spalart-Allmaras";
+
+            if (!saParsedOptions.ft2) cout << "-noft2";
+            if (saParsedOptions.rot) cout << "-R";
+            if (saParsedOptions.comp) cout << "-comp";
+            if (saParsedOptions.qcr2000) cout << "-QCR2000";
+            if (saParsedOptions.bc) cout << "-BCM";
+            cout << endl;
+            break;
           case TURB_MODEL::SST:
             cout << "Menter's k-omega SST";
             if (sstParsedOptions.version == SST_OPTIONS::V1994) cout << "-1994";
@@ -5891,11 +5917,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
             cout << "." << endl;
             break;
         }
-        if (QCR) cout << "Using Quadratic Constitutive Relation, 2000 version (QCR2000)" << endl;
-        if (Kind_Trans_Model == TURB_TRANS_MODEL::BC) cout << "Using the revised BC transition model (2020)" << endl;
         cout << "Hybrid RANS/LES: ";
-        switch (Kind_HybridRANSLES){
-          case NO_HYBRIDRANSLES: cout <<  "No Hybrid RANS/LES" << endl; break;
+        switch (Kind_HybridRANSLES) {
+          case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
           case SA_DES:   cout << "Detached Eddy Simulation (DES97) " << endl; break;
           case SA_DDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
           case SA_ZDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 9fb9b298537..9be7239ba7c 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -71,8 +71,7 @@ class CSourceBase_TurbSA : public CNumerics {
   su2double Jacobian_Buffer; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */
 
   const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
-  const bool rotating_frame = false;
-  const bool transition = false;
+  const SA_ParsedOptions options; /*!< \brief Struct with SA options. */
 
  public:
   /*!
@@ -80,11 +79,10 @@ class CSourceBase_TurbSA : public CNumerics {
    * \param[in] nDim - Number of dimensions of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  CSourceBase_TurbSA(unsigned short nDim, unsigned short, const CConfig* config)
+  CSourceBase_TurbSA(unsigned short nDim, const CConfig* config)
       : CNumerics(nDim, 1, config),
         idx(nDim, config->GetnSpecies()),
-        rotating_frame(config->GetRotating_Frame()),
-        transition(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+        options(config->GetSAParsedOptions()) {
     /*--- Setup the Jacobian pointer, we need to return su2double** but we know
      * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/
     Jacobian_i = &Jacobian_Buffer;
@@ -121,9 +119,8 @@ class CSourceBase_TurbSA : public CNumerics {
 
     Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var);
 
-    /*--- Dacles-Mariani et. al. rotation correction ("-R"), this is applied by
-     * default for rotating frame, but should be controled in the config. ---*/
-    if (rotating_frame) {
+    /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/
+    if (options.rot) {
       var.Omega += 2.0 * min(0.0, StrainMag_i - var.Omega);
     }
 
@@ -173,7 +170,7 @@ class CSourceBase_TurbSA : public CNumerics {
 
       var.norm2_Grad = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0]);
 
-      if (transition) {
+      if (options.bc) {
         /*--- BC transition model (2020 revision). This should only be used with SA-noft2.
          * TODO: Consider making this part of the "SourceTerms" template. ---*/
         const su2double chi_1 = 0.002;
@@ -466,8 +463,8 @@ class CCompressibilityCorrection final : public ParentClass {
    * \param[in] nDim - Number of dimensions of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  CCompressibilityCorrection(unsigned short nDim, unsigned short, const CConfig* config)
-      : ParentClass(nDim, 0, config) {}
+  CCompressibilityCorrection(unsigned short nDim, const CConfig* config)
+      : ParentClass(nDim, config) {}
 
   /*!
    * \brief Residual for source term integration.
@@ -498,58 +495,67 @@ class CCompressibilityCorrection final : public ParentClass {
   }
 };
 
-}  // namespace detail
-
 /* =============================================================================
- * SPALART-ALLMARAS CLASSES
+ * HELPERS TO INSTANTIATE THE SA BASE CLASS
  * ============================================================================*/
 
-/// TODO: Factory method to create combinations of the different variations based on the config.
-/// See PR #1413, the combinations exposed should follow https://turbmodels.larc.nasa.gov/spalart.html
-
-/*!
- * \class CSourcePieceWise_TurbSA
- * \brief Class for integrating the source terms of the Spalart-Allmaras turbulence model equation.
- */
 template <class FlowIndices>
-using CSourcePieceWise_TurbSA = CSourceBase_TurbSA<FlowIndices, detail::Omega::Bsl, detail::ft2::Zero,
-                                                   detail::ModVort::Bsl, detail::r::Bsl, detail::SourceTerms::Bsl>;
+struct BaselineSA {
+  template <class Ft2>
+  using type = CSourceBase_TurbSA<FlowIndices, Omega::Bsl, Ft2, ModVort::Bsl, r::Bsl, SourceTerms::Bsl>;
+};
 
-/*!
- * \class CSourcePieceWise_TurbSA_COMP
- * \brief Class for integrating the source terms of the Spalart-Allmaras model with compressibility correction.
- */
 template <class FlowIndices>
-using CSourcePieceWise_TurbSA_COMP = detail::CCompressibilityCorrection<CSourcePieceWise_TurbSA<FlowIndices>>;
+struct NegativeSA {
+  template <class Ft2>
+  using type = CSourceBase_TurbSA<FlowIndices, Omega::Bsl, Ft2, ModVort::Neg, r::Bsl, SourceTerms::Neg>;
+};
 
-/*!
- * \class CSourcePieceWise_TurbSA_E
- * \brief Class for integrating the source terms of the Spalart-Allmaras model with Edwards modification.
- * \note From NASA Turbulence model site. http://turbmodels.larc.nasa.gov/spalart.html
- * This form was developed primarily to improve the near-wall numerical behavior of the model (i.e. the goal was to
- * improve the convergence behavior). The reference is: Edwards, J. R. and Chandra, S. "Comparison of Eddy
- * Viscosity-Transport Turbulence Models for Three-Dimensional, Shock-Separated Flowfields," AIAA Journal, Vol. 34, No.
- * 4, 1996, pp. 756-763.
- */
 template <class FlowIndices>
-using CSourcePieceWise_TurbSA_E = CSourceBase_TurbSA<FlowIndices, detail::Omega::Edw, detail::ft2::Zero,
-                                                     detail::ModVort::Edw, detail::r::Edw, detail::SourceTerms::Bsl>;
+struct EdwardsSA {
+  template <class Ft2>
+  using type = CSourceBase_TurbSA<FlowIndices, Omega::Edw, Ft2, ModVort::Edw, r::Edw, SourceTerms::Bsl>;
+};
 
-/*!
- * \class CSourcePieceWise_TurbSA_E_COMP
- * \brief Class for integrating the source terms of the Spalart-Allmaras model with Edwards modification
- *        and compressibility correction.
- */
-template <class FlowIndices>
-using CSourcePieceWise_TurbSA_E_COMP = detail::CCompressibilityCorrection<CSourcePieceWise_TurbSA_E<FlowIndices>>;
+template <class BaseType, class... Ts>
+CNumerics* AddCompressibilityCorrection(bool use_comp, Ts... args) {
+  if (use_comp) {
+    return new detail::CCompressibilityCorrection<BaseType>(args...);
+  } else {
+    return new BaseType(args...);
+  }
+}
+
+template <class BaseType, class... Ts>
+CNumerics* SAFactoryImpl(bool use_ft2, Ts... args) {
+  if (use_ft2) {
+    return AddCompressibilityCorrection<typename BaseType::template type<ft2::Nonzero>>(args...);
+  } else {
+    return AddCompressibilityCorrection<typename BaseType::template type<ft2::Zero>>(args...);
+  }
+}
+
+}  // namespace detail
 
 /*!
- * \class CSourcePieceWise_TurbSA_Neg
- * \brief Class for integrating the source terms of the negative Spalart-Allmaras model.
+ * \brief Creates an SA source based on the version and modifications/correction in the config.
  */
 template <class FlowIndices>
-using CSourcePieceWise_TurbSA_Neg = CSourceBase_TurbSA<FlowIndices, detail::Omega::Bsl, detail::ft2::Zero,
-                                                       detail::ModVort::Neg, detail::r::Bsl, detail::SourceTerms::Neg>;
+CNumerics* SAFactory(unsigned short nDim, const CConfig* config) {
+  const auto options = config->GetSAParsedOptions();
+
+  switch (options.version) {
+    case SA_OPTIONS::NONE:
+      return detail::SAFactoryImpl<detail::BaselineSA<FlowIndices>>(options.ft2, options.comp, nDim, config);
+    case SA_OPTIONS::NEG:
+      return detail::SAFactoryImpl<detail::NegativeSA<FlowIndices>>(options.ft2, options.comp, nDim, config);
+    case SA_OPTIONS::EDW:
+      return detail::SAFactoryImpl<detail::EdwardsSA<FlowIndices>>(options.ft2, options.comp, nDim, config);
+    default:
+      break;
+  }
+  return nullptr;
+}
 
 /*!
  * \class CSourcePieceWise_TurbSST
diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp
index d0a73ea43fa..e0a89838b21 100644
--- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp
+++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp
@@ -99,7 +99,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
     prandtlTurb(config.GetPrandtl_Turb()),
     cp(gamma * gasConst / (gamma - 1)),
     correct(iMesh == MESH_0),
-    useSA_QCR(config.GetQCR()),
+    useSA_QCR(config.GetSAParsedOptions().qcr2000),
     wallFun(config.GetWall_Functions()),
     uq(config.GetSSTParsedOptions().uq),
     uq_permute(config.GetUQ_Permute()),
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 68efa356d5e..bf589df95a7 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -1482,7 +1482,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Fluid_Interface(CGeometry* geometry,
 
               /*--- Turbulent kinetic energy ---*/
 
-              if (config->GetKind_Turb_Model() == TURB_MODEL::SST) 
+              if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
                 visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0),
                                                     solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0));
 
@@ -2459,7 +2459,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   const su2double Prandtl_Lam = config->GetPrandtl_Lam();
   const bool energy = config->GetEnergy_Equation();
-  const bool QCR = config->GetQCR();
+  const bool QCR = config->GetSAParsedOptions().qcr2000;
   const bool axisymmetric = config->GetAxisymmetric();
   const bool roughwall = (config->GetnRoughWall() > 0);
   const bool nemo = config->GetNEMOProblem();
@@ -2585,7 +2585,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /*--- Compute non-dimensional velocity and y+ ---*/
 
       FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
-      
+
       if (!wallfunctions && (MGLevel == MESH_0 || geometry->nodes->GetDomain(iPoint))) {
         // for CMultiGridGeometry, the normal neighbor of halo nodes in not set
         iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
@@ -2749,7 +2749,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   AllBoundViscCoeff.CEff = AllBoundViscCoeff.CL / (AllBoundViscCoeff.CD + EPS);
   AllBoundViscCoeff.CMerit = AllBoundViscCoeff.CT / (AllBoundViscCoeff.CQ + EPS);
-  AllBound_MaxHF_Visc = pow(AllBound_MaxHF_Visc, 1.0 / MaxNorm);
 
 #ifdef HAVE_MPI
 
@@ -2784,7 +2783,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     AllBoundViscCoeff.CMerit = AllBoundViscCoeff.CT / (AllBoundViscCoeff.CQ + EPS);
 
     AllBound_HF_Visc = Allreduce(AllBound_HF_Visc);
-    AllBound_MaxHF_Visc = pow(Allreduce(pow(AllBound_MaxHF_Visc, MaxNorm)), 1.0 / MaxNorm);
+    AllBound_MaxHF_Visc = Allreduce(AllBound_MaxHF_Visc);
   }
 
   /*--- Add the forces on the surfaces using all the nodes ---*/
@@ -2826,6 +2825,13 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
 #endif
 
+  /*--- Complete the calculation of maximum heat flux. ---*/
+
+  for (auto& hf : Surface_MaxHF_Visc) {
+    hf = pow(hf, 1.0 / MaxNorm);
+  }
+  AllBound_MaxHF_Visc = pow(AllBound_MaxHF_Visc, 1.0 / MaxNorm);
+
   /*--- Update the total coefficients (note that all the nodes have the same value)---*/
 
   TotalCoeff.CD += AllBoundViscCoeff.CD;
@@ -2863,7 +2869,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     SurfaceCoeff.CMz[iMarker_Monitoring] += SurfaceViscCoeff.CMz[iMarker_Monitoring];
   }
 
-
   Buffet_Monitoring(geometry, config);
 
 }
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index f945e9df183..76dcff59cb9 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -1225,21 +1225,20 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset,
   const int conv_bound_term = CONV_BOUND_TERM + offset;
   const int visc_bound_term = VISC_BOUND_TERM + offset;
 
-  bool spalart_allmaras, neg_spalart_allmaras, e_spalart_allmaras, comp_spalart_allmaras, e_comp_spalart_allmaras, menter_sst;
-  spalart_allmaras = neg_spalart_allmaras = e_spalart_allmaras = comp_spalart_allmaras = e_comp_spalart_allmaras = menter_sst = false;
-
   /*--- Assign turbulence model booleans ---*/
 
+  bool spalart_allmaras = false, menter_sst = false;
+
   switch (config->GetKind_Turb_Model()) {
     case TURB_MODEL::NONE:
       SU2_MPI::Error("No turbulence model selected.", CURRENT_FUNCTION);
       break;
-    case TURB_MODEL::SA:        spalart_allmaras = true;        break;
-    case TURB_MODEL::SA_NEG:    neg_spalart_allmaras = true;    break;
-    case TURB_MODEL::SA_E:      e_spalart_allmaras = true;      break;
-    case TURB_MODEL::SA_COMP:   comp_spalart_allmaras = true;   break;
-    case TURB_MODEL::SA_E_COMP: e_comp_spalart_allmaras = true; break;
-    case TURB_MODEL::SST:       menter_sst = true;              break;
+    case TURB_MODEL::SA:
+      spalart_allmaras = true;
+      break;
+    case TURB_MODEL::SST:
+      menter_sst = true;
+      break;
   }
 
   /*--- If the Menter SST model is used, store the constants of the model and determine the
@@ -1262,10 +1261,11 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset,
       break;
     case SPACE_UPWIND :
       for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) {
-        if (spalart_allmaras || neg_spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) {
+        if (spalart_allmaras) {
           numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSA<Indices>(nDim, nVar_Turb, config);
         }
-        else if (menter_sst) numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSST<Indices>(nDim, nVar_Turb, config);
+        else if (menter_sst)
+          numerics[iMGlevel][TURB_SOL][conv_term] = new CUpwSca_TurbSST<Indices>(nDim, nVar_Turb, config);
       }
       break;
     default:
@@ -1276,11 +1276,13 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset,
   /*--- Definition of the viscous scheme for each equation and mesh level ---*/
 
   for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) {
-    if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) {
-      numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA<Indices>(nDim, nVar_Turb, true, config);
+    if (spalart_allmaras) {
+      if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) {
+        numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg<Indices>(nDim, nVar_Turb, true, config);
+      } else {
+        numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA<Indices>(nDim, nVar_Turb, true, config);
+      }
     }
-    else if (neg_spalart_allmaras)
-      numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSA_Neg<Indices>(nDim, nVar_Turb, true, config);
     else if (menter_sst)
       numerics[iMGlevel][TURB_SOL][visc_term] = new CAvgGrad_TurbSST<Indices>(nDim, nVar_Turb, constants, true, config);
   }
@@ -1291,15 +1293,7 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset,
     auto& turb_source_first_term = numerics[iMGlevel][TURB_SOL][source_first_term];
 
     if (spalart_allmaras)
-      turb_source_first_term = new CSourcePieceWise_TurbSA<Indices>(nDim, nVar_Turb, config);
-    else if (e_spalart_allmaras)
-      turb_source_first_term = new CSourcePieceWise_TurbSA_E<Indices>(nDim, nVar_Turb, config);
-    else if (comp_spalart_allmaras)
-      turb_source_first_term = new CSourcePieceWise_TurbSA_COMP<Indices>(nDim, nVar_Turb, config);
-    else if (e_comp_spalart_allmaras)
-      turb_source_first_term = new CSourcePieceWise_TurbSA_E_COMP<Indices>(nDim, nVar_Turb, config);
-    else if (neg_spalart_allmaras)
-      turb_source_first_term = new CSourcePieceWise_TurbSA_Neg<Indices>(nDim, nVar_Turb, config);
+      turb_source_first_term = SAFactory<Indices>(nDim, config);
     else if (menter_sst)
       turb_source_first_term = new CSourcePieceWise_TurbSST<Indices>(nDim, nVar_Turb, constants, kine_Inf, omega_Inf,
                                                                      config);
@@ -1310,13 +1304,14 @@ void CDriver::InstantiateTurbulentNumerics(unsigned short nVar_Turb, int offset,
   /*--- Definition of the boundary condition method ---*/
 
   for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) {
-    if (spalart_allmaras || e_spalart_allmaras || comp_spalart_allmaras || e_comp_spalart_allmaras) {
+    if (spalart_allmaras) {
       numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA<Indices>(nDim, nVar_Turb, config);
-      numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA<Indices>(nDim, nVar_Turb, false, config);
-    }
-    else if (neg_spalart_allmaras) {
-      numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSA<Indices>(nDim, nVar_Turb, config);
-      numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg<Indices>(nDim, nVar_Turb, false, config);
+
+      if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) {
+        numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA_Neg<Indices>(nDim, nVar_Turb, false, config);
+      } else {
+        numerics[iMGlevel][TURB_SOL][visc_bound_term] = new CAvgGrad_TurbSA<Indices>(nDim, nVar_Turb, false, config);
+      }
     }
     else if (menter_sst) {
       numerics[iMGlevel][TURB_SOL][conv_bound_term] = new CUpwSca_TurbSST<Indices>(nDim, nVar_Turb, config);
diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
index 789e2423118..d95583a52ad 100644
--- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
@@ -436,7 +436,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config)
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
                   Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
-  if (config->GetQCR()) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
+  if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
   SetHeatFluxVector(Mean_GradPrimVar, Mean_Laminar_Viscosity,
@@ -620,7 +620,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi
   /*--- Get projected flux tensor (viscous residual) ---*/
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
                   Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
-  if (config->GetQCR()) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
+  if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
   GetViscousIncProjFlux(Mean_GradPrimVar, Normal, Mean_Thermal_Conductivity);
@@ -951,7 +951,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
                   Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
-  if (config->GetQCR()) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
+  if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
   SetHeatFluxVector(Mean_GradPrimVar, Mean_Laminar_Viscosity,
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
index 01f5b1ddc46..3ba641f2f88 100644
--- a/SU2_CFD/src/output/CFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -327,7 +327,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
 
   }
 
-  if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC){
+  if (config->GetSAParsedOptions().bc) {
     AddVolumeOutput("INTERMITTENCY", "gamma_BC", "INTERMITTENCY", "Intermittency");
   }
 
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 99bca4bd7a6..027c9a20699 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -977,7 +977,7 @@ void CFlowOutput::SetVolumeOutputFields_ScalarLimiter(const CConfig* config) {
     AddVolumeOutput("EDDY_VISCOSITY", "Eddy_Viscosity", "PRIMITIVE", "Turbulent eddy viscosity");
   }
 
-  if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+  if (config->GetSAParsedOptions().bc) {
     AddVolumeOutput("INTERMITTENCY", "gamma_BC", "INTERMITTENCY", "Intermittency");
   }
 
@@ -1047,7 +1047,7 @@ void CFlowOutput::LoadVolumeData_Scalar(const CConfig* config, const CSolver* co
     SetVolumeOutputValue("EDDY_VISCOSITY", iPoint, Node_Flow->GetEddyViscosity(iPoint));
   }
 
-  if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+  if (config->GetSAParsedOptions().bc) {
     SetVolumeOutputValue("INTERMITTENCY", iPoint, Node_Turb->GetGammaBC(iPoint));
   }
 
@@ -2082,26 +2082,15 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo
       switch (Kind_Turb_Model) {
         case TURB_MODEL::NONE: break;
         case TURB_MODEL::SA:
+          /// TODO: add the submodels here
           file << "Spalart Allmaras\n";
           break;
-        case TURB_MODEL::SA_NEG:
-          file << "Negative Spalart Allmaras\n";
-          break;
-        case TURB_MODEL::SA_E:
-          file << "Edwards Spalart Allmaras\n";
-          break;
-        case TURB_MODEL::SA_COMP:
-          file << "Compressibility Correction Spalart Allmaras\n";
-          break;
-        case TURB_MODEL::SA_E_COMP:
-          file << "Compressibility Correction Edwards Spalart Allmaras\n";
-          break;
         case TURB_MODEL::SST:
-          file << "Menter's SST\n";
-          // add the submodels here
-          //file << "Menter's SST with sustaining terms\n";
+          /// TODO: add the submodels here
           if (config->GetSSTParsedOptions().sust)
-            cout << "Menter's SST with sustaining terms" << endl; 
+            file << "Menter's SST with sustaining terms\n";
+          else
+            file << "Menter's SST\n";
          break;
       }
       break;
@@ -3279,7 +3268,7 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo
 }
 
 bool CFlowOutput::WriteVolume_Output(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile){
-  
+
   bool writeRestart = false;
   auto FileFormat = config->GetVolumeOutputFiles();
 
@@ -3289,9 +3278,9 @@ bool CFlowOutput::WriteVolume_Output(CConfig *config, unsigned long Iter, bool f
       return true;
     }
 
-    /* check if we want to write a restart file*/ 
+    /* check if we want to write a restart file*/
     if (FileFormat[iFile] == OUTPUT_TYPE::RESTART_ASCII || FileFormat[iFile] == OUTPUT_TYPE::RESTART_BINARY || FileFormat[iFile] == OUTPUT_TYPE::CSV) {
-      writeRestart = true;  
+      writeRestart = true;
     }
 
     /* only write 'double' files for the restart files */
diff --git a/SU2_CFD/src/output/output_physics.cpp b/SU2_CFD/src/output/output_physics.cpp
index 6bad7f8bb74..60909d216ba 100644
--- a/SU2_CFD/src/output/output_physics.cpp
+++ b/SU2_CFD/src/output/output_physics.cpp
@@ -142,7 +142,7 @@ void COutputLegacy::ComputeTurboPerformance(CSolver *solver_container, CGeometry
         else{
           nu = solver_container->GetNuIn(iMarkerTP, iSpan);
           NuFactorIn[iMarkerTP][iSpan]          = nu*DensityIn[iMarkerTP][iSpan]/muLam;
-          if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+          if (config->GetSAParsedOptions().bc) {
             NuFactorIn[iMarkerTP][iSpan]        = nu*DensityIn[iMarkerTP][iSpan]/muLam/0.005;
           }
         }
@@ -229,7 +229,7 @@ void COutputLegacy::ComputeTurboPerformance(CSolver *solver_container, CGeometry
         else{
           nu = solver_container->GetNuOut(iMarkerTP, iSpan);
           NuFactorOut[iMarkerTP][iSpan]          = nu*DensityOut[iMarkerTP][iSpan]/muLam;
-          if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+          if (config->GetSAParsedOptions().bc) {
             NuFactorOut[iMarkerTP][iSpan]        = nu*DensityOut[iMarkerTP][iSpan]/muLam/0.005;
           }
         }
diff --git a/SU2_CFD/src/output/output_structure_legacy.cpp b/SU2_CFD/src/output/output_structure_legacy.cpp
index e11a09a861a..560b5c8483a 100644
--- a/SU2_CFD/src/output/output_structure_legacy.cpp
+++ b/SU2_CFD/src/output/output_structure_legacy.cpp
@@ -476,7 +476,7 @@ void COutputLegacy::SetConvHistory_Header(ofstream *ConvHist_file, CConfig *conf
   char flow_resid[]= ",\"Res_Flow[0]\",\"Res_Flow[1]\",\"Res_Flow[2]\",\"Res_Flow[3]\",\"Res_Flow[4]\"";
   char adj_flow_resid[]= ",\"Res_AdjFlow[0]\",\"Res_AdjFlow[1]\",\"Res_AdjFlow[2]\",\"Res_AdjFlow[3]\",\"Res_AdjFlow[4]\"";
   switch (config->GetKind_Turb_Model()) {
-    case TURB_MODEL::SA:case TURB_MODEL::SA_NEG:case TURB_MODEL::SA_E: case TURB_MODEL::SA_COMP: case TURB_MODEL::SA_E_COMP:
+    case TURB_MODEL::SA:
       SPRINTF (turb_resid, ",\"Res_Turb[0]\"");
       break;
     case TURB_MODEL::SST:
@@ -485,7 +485,7 @@ void COutputLegacy::SetConvHistory_Header(ofstream *ConvHist_file, CConfig *conf
     default: break;
   }
   switch (config->GetKind_Turb_Model()) {
-    case TURB_MODEL::SA:case TURB_MODEL::SA_NEG:case TURB_MODEL::SA_E: case TURB_MODEL::SA_COMP: case TURB_MODEL::SA_E_COMP:
+    case TURB_MODEL::SA:
       SPRINTF (adj_turb_resid, ",\"Res_AdjTurb[0]\"");
       break;
     case TURB_MODEL::SST:
@@ -862,7 +862,7 @@ void COutputLegacy::SetConvHistory_Body(ofstream *ConvHist_file,
     if (compressible) nVar_Flow = nDim+2; else nVar_Flow = nDim+2;
     if (turbulent) {
       switch (config[val_iZone]->GetKind_Turb_Model()) {
-        case TURB_MODEL::SA: case TURB_MODEL::SA_NEG: case TURB_MODEL::SA_E: case TURB_MODEL::SA_E_COMP: case TURB_MODEL::SA_COMP: nVar_Turb = 1; break;
+        case TURB_MODEL::SA: nVar_Turb = 1; break;
         case TURB_MODEL::SST: nVar_Turb = 2; break;
         default: break;
       }
@@ -884,7 +884,7 @@ void COutputLegacy::SetConvHistory_Body(ofstream *ConvHist_file,
     if (compressible) nVar_AdjFlow = nDim+2; else nVar_AdjFlow = nDim+2;
     if (turbulent) {
       switch (config[val_iZone]->GetKind_Turb_Model()) {
-        case TURB_MODEL::SA: case TURB_MODEL::SA_NEG: case TURB_MODEL::SA_E: case TURB_MODEL::SA_E_COMP: case TURB_MODEL::SA_COMP: nVar_AdjTurb = 1; break;
+        case TURB_MODEL::SA: nVar_AdjTurb = 1; break;
         case TURB_MODEL::SST: nVar_AdjTurb = 2; break;
         default: break;
       }
@@ -1782,7 +1782,7 @@ void COutputLegacy::SetConvHistory_Body(ofstream *ConvHist_file,
             else cout << "      Res[Rho]";//, cout << "     Res[RhoE]";
 
             switch (config[val_iZone]->GetKind_Turb_Model()) {
-              case TURB_MODEL::SA: case TURB_MODEL::SA_NEG: case TURB_MODEL::SA_E: case TURB_MODEL::SA_E_COMP: case TURB_MODEL::SA_COMP:        cout << "       Res[nu]"; break;
+              case TURB_MODEL::SA: cout << "       Res[nu]"; break;
               case TURB_MODEL::SST: cout << "     Res[kine]" << "    Res[omega]"; break;
               default: break;
             }
@@ -2820,11 +2820,7 @@ void COutputLegacy::SpecialOutput_ForcesBreakdown(CSolver *****solver, CGeometry
         Breakdown_file << "Turbulence model: ";
         switch (Kind_Turb_Model) {
           case TURB_MODEL::SA:        Breakdown_file << "Spalart Allmaras" << "\n"; break;
-          case TURB_MODEL::SA_NEG:    Breakdown_file << "Negative Spalart Allmaras" << "\n"; break;
-          case TURB_MODEL::SA_E:      Breakdown_file << "Edwards Spalart Allmaras" << "\n"; break;
-          case TURB_MODEL::SA_COMP:   Breakdown_file << "Compressibility Correction Spalart Allmaras" << "\n"; break;
-          case TURB_MODEL::SA_E_COMP: Breakdown_file << "Compressibility Correction Edwards Spalart Allmaras" << "\n"; break;
-          case TURB_MODEL::SST:       Breakdown_file << "Menter's TURB_MODEL::SST"     << "\n"; break;
+          case TURB_MODEL::SST:       Breakdown_file << "Menter's SST"     << "\n"; break;
           default: break;
         }
         break;
diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp
index be10e5cb275..a2f57f7cc64 100644
--- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp
@@ -774,7 +774,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
 
             /*--- Turbulent kinetic energy ---*/
             /// TODO: This does not seem to be consistent with the primal treatment.
-            if (config->GetKind_Turb_Model() == TURB_MODEL::SST) 
+            if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
               val_turb_ke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0);
             else
               val_turb_ke = 0.0;
diff --git a/SU2_CFD/src/solvers/CBaselineSolver.cpp b/SU2_CFD/src/solvers/CBaselineSolver.cpp
index ea78922cf71..6dc757a45fc 100644
--- a/SU2_CFD/src/solvers/CBaselineSolver.cpp
+++ b/SU2_CFD/src/solvers/CBaselineSolver.cpp
@@ -421,7 +421,7 @@ void CBaselineSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf
         /*--- First, remove any variables for the turbulence model that
          appear in the restart file before the grid velocities. ---*/
 
-        if (turb_model == TURB_MODEL::SA || turb_model == TURB_MODEL::SA_NEG) {
+        if (turb_model == TURB_MODEL::SA) {
           index++;
         } else if (turb_model == TURB_MODEL::SST) {
           index+=2;
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index 05ec686661e..b0d76b496c8 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -4499,7 +4499,7 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container,
 
         /*--- Turbulent kinetic energy ---*/
 
-        if (config->GetKind_Turb_Model() == TURB_MODEL::SST) 
+        if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0),
                                               solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0));
 
@@ -6455,19 +6455,17 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
                             CConfig *config, unsigned short val_marker) {
   unsigned short iDim;
   unsigned long iVertex, iPoint;
-  su2double P_Total, T_Total, Velocity[3], Velocity2, H_Total, Temperature, Riemann,
+  su2double P_Total, T_Total, Velocity[MAXNDIM], Velocity2, H_Total, Temperature, Riemann,
   Pressure, Density, Energy, *Flow_Dir, Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag,
-  alpha, aa, bb, cc, dd, Area, UnitNormal[3];
+  alpha, aa, bb, cc, dd, Area, UnitNormal[MAXNDIM], Normal[MAXNDIM] = {0.0};
   su2double *V_inlet, *V_domain;
 
-  bool implicit             = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-  su2double Two_Gamma_M1       = 2.0/Gamma_Minus_One;
-  su2double Gas_Constant       = config->GetGas_ConstantND();
-  INLET_TYPE Kind_Inlet= config->GetKind_Inlet();
-  string Marker_Tag         = config->GetMarker_All_TagBound(val_marker);
-  bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
-
-  su2double *Normal = new su2double[nDim];
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const su2double Two_Gamma_M1 = 2.0 / Gamma_Minus_One;
+  const su2double Gas_Constant = config->GetGas_ConstantND();
+  const auto Kind_Inlet = config->GetKind_Inlet();
+  const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+  const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
 
   /*--- Loop over all the vertices on this boundary marker ---*/
 
@@ -6744,10 +6742,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
   }
   END_SU2_OMP_FOR
 
-  /*--- Free locally allocated memory ---*/
-
-  delete [] Normal;
-
 }
 
 void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
@@ -6927,147 +6921,89 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
 void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_container,
                                        CNumerics *conv_numerics, CNumerics *visc_numerics,
                                        CConfig *config, unsigned short val_marker) {
-  unsigned short iDim;
-  unsigned long iVertex, iPoint;
-  su2double *V_inlet, *V_domain;
-
-  su2double Density, Energy, Velocity2;
-  su2double Gas_Constant = config->GetGas_ConstantND();
-
-  bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-  string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
-  bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
-  su2double *Normal = new su2double[nDim];
-  su2double *Velocity = new su2double[nDim];
+  const su2double Gas_Constant = config->GetGas_ConstantND();
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+  const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
 
   /*--- Supersonic inlet flow: there are no outgoing characteristics,
    so all flow variables can be imposed at the inlet.
    First, retrieve the specified values for the primitive variables. ---*/
 
-  auto Temperature = config->GetInlet_Temperature(Marker_Tag);
-  auto Pressure    = config->GetInlet_Pressure(Marker_Tag);
-  auto Vel         = config->GetInlet_Velocity(Marker_Tag);
-
-  /*--- Non-dim. the inputs if necessary. ---*/
+  const su2double Temperature = config->GetInlet_Temperature(Marker_Tag) / config->GetTemperature_Ref();
+  const su2double Pressure = config->GetInlet_Pressure(Marker_Tag) / config->GetPressure_Ref();
+  const auto* Vel = config->GetInlet_Velocity(Marker_Tag);
 
-  Temperature /= config->GetTemperature_Ref();
-  Pressure    /= config->GetPressure_Ref();
-  for (iDim = 0; iDim < nDim; iDim++)
+  su2double Velocity[MAXNDIM] = {0.0};
+  for (unsigned short iDim = 0; iDim < nDim; iDim++)
     Velocity[iDim] = Vel[iDim] / config->GetVelocity_Ref();
 
   /*--- Density at the inlet from the gas law ---*/
 
-  Density = Pressure/(Gas_Constant*Temperature);
+  const su2double Density = Pressure / (Gas_Constant * Temperature);
 
   /*--- Compute the energy from the specified state ---*/
 
-  Velocity2 = 0.0;
-  for (iDim = 0; iDim < nDim; iDim++)
-    Velocity2 += Velocity[iDim]*Velocity[iDim];
-  Energy = Pressure/(Density*Gamma_Minus_One)+0.5*Velocity2;
+  const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity);
+  su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2;
   if (tkeNeeded) Energy += GetTke_Inf();
 
   /*--- Loop over all the vertices on this boundary marker ---*/
 
   SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
-  for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+  for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+    const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
-    /*--- Allocate the value at the outlet ---*/
+    if (!geometry->nodes->GetDomain(iPoint)) continue;
 
-    V_inlet = GetCharacPrimVar(val_marker, iVertex);
+    /*--- Allocate the value at the inlet ---*/
+
+    auto* V_inlet = GetCharacPrimVar(val_marker, iVertex);
 
     /*--- Primitive variables, using the derived quantities ---*/
 
-    V_inlet[0] = Temperature;
-    for (iDim = 0; iDim < nDim; iDim++)
-      V_inlet[iDim+1] = Velocity[iDim];
-    V_inlet[nDim+1] = Pressure;
-    V_inlet[nDim+2] = Density;
-    V_inlet[nDim+3] = Energy + Pressure/Density;
+    V_inlet[prim_idx.Temperature()] = Temperature;
+    V_inlet[prim_idx.Pressure()] = Pressure;
+    V_inlet[prim_idx.Density()] = Density;
+    V_inlet[prim_idx.Enthalpy()] = Energy + Pressure / Density;
+    for (unsigned short iDim = 0; iDim < nDim; iDim++)
+      V_inlet[iDim+prim_idx.Velocity()] = Velocity[iDim];
 
-    iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
+    /*--- Current solution at this boundary node ---*/
 
-    /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
+    auto* V_domain = nodes->GetPrimitive(iPoint);
 
-    if (geometry->nodes->GetDomain(iPoint)) {
+    /*--- Normal vector for this vertex (negate for outward convention) ---*/
 
-      /*--- Current solution at this boundary node ---*/
+    su2double Normal[MAXNDIM] = {0.0};
+    geometry->vertex[val_marker][iVertex]->GetNormal(Normal);
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
 
-      V_domain = nodes->GetPrimitive(iPoint);
+    /*--- Set various quantities in the solver class ---*/
 
-      /*--- Normal vector for this vertex (negate for outward convention) ---*/
+    conv_numerics->SetNormal(Normal);
+    conv_numerics->SetPrimitive(V_domain, V_inlet);
 
-      geometry->vertex[val_marker][iVertex]->GetNormal(Normal);
-      for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
-
-      /*--- Set various quantities in the solver class ---*/
+    if (dynamic_grid)
+      conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
+                                geometry->nodes->GetGridVel(iPoint));
 
-      conv_numerics->SetNormal(Normal);
-      conv_numerics->SetPrimitive(V_domain, V_inlet);
+    /*--- Compute the residual using an upwind scheme ---*/
 
-      if (dynamic_grid)
-        conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
-                                  geometry->nodes->GetGridVel(iPoint));
+    auto residual = conv_numerics->ComputeResidual(config);
 
-      /*--- Compute the residual using an upwind scheme ---*/
+    LinSysRes.AddBlock(iPoint, residual);
 
-      auto residual = conv_numerics->ComputeResidual(config);
+    /*--- Jacobian contribution for implicit integration ---*/
 
-      LinSysRes.AddBlock(iPoint, residual);
-
-      /*--- Jacobian contribution for implicit integration ---*/
-
-      if (implicit)
-        Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+    if (implicit)
+      Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
 
-//      /*--- Viscous contribution, commented out because serious convergence problems ---*/
-//
-//      if (viscous) {
-//
-//        /*--- Set laminar and eddy viscosity at the infinity ---*/
-//
-//        V_inlet[nDim+5] = nodes->GetLaminarViscosity(iPoint);
-//        V_inlet[nDim+6] = nodes->GetEddyViscosity(iPoint);
-//
-//        /*--- Set the normal vector and the coordinates ---*/
-//
-//        visc_numerics->SetNormal(Normal);
-//        su2double Coord_Reflected[MAXNDIM];
-//        GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal),
-//                                                 geometry->nodes->GetCoord(iPoint), Coord_Reflected);
-//        visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected);
-//
-//        /*--- Primitive variables, and gradient ---*/
-//
-//        visc_numerics->SetPrimitive(V_domain, V_inlet);
-//        visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
-//
-//        /*--- Turbulent kinetic energy ---*/
-//
-//        if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
-//          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0),
-//                                              solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0));
-//
-//        /*--- Compute and update residual ---*/
-//
-//        auto residual = visc_numerics->ComputeResidual(config);
-//        LinSysRes.SubtractBlock(iPoint, residual);
-//
-//        /*--- Jacobian contribution for implicit integration ---*/
-//
-//        if (implicit)
-//          Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i);
-//      }
+    /*--- Viscous contribution, omited to improve convergence. ---*/
 
-    }
   }
   END_SU2_OMP_FOR
 
-  /*--- Free locally allocated memory ---*/
-
-  delete [] Normal;
-  delete [] Velocity;
-
 }
 
 void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_container,
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index d6fc198097a..42c8f32bed4 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -110,7 +110,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
 
   Factor_nu_Inf = config->GetNuFactor_FreeStream();
   su2double nu_tilde_Inf  = Factor_nu_Inf*Viscosity_Inf/Density_Inf;
-  if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+  if (config->GetSAParsedOptions().bc) {
     nu_tilde_Inf  = 0.005*Factor_nu_Inf*Viscosity_Inf/Density_Inf;
   }
 
@@ -119,7 +119,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
   /*--- Factor_nu_Engine ---*/
   Factor_nu_Engine = config->GetNuFactor_Engine();
   nu_tilde_Engine  = Factor_nu_Engine*Viscosity_Inf/Density_Inf;
-  if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC) {
+  if (config->GetSAParsedOptions().bc) {
     nu_tilde_Engine  = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf;
   }
 
@@ -217,7 +217,7 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain
 
   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() == TURB_MODEL::SA_NEG);
+  const bool neg_spalart_allmaras = config->GetSAParsedOptions().version == SA_OPTIONS::NEG;
 
   /*--- Compute eddy viscosity ---*/
 
@@ -226,22 +226,20 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) {
 
-    su2double rho = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-    su2double mu  = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
+    const su2double rho = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
+    const su2double mu = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
 
-    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;
+    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;
 
-    su2double Ji   = nu_hat/nu ;
+    su2double Ji = nu_hat/nu;
     if (roughness > 1.0e-10)
-      Ji+= cR1*roughness/(dist+EPS);
+      Ji += cR1*roughness/(dist+EPS);
 
-    su2double Ji_3 = Ji*Ji*Ji;
-    su2double fv1  = Ji_3/(Ji_3+cv1_3);
+    const su2double Ji_3 = Ji*Ji*Ji;
+    const su2double fv1  = Ji_3/(Ji_3+cv1_3);
 
     su2double muT = rho*fv1*nu_hat;
 
@@ -261,8 +259,7 @@ void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, C
   /*--- Define an object to set solver specific numerics contribution. ---*/
   auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {
     /*--- Roughness heights. ---*/
-    if (config->GetKind_Turb_Model() == TURB_MODEL::SA)
-      numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), geometry->nodes->GetRoughnessHeight(jPoint));
+    numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), geometry->nodes->GetRoughnessHeight(jPoint));
   };
 
   /*--- Now instantiate the generic implementation with the functor above. ---*/
@@ -275,7 +272,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);
-  const bool transition_BC = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::BC);
+  const bool transition_BC = config->GetSAParsedOptions().bc;
 
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
@@ -1573,7 +1570,7 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) {
    SA_NEG model is more robust due to allowing for negative nu_tilde,
    so the under-relaxation is not applied to that variant. */
 
-  if (config->GetKind_Turb_Model() == TURB_MODEL::SA_NEG) return;
+  if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) return;
 
   /* Loop over the solution update given by relaxing the linear
    system for this nonlinear iteration. */
diff --git a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg
index bcefda6e8ed..e240c75f768 100644
--- a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg
+++ b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg
@@ -12,19 +12,17 @@
 
 % ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
 %
-% Physical governing equations (EULER, NAVIER_STOKES,
-%                               TNE2_EULER, TNE2_NAVIER_STOKES,
-%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
-%                               POISSON_EQUATION)                             
+% Physical governing equations                             
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_NEG
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= NEGATIVE, EXPERIMENTAL
 %
-% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
+% Mathematical problem
 MATH_PROBLEM= DISCRETE_ADJOINT
 %
-% Restart solution (NO, YES)
+% Restart solution
 RESTART_SOL= NO
 
 % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa.cfg
index d9119125c7e..cc30febbadb 100644
--- a/TestCases/rans/naca0012/turb_NACA0012_sa.cfg
+++ b/TestCases/rans/naca0012/turb_NACA0012_sa.cfg
@@ -17,8 +17,9 @@
 %                               POISSON_EQUATION)                             
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_NEG
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= NEGATIVE, EXPERIMENTAL
 %
 % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
 MATH_PROBLEM= DIRECT
diff --git a/TestCases/transition/E387_Airfoil/transitional_BC_model_ConfigFile.cfg b/TestCases/transition/E387_Airfoil/transitional_BC_model_ConfigFile.cfg
index efd875c9f15..b90512ff607 100644
--- a/TestCases/transition/E387_Airfoil/transitional_BC_model_ConfigFile.cfg
+++ b/TestCases/transition/E387_Airfoil/transitional_BC_model_ConfigFile.cfg
@@ -17,10 +17,9 @@
 %                               POISSON_EQUATION)
 SOLVER= INC_RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
+% Specify turbulence model
 KIND_TURB_MODEL= SA
-% Specify transition model (NONE, LM, BC)
-KIND_TRANS_MODEL= BC
+SA_OPTIONS= BCM
 FREESTREAM_TURBULENCEINTENSITY = 0.1
 %
 % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
diff --git a/TestCases/transition/Schubauer_Klebanoff/transitional_BC_model_ConfigFile.cfg b/TestCases/transition/Schubauer_Klebanoff/transitional_BC_model_ConfigFile.cfg
index bbe4f98be94..a5aae5dcff7 100644
--- a/TestCases/transition/Schubauer_Klebanoff/transitional_BC_model_ConfigFile.cfg
+++ b/TestCases/transition/Schubauer_Klebanoff/transitional_BC_model_ConfigFile.cfg
@@ -2,11 +2,11 @@
 %                                                                              %
 % SU2 configuration file                                                       %
 % Case description: Schubauer-Klebanoff Flat Plate Natural Transition          %
-% Author: Samet Cakmakcioglu	                                                 %
-% Institution: TOBB University of Economics and Technology             		     %	
-%			   TAI-TUSAS Turkish Aerospace Industries 						                   %
+% Author: Samet Cakmakcioglu                                                   %
+% Institution: TOBB University of Economics and Technology                     %	
+%              TAI-TUSAS Turkish Aerospace Industries                          %
 % Date: Oct 10th, 2016                                                         %
-% File Version 7.3.1 "Blackbird"                                                %
+% File Version 7.3.1 "Blackbird"                                               %
 %                                                                              %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -17,10 +17,9 @@
 %                               POISSON_EQUATION)
 SOLVER= INC_RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
+% Specify turbulence model
 KIND_TURB_MODEL= SA
-% Specify transition model (NONE, LM, BC)
-KIND_TRANS_MODEL= BC
+SA_OPTIONS= BCM
 FREESTREAM_TURBULENCEINTENSITY = 0.18
 %
 % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
diff --git a/TestCases/transition/T3A_FlatPlate/transitional_BC_model_ConfigFile.cfg b/TestCases/transition/T3A_FlatPlate/transitional_BC_model_ConfigFile.cfg
index 87cfc745ab8..be20294ead5 100644
--- a/TestCases/transition/T3A_FlatPlate/transitional_BC_model_ConfigFile.cfg
+++ b/TestCases/transition/T3A_FlatPlate/transitional_BC_model_ConfigFile.cfg
@@ -17,10 +17,9 @@
 %                               POISSON_EQUATION)
 SOLVER= INC_RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
+% Specify turbulence model
 KIND_TURB_MODEL= SA
-% Specify transition model (NONE, LM, BC)
-KIND_TRANS_MODEL= BC
+SA_OPTIONS= BCM
 FREESTREAM_TURBULENCEINTENSITY = 3.0
 %
 % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
diff --git a/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_EDW_RAE2822.cfg b/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_EDW_RAE2822.cfg
index e6ee9301713..fb7b111b3d3 100644
--- a/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_EDW_RAE2822.cfg
+++ b/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_EDW_RAE2822.cfg
@@ -3,8 +3,9 @@
 % Physical governing equations                    
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_E_COMP
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= EDWARDS, COMPRESSIBILITY, EXPERIMENTAL
 %
 % Restart solution (NO, YES)
 RESTART_SOL= NO
diff --git a/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_RAE2822.cfg b/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_RAE2822.cfg
index 31aeaf9c3fa..aeb8e3dbbf1 100644
--- a/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_RAE2822.cfg
+++ b/TestCases/turbulence_models/sa/rae2822/turb_SA_COMP_RAE2822.cfg
@@ -3,8 +3,9 @@
 % Physical governing equations                    
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_COMP
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= COMPRESSIBILITY
 %
 % Restart solution (NO, YES)
 RESTART_SOL= NO
diff --git a/TestCases/turbulence_models/sa/rae2822/turb_SA_EDW_RAE2822.cfg b/TestCases/turbulence_models/sa/rae2822/turb_SA_EDW_RAE2822.cfg
index a8108cacf97..70c20176b7f 100644
--- a/TestCases/turbulence_models/sa/rae2822/turb_SA_EDW_RAE2822.cfg
+++ b/TestCases/turbulence_models/sa/rae2822/turb_SA_EDW_RAE2822.cfg
@@ -3,8 +3,9 @@
 % Physical governing equations                    
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_E
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= EDWARDS
 %
 % Restart solution (NO, YES)
 RESTART_SOL= NO
diff --git a/TestCases/turbulence_models/sa/rae2822/turb_SA_NEG_RAE2822.cfg b/TestCases/turbulence_models/sa/rae2822/turb_SA_NEG_RAE2822.cfg
index 62ad6723e87..039a0daffb8 100644
--- a/TestCases/turbulence_models/sa/rae2822/turb_SA_NEG_RAE2822.cfg
+++ b/TestCases/turbulence_models/sa/rae2822/turb_SA_NEG_RAE2822.cfg
@@ -19,8 +19,9 @@
 % Physical governing equations                    
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
-KIND_TURB_MODEL= SA_NEG
+% Specify turbulence model
+KIND_TURB_MODEL= SA
+SA_OPTIONS= NEGATIVE, EXPERIMENTAL
 %
 % Restart solution (NO, YES)
 RESTART_SOL= NO
diff --git a/TestCases/turbulence_models/sa/rae2822/turb_SA_QCR_RAE2822.cfg b/TestCases/turbulence_models/sa/rae2822/turb_SA_QCR_RAE2822.cfg
index f6246965822..05b09d9a911 100644
--- a/TestCases/turbulence_models/sa/rae2822/turb_SA_QCR_RAE2822.cfg
+++ b/TestCases/turbulence_models/sa/rae2822/turb_SA_QCR_RAE2822.cfg
@@ -3,14 +3,13 @@
 % Physical governing equations                    
 SOLVER= RANS
 %
-% Specify turbulent model (NONE, SA, SA_NEG, SST)
+% Specify turbulence model
 KIND_TURB_MODEL= SA
-%
-% Activate SA Quadratic Constitutive Relation, 2000 version
-SA_QCR= YES
+SA_OPTIONS= QCR2000
 %
 % Restart solution (NO, YES)
 RESTART_SOL= NO
+
 % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
 %
 % Mach number (non-dimensional, based on the free-stream values)
diff --git a/config_template.cfg b/config_template.cfg
index 78e2778b6fc..70e3e03f52a 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -18,13 +18,16 @@
 %              HEAT_EQUATION_FVM, ELASTICITY)
 SOLVER= EULER
 %
-% Specify turbulence model (NONE, SA, SA_NEG, SST, SA_E, SA_COMP, SA_E_COMP)
+% Specify turbulence model (NONE, SA, SST)
 KIND_TURB_MODEL= NONE
 %
-% Specify versions/corrections to SST model (V2003m, V1994m, VORTICITY, KATO_LAUNDER, UQ, SUSTAINING)
+% Specify versions/corrections of the SST model (V2003m, V1994m, VORTICITY, KATO_LAUNDER, UQ, SUSTAINING)
 SST_OPTIONS= NONE
 %
-% Transition model (NONE, BC)
+% Specify versions/corrections of the SA model (NEGATIVE, EDWARDS, WITHFT2, QCR2000, COMPRESSIBILITY, ROTATION, BCM, EXPERIMENTAL)
+SA_OPTIONS= NONE
+%
+% Transition model (NONE, LM)
 KIND_TRANS_MODEL= NONE
 %
 % Specify subgrid scale model(NONE, IMPLICIT_LES, SMAGORINSKY, WALE, VREMAN)