diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index bc32b5407b1..94a2d390e81 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1247,6 +1247,23 @@ class CConfig {
   string* user_scalar_names;          /*!< \brief Names of the user defined (auxiliary) transported scalars .*/
   string* user_source_names;          /*!< \brief Names of the source terms for the user defined transported scalars. */
 
+  /*--- Vortex Generator model config options ---*/  // Added by Max
+  string vg_filename;
+  su2double vg_constant;
+  ENUM_VG_MODEL vg_bay;
+
+  unsigned short nVgs=0;
+  su2double ***Vg_Coordinates=nullptr;
+  su2double **Vg_SurfaceNormalDirection=nullptr,
+            **Vg_SurfaceTangentialDirection=nullptr,
+            **Vg_SurfaceCrossFlowDirection=nullptr,
+            **Vg_StreamwiseDirection=nullptr;
+  su2double *Vg_SurfaceArea=nullptr,
+            *Vg_Angle=nullptr;
+  unsigned short nPointsVG = 4;
+
+  // End by Max
+
   /*!
    * \brief Set the default values of config options not set in the config file using another config object.
    * \param config - Config object to use the default values from.
@@ -1379,6 +1396,11 @@ class CConfig {
                            unsigned short& nMarker_ActDiskBemInlet, unsigned short& nMarker_ActDiskBemOutlet, string*& Marker_ActDiskBemInlet, string*& Marker_ActDiskBemOutlet,
                            su2double**& ActDiskBem_X, su2double**& ActDiskBem_Y, su2double**& ActDiskBem_Z);
 
+  //added by max
+  // void addVgOption(const string& name, unsigned short& nVgs, su2double***& Vg_Coordinates,
+  //                su2double**& Vg_SurfaceNormalDirection, su2double**& Vg_SurfaceTangentialDirection, su2double**& Vg_SurfaceCrossFlowDirection,su2double*& Vg_SurfaceArea,ENUM_VG_MODEL& bayModel);
+  //emd added by max
+
   void addWallFunctionOption(const string &name,               unsigned short &list_size,
                              string* &string_field,            WALL_FUNCTIONS* &val_Kind_WF,
                              unsigned short** &val_IntInfo_WF, su2double** &val_DoubleInfo_WF);
@@ -9859,4 +9881,126 @@ class CConfig {
    */
   LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; }
 
+  //Added by Max
+  /*!
+   * \brief Check which vortex generator model is in use.
+   * \return Vortex Genrator Model.
+   */
+  ENUM_VG_MODEL GetVGModel(void) const {return vg_bay; }
+
+  /*!
+   * \brief Get VG Model calibration constant.
+   * \return Bay model calibration constant.
+   */
+  su2double GetVGConstant(void) const {return vg_constant;}
+
+  /*!
+   * \brief Set vortex generator height vector.
+   * \param[in] b - height vector.
+   */
+  void Set_bVG(su2double** b) { Vg_SurfaceNormalDirection = b; }
+
+  /*!
+   * \brief Ger vortex generator height vector.
+   * \param[in] iVg - Vortex generator index
+   * \return height vector.
+   */
+  su2double* Get_bVG(unsigned short iVg) const { return Vg_SurfaceNormalDirection[iVg]; }
+
+  /*!
+   * \brief Set vortex generator normal vector.
+   * \param[in] n - height vector.
+   */
+  void Set_nVG(su2double** n)  { Vg_SurfaceCrossFlowDirection = n; }
+
+  /*!
+   * \brief Get vortex generator normal vector.
+   * \param[in] iVg - Vortex generator index
+   * \return height vector.
+   */
+  su2double* Get_nVG(unsigned short iVg) const { return Vg_SurfaceCrossFlowDirection[iVg]; }
+
+  /*!
+   * \brief Set vortex generator tangential vector.
+   * \param[in] t - tangential vector.
+   */
+  void Set_tVG(su2double** t) { Vg_SurfaceTangentialDirection = t; }
+
+   /*!
+   * \brief Get vortex generator tangential vector.
+   * \param[in] iVg - Vortex generator index
+   * \return tangential vector.
+   */
+  su2double* Get_tVG(unsigned short iVg) const { return Vg_SurfaceTangentialDirection[iVg]; }
+
+  /*!
+   * \brief Get vortex generator coordinates.
+   * \param[in] iVg - Vortex generator index
+   * \param[in] vg_coord - vortex generator coordinates.
+   */
+  su2double** GetVGcoord(unsigned short iVg) const {return Vg_Coordinates[iVg];}
+
+  /*!
+   * \brief Set vortex generator coordinates.
+   * \param[in] vg_coord - vortex generator coordinates.
+   */
+  void SetVGCoord(su2double ***vg_coord){Vg_Coordinates=vg_coord;}
+  /*!
+   * \brief Get number of vortex generators for a specific problem.
+   * \return - Number of vortex genrators for a specific zone.
+   */
+  unsigned short Get_nVGs(void) const {return nVgs;}
+
+  /*!
+   * \brief Set number of vortex generators for a specific problem.
+   * \param[in] nVgs_zone - Number of vortex genrators for a specific zone.
+   */
+  void Set_nVGs(unsigned short nVgs_zone){nVgs=nVgs_zone;}
+
+  /*!
+   * \brief Get number of vortex generators for a specific problem.
+   * \return - Number of vortex genrators for a specific zone.
+   */
+  unsigned short Get_nPointsVg() const {return nPointsVG;};
+  /*!
+   * \brief Get vortex generators area.
+   * \param[in] iVg - Vortex generator index
+   * \return Area of the vortex generators plane.
+   */
+  su2double Get_Svg(unsigned short iVg) const {return Vg_SurfaceArea[iVg];}
+
+  /*!
+   * \brief Set vortex generators area.
+   * \param[in] Svg - Areas of the vortex generators planes.
+   */
+  void Set_Svg(su2double* Svg){Vg_SurfaceArea=Svg;}
+
+  /*!
+   * \brief Set vortex generators streawise direction.
+   * \param[in] Vg_StreamwiseDirection - Vortex generators streawise direction.
+   */
+  void Set_uhatVg(su2double** uhat){Vg_StreamwiseDirection=uhat;}
+
+  /*!
+   * \brief Get vortex generators streawise direction.
+   * \param[in] iVg - Vortex generator index
+   * \return Vortex generators streawise direction.
+   */
+  su2double* Get_uhatVg(unsigned short iVG) const {return Vg_StreamwiseDirection[iVG];}
+  /*!
+   * \brief Set vortex generators angles.
+   * \param[in] beta - Vortex generators angles.
+   */
+  void Set_betaVg(su2double *beta){Vg_Angle=beta;}
+
+  /*!
+   * \brief Get vortex generators angles.
+   * \param[in] iVg - Vortex generator index
+   * \return Vortex generator angle.
+   */
+  su2double Get_betaVg(unsigned short iVg) const {return Vg_Angle[iVg];}
+
+
+  const string& GetVGFileName() const {return vg_filename;}
+  //End added by max
 };
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 919aefa65b9..6d084bc24fc 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -2571,6 +2571,19 @@ static const MapType<std::string, ENUM_PROJECTION_FUNCTION> Projection_Function_
   MakePair("HEAVISIDE_UP", ENUM_PROJECTION_FUNCTION::HEAVISIDE_UP)
   MakePair("HEAVISIDE_DOWN", ENUM_PROJECTION_FUNCTION::HEAVISIDE_DOWN)
 };
+//Added by max
+
+enum class ENUM_VG_MODEL{
+  NONE,
+  BAY,
+  JBAY
+};
+static const MapType<std::string,ENUM_VG_MODEL> VgModel_Map={
+  MakePair("NONE",ENUM_VG_MODEL::NONE)
+  MakePair("BAY",ENUM_VG_MODEL::BAY)
+  MakePair("JBAY",ENUM_VG_MODEL::JBAY)
+};
+//End added by max
 
 /*!
  * \brief the different validation solution
diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp
index ce0baf55bf6..f375581ff20 100644
--- a/Common/include/toolboxes/geometry_toolbox.hpp
+++ b/Common/include/toolboxes/geometry_toolbox.hpp
@@ -216,5 +216,100 @@ inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector,
 
   for (Int iDim = 0; iDim < nDim; iDim++) proj[iDim] -= normalProj * vector[iDim];
 }
+// Added by max
+
+/*!
+ * \brief Check if the edge intersects the plane
+ * \param[in] nDim - Number of dimensions of the particular problem.
+ * \param[in] p0 - Point defining the plane.
+ * \param[in] n - Plane normal vector.
+ * \param[in] p1,p2 - Points defining the edge
+ * \return Returns <code>TRUE</code> if the edge intersect the plane, <code>FALSE</code> it doesn't.
+ */
+template <class T, class Int>
+inline bool IntersectEdge(Int nDim, const T* p0, const T* n, const T* p1, const T* p2) {
+  T d = -DotProduct(nDim, n, p0);
+  auto a = DotProduct(nDim, n, p1) + d;
+  auto b = DotProduct(nDim, n, p2) + d;
+  if (a * b <= 0 && (a != 0 || b != 0))
+    return true;
+  else
+    return false;
+}
+
+/*!
+ * \brief Check if a point is inside a polygon
+ * \param[in] nDim - Number of dimensions of the particular problem.
+ * \param[in] vector - Vector to normalize.
+ */
+template <class T, class Int>
+inline void NormalizeVector(Int nDim, T* vector) {
+  T norm = Norm(nDim, vector);
+  for (Int i = 0; i < nDim; i++) vector[i] /= norm;
+}
+
+/*!
+ * \brief Check if a point is inside a polygon
+ * \param[in] nDim - Number of dimensions of the particular problem.
+ * \param[in] pVert - Point defining the polygon.
+ * \param[in] p0 - Point to test
+ * \param[in] nVert - Number of points that form the polygon.
+ * \return Returns TRUE if the edge intersect the plane, FALSE it doesn't.
+ */
+template <class Mat, class T, class Int>
+inline bool PointInConvexPolygon(Int nDim, const Mat& pVert, const T* p0, int nVert) {
+  unsigned short i = 0, j = 1;
+  unsigned short idxPoint1 = 0, idxPoint2 = 1;
+  unsigned short nIntersections = 0;
+  unsigned short proj_index = 0;
+  T max_proj = 0;
+
+  /* Check which of the x,y,z planes is the most parallel to the polygon */
+
+  if (nDim == 3) {
+    T plane_norm[3];
+    TriangleNormal(pVert, plane_norm);
+    NormalizeVector(3, plane_norm);
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      if (abs(plane_norm[iDim]) > max_proj) {
+        proj_index = iDim;
+        max_proj = abs(plane_norm[iDim]);
+      }
+    }
+
+    switch (proj_index) {
+      case 0:
+        i = 1;
+        j = 2;
+        break;
+      case 1:
+        i = 2;
+        j = 0;
+        break;
+      default:
+        i = 0;
+        j = 1;
+        break;
+    }
+  }
+
+  /* Loop across the polygon edges and check if the ray intersect a vertex */
+
+  for (unsigned short iVert = 1; iVert < nVert + 1; iVert++) {
+    idxPoint2 = iVert % nVert;
+
+    if (((p0[j] < pVert[idxPoint1][j]) != (p0[j] < pVert[idxPoint2][j])) &&
+        (p0[i] < pVert[idxPoint1][i] + ((p0[j] - pVert[idxPoint1][j]) / (pVert[idxPoint2][j] - pVert[idxPoint1][j])) *
+                                           (pVert[idxPoint2][i] - pVert[idxPoint1][i]))) {
+      nIntersections++;
+    }
+    idxPoint1 = idxPoint2;
+  }
+
+  return nIntersections % 2 == 1;
+}
+
+// end added by max
 /// @}
 }  // namespace GeometryToolbox
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 63810ab3e03..bb597929bf1 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -570,6 +570,18 @@ void CConfig::addActDiskBemOption(const string& name,
   option_map.insert(pair<string, COptionBase *>(name, val));
 }
 
+// // added by max
+// void CConfig::addVgOption(const string& name, unsigned short& nVgs,
+//                           su2double***& Vg_Coordinates, su2double**& Vg_SurfaceNormalDirection,
+//                           su2double**& Vg_SurfaceTangentialDirection, su2double**& Vg_SurfaceCrossFlowDirection,su2double*& Vg_SurfaceArea,ENUM_VG_MODEL& bayModel) {
+//   assert(option_map.find(name) == option_map.end());
+//   all_options.insert(pair<string, bool>(name, true));
+//   COptionBase* val = new COptionVGmodel(name, nVgs, Vg_Coordinates, Vg_SurfaceNormalDirection,
+//                                         Vg_SurfaceTangentialDirection, Vg_SurfaceCrossFlowDirection,Vg_SurfaceArea,bayModel);
+//   option_map.insert(pair<string, COptionBase*>(name, val));
+// }
+// // end added by max
+
 void CConfig::addWallFunctionOption(const string &name, unsigned short &list_size, string* &string_field,
                                     WALL_FUNCTIONS* &val_Kind_WF, unsigned short** &val_IntInfo_WF,
                                     su2double** &val_DoubleInfo_WF) {
@@ -2961,6 +2973,20 @@ void CConfig::SetConfig_Options() {
   /*!\brief ROM_SAVE_FREQ \n DESCRIPTION: How often to save snapshots for unsteady problems.*/
   addUnsignedShortOption("ROM_SAVE_FREQ", rom_save_freq, 1);
 
+//Added by Max
+
+  /*!\brief VG_CONST \n DESCRIPTION: Calibration constant for VG model.*/
+  addEnumOption("VG_MODEL", vg_bay,VgModel_Map, ENUM_VG_MODEL::NONE);
+
+  /*!\brief VG_CONST \n DESCRIPTION: Calibration constant for VG model.*/
+  addDoubleOption("VG_CONST", vg_constant, 10);
+
+  /*!\brief VG_CONFIG \n DESCRIPTION: VG configuration file name*/
+  addStringOption("VG_CONFIG", vg_filename, " ");
+  // addVgOption("VG_CONFIG",nVgs,Vg_Coordinates,Vg_SurfaceNormalDirection,Vg_SurfaceTangentialDirection,Vg_SurfaceCrossFlowDirection,Vg_SurfaceArea,vg_bay);
+
+  //End Max
+
   /* END_CONFIG_OPTIONS */
 
 }
@@ -8220,6 +8246,29 @@ CConfig::~CConfig() {
 
   delete [] nBlades;
   delete [] FreeStreamTurboNormal;
+  //Added by max
+
+  /*--- Free variables used for the VG model ---*/
+
+  for(unsigned short iVG=0;iVG<nVgs;iVG++){
+    for(unsigned short iPoint=0;iPoint<nPointsVG;iPoint++){
+      delete [] Vg_Coordinates[iVG][iPoint];
+    }
+    delete [] Vg_Coordinates[iVG];
+    delete [] Vg_SurfaceNormalDirection[iVG];
+    delete [] Vg_SurfaceTangentialDirection[iVG];
+    delete [] Vg_SurfaceCrossFlowDirection[iVG];
+    delete [] Vg_StreamwiseDirection[iVG];
+  }
+  delete [] Vg_Coordinates;
+  delete [] Vg_SurfaceNormalDirection;
+  delete [] Vg_SurfaceTangentialDirection;
+  delete [] Vg_SurfaceCrossFlowDirection;
+  delete [] Vg_StreamwiseDirection;
+  delete [] Vg_SurfaceArea;
+  delete [] Vg_Angle;
+
+//End added by max
 }
 
 string CConfig::GetFilename(string filename, const string& ext, int timeIter) const {
diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp
index 5b09a9cdd71..8ece513941b 100644
--- a/Common/src/grid_movement/CVolumetricMovement.cpp
+++ b/Common/src/grid_movement/CVolumetricMovement.cpp
@@ -1932,6 +1932,7 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry* geometry, CConfig* config, u
   su2double cosPhi, sinPhi, cosPsi, sinPsi;
   bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);
   bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint());
+  bool bay = config->GetVGModel() != ENUM_VG_MODEL::NONE;
 
   /*--- Problem dimension and physical time step ---*/
   nDim = geometry->GetnDim();
@@ -2061,6 +2062,34 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry* geometry, CConfig* config, u
     config->SetRefOriginMoment_Z(jMarker, Center[2] + rotCoord[2]);
   }
 
+  /*--- Update Vortex Generator Model Variables ---*/
+
+  if (bay) {
+    for (unsigned short iVG = 0; iVG < config->Get_nVGs(); iVG++) {
+      auto* t = config->Get_tVG(iVG);
+      auto* b = config->Get_bVG(iVG);
+      auto* n = config->Get_nVG(iVG);
+      auto* uhat = config->Get_uhatVg(iVG);
+
+      for (iDim = 0; iDim < nDim; iDim++) {
+        t[iDim] = t[0] * rotMatrix[iDim][0] + t[1] * rotMatrix[iDim][1] + t[2] * rotMatrix[iDim][2];
+        b[iDim] = b[0] * rotMatrix[iDim][0] + b[1] * rotMatrix[iDim][1] + b[2] * rotMatrix[iDim][2];
+        n[iDim] = n[0] * rotMatrix[iDim][0] + n[1] * rotMatrix[iDim][1] + n[2] * rotMatrix[iDim][2];
+        uhat[iDim] = uhat[0] * rotMatrix[iDim][0] + uhat[1] * rotMatrix[iDim][1] + uhat[2] * rotMatrix[iDim][2];
+      }
+      auto coords_vg = config->GetVGcoord(iVG);
+      for (unsigned short iPoint = 0; iPoint < config->Get_nPointsVg(); iPoint++) {
+        for (iDim = 0; iDim < nDim; iDim++) r[iDim] = (coords_vg[iPoint][iDim] - Center[iDim]) / Lref;
+        rotCoord[0] = rotMatrix[0][0] * r[0] + rotMatrix[0][1] * r[1] + rotMatrix[0][2] * r[2];
+
+        rotCoord[1] = rotMatrix[1][0] * r[0] + rotMatrix[1][1] * r[1] + rotMatrix[1][2] * r[2];
+
+        rotCoord[2] = rotMatrix[2][0] * r[0] + rotMatrix[2][1] * r[1] + rotMatrix[2][2] * r[2];
+        for (iDim = 0; iDim < nDim; iDim++) coords_vg[iPoint][iDim] = rotCoord[iDim] + Center[iDim];
+      }
+    }
+  }
+
   /*--- After moving all nodes, update geometry class ---*/
 
   UpdateDualGrid(geometry, config);
@@ -2082,6 +2111,7 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry* geometry, CConfig* config, u
   unsigned long iPoint;
   bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);
   bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint());
+  bool bay = config->GetVGModel() != ENUM_VG_MODEL::NONE;
 
   /*--- Retrieve values from the config file ---*/
   deltaT = config->GetDelta_UnstTimeND();
@@ -2205,7 +2235,39 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry* geometry, CConfig* config, u
     }
   }
 
-  /*--- For pitching we don't update the motion origin and moment reference origin. ---*/
+  // Added by max
+
+  /*--- Update Vortex Generator Model Variables ---*/
+
+  if (bay) {
+    for (unsigned short iVG = 0; iVG < config->Get_nVGs(); iVG++) {
+      auto* t = config->Get_tVG(iVG);
+      auto* b = config->Get_bVG(iVG);
+      auto* n = config->Get_nVG(iVG);
+      auto* uhat = config->Get_uhatVg(iVG);
+
+      for (iDim = 0; iDim < nDim; iDim++) {
+        t[iDim] = t[0] * rotMatrix[iDim][0] + t[1] * rotMatrix[iDim][1] + t[2] * rotMatrix[iDim][2];
+        b[iDim] = b[0] * rotMatrix[iDim][0] + b[1] * rotMatrix[iDim][1] + b[2] * rotMatrix[iDim][2];
+        n[iDim] = n[0] * rotMatrix[iDim][0] + n[1] * rotMatrix[iDim][1] + n[2] * rotMatrix[iDim][2];
+        uhat[iDim] = uhat[0] * rotMatrix[iDim][0] + uhat[1] * rotMatrix[iDim][1] + uhat[2] * rotMatrix[iDim][2];
+      }
+      auto coords_vg = config->GetVGcoord(iVG);
+      for (unsigned short iPoint = 0; iPoint < config->Get_nPointsVg(); iPoint++) {
+        for (iDim = 0; iDim < nDim; iDim++) r[iDim] = (coords_vg[iPoint][iDim] - Center[iDim]) / Lref;
+        rotCoord[0] = rotMatrix[0][0] * r[0] + rotMatrix[0][1] * r[1] + rotMatrix[0][2] * r[2];
+
+        rotCoord[1] = rotMatrix[1][0] * r[0] + rotMatrix[1][1] * r[1] + rotMatrix[1][2] * r[2];
+
+        rotCoord[2] = rotMatrix[2][0] * r[0] + rotMatrix[2][1] * r[1] + rotMatrix[2][2] * r[2];
+        for (iDim = 0; iDim < nDim; iDim++) coords_vg[iPoint][iDim] = rotCoord[iDim] + Center[iDim];
+      }
+    }
+  }
+
+  // end added by max
+
+  /*--- For pitching we   don't update the motion origin and moment reference origin. ---*/
 
   /*--- After moving all nodes, update geometry class ---*/
 
@@ -2222,6 +2284,7 @@ void CVolumetricMovement::Rigid_Plunging(CGeometry* geometry, CConfig* config, u
   unsigned long iPoint;
   bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);
   bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint());
+  bool bay = config->GetVGModel() != ENUM_VG_MODEL::NONE;
 
   /*--- Retrieve values from the config file ---*/
   deltaT = config->GetDelta_UnstTimeND();
@@ -2335,6 +2398,17 @@ void CVolumetricMovement::Rigid_Plunging(CGeometry* geometry, CConfig* config, u
     config->SetRefOriginMoment_Z(jMarker, Center[2]);
   }
 
+  /*--- Update Vortex Generator Model Variables ---*/
+
+  if (bay) {
+    for (unsigned short iVG = 0; iVG < config->Get_nVGs(); iVG++) {
+      auto coords_vg = config->GetVGcoord(iVG);
+      for (unsigned short iPoint = 0; iPoint < config->Get_nPointsVg(); iPoint++) {
+        for (iDim = 0; iDim < nDim; iDim++) coords_vg[iPoint][iDim] = coords_vg[iPoint][iDim] + deltaX[iDim];
+      }
+    }
+  }
+
   /*--- After moving all nodes, update geometry class ---*/
 
   UpdateDualGrid(geometry, config);
@@ -2350,6 +2424,7 @@ void CVolumetricMovement::Rigid_Translation(CGeometry* geometry, CConfig* config
   unsigned long iPoint;
   bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);
   bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint());
+  bool bay = config->GetVGModel() != ENUM_VG_MODEL::NONE;
 
   /*--- Retrieve values from the config file ---*/
   deltaT = config->GetDelta_UnstTimeND();
@@ -2448,6 +2523,17 @@ void CVolumetricMovement::Rigid_Translation(CGeometry* geometry, CConfig* config
   /*--- After moving all nodes, update geometry class ---*/
 
   UpdateDualGrid(geometry, config);
+
+  /*--- Update Vortex Generator Model Variables ---*/
+
+  if (bay) {
+    for (unsigned short iVG = 0; iVG < config->Get_nVGs(); iVG++) {
+      auto coords_vg = config->GetVGcoord(iVG);
+      for (unsigned short iPoint = 0; iPoint < config->Get_nPointsVg(); iPoint++) {
+        for (iDim = 0; iDim < nDim; iDim++) coords_vg[iPoint][iDim] = coords_vg[iPoint][iDim] + deltaX[iDim];
+      }
+    }
+  }
 }
 
 void CVolumetricMovement::SetVolume_Scaling(CGeometry* geometry, CConfig* config, bool UpdateGeo) {
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index 450112df8f7..0c9501fec57 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -189,7 +189,10 @@ class CNumerics {
   bool nemo;                      /*!< \brief Flag for NEMO problems  */
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
-
+  //Added by max
+  unsigned long Point_i,Point_j; /*!< \brief Points indexes */
+  unsigned long Edge; /*!< \brief Edge indexe */
+  //End added by max
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -847,6 +850,25 @@ class CNumerics {
     Coord_i = val_coord_i;
     Coord_j = val_coord_j;
   }
+  //Added by max
+  /*!
+   * \brief Set index of the points.
+   * \param[in] val_coord_i - Index of the point i.
+   * \param[in] val_coord_j - Index of the point j.
+   */
+  inline void SetIndex(const unsigned long idxP1, const unsigned long idxP2) {
+    Point_i = idxP1;
+    Point_j = idxP2;
+  }
+
+  /*!
+   * \brief Set index of the edge.
+   * \param[in] idx_edge - Index of the edge.
+   */
+  inline void SetEdge(const unsigned long idx_edge){
+    Edge=idx_edge;
+  }
+  //End added by max
 
   /*!
    * \brief Set the velocity of the computational grid.
@@ -1605,6 +1627,13 @@ class CNumerics {
    * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport
    */
   inline bool GetBoundedScalar() const { return bounded_scalar;}
+  //Added by max
+   /*!
+   * \brief Updat the source container for non-uniform sources
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void UpdateSource(const CConfig* config){};
+  //end added by max
 };
 
 /*!
diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp
index ab9bedc8037..9287afc18b4 100644
--- a/SU2_CFD/include/numerics/flow/flow_sources.hpp
+++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp
@@ -467,3 +467,142 @@ class CSourceRadiation : public CSourceBase_Flow {
   ResidualType<> ComputeResidual(const CConfig* config) override;
 
 };
+
+/*!
+ * \class CSourceBAYModel
+ * \brief Source term for BAY model for Vortex Generators
+ * \ingroup SourceDiscr
+ * \author M. Lagorio
+ */
+class CSourceBAYModel : public CSourceBase_Flow {
+ private:
+  bool implicit;                 /*!< \brief Implicit Calculation */
+  bool first_computation{true}; /*!< \brief Flag for MPI comunication*/
+  su2double calibrationConstant; /*!< \bried Calibration constant for the BAY model force calculation*/
+  unsigned short nVgs{0}; /*!< \brief Number of Vortex Generators */
+
+  /*!
+   *\class Edge_info_VGModel
+   *\brief Data structure to store the edge information for the BAY model
+   */
+  struct Edge_info_VGModel {
+    unsigned long iPoint, jPoint;
+    su2double iDistance, jDistance, cells_volume;
+};
+  /*!
+   *\class Vortex_Generator
+   *\brief Data structure to store the Vortex Generator information for the BAY model
+   */
+  public:
+  class Vortex_Generator {
+  private:
+
+    public:
+    su2double** coords_vg = nullptr; /*!< \brief Data structure to store the points defining the VG */
+    su2double* t=nullptr; /*!< \brief Vectors defining the VG directions */
+    su2double* b=nullptr;
+    su2double* n=nullptr;
+
+    su2double Vtot{0.0};        /*!< \brief Total Volume of the cells defining the Vortex Generator Volume*/
+    unsigned short nPoints;
+
+
+    map<unsigned long, Edge_info_VGModel*>
+        EdgesBay; /*!< \brief Data structure to store edge informations defining the BAY model */
+
+    map<unsigned long, unsigned long> pointsBAY; /*!< \brief Structure to map the mesh nodes to edges */
+
+    su2double Svg; /*!< \brief Platform area of the VG */
+
+    /*!
+     * \brief VG constructor
+     * \param[in] config - Definition of the particular problem
+     * \param[in] iVG - ndex of the vortex generator 
+     */
+    Vortex_Generator(const CConfig* config, const unsigned short iVG);
+
+    /*!
+     * \brief Class deconstructor
+     */
+     ~Vortex_Generator();
+
+    /*!
+     * \brief Add cell volume to total volume of the VG
+     * \param[in] Vcell - Cell volume
+     */
+    void addVGcellVolume(const su2double Vcell) { Vtot += Vcell; }
+
+    /*!
+     * \brief Checks if the point is alredy selected by anothe edge
+     * \param[in] Point - Point to check
+     * \param[in] jEdge - Variable to store the previous edge
+     * \param[in] distanceOld - Variable to store the distance between the Vg and the previous point
+     */
+    bool Check_edge_map(const unsigned long Point, unsigned long &jEdge,su2double &distanceOld);
+
+    /*!
+     * \brief Checks if the edge intersects the vortex generator
+     * \param[in] Coord_i - Coordinates of the first point
+     * \param[in] Coord_j - Coodinates of the second point
+     * \param[in] Normal - Normal vector
+     * \param[out] distanceToVg - Distance between the edge and the VG
+     */
+    bool EdgeIntersectsVG(const su2double *Coord_i,const su2double *Coord_j, const su2double *Normal,su2double &distanceToVg);
+
+    /*!
+     * \brief Method returns if a point is inside the VG
+     * \param[in] vgModel - Type of VG model
+     * \param[in] Point_i - Index of the point
+     * \param[in] Edge - Index of the edge
+     * \return True if the point is inside the VG
+     */
+    const bool PointInVG(const ENUM_VG_MODEL vgModel,const unsigned long& Point_i,const unsigned long& Edge);
+
+    /*!
+     * \brief Method returns interpolation coefficient for obtaining velocity at VG plane
+     * \param[in] vgModel - Type of VG model
+     * \param[in] Edge - Index of the edge
+     * \param[in] Point_i - Index of the point
+     * \return Interpolation coefficient
+     */
+    const su2double Get_interpolationCoefficient(ENUM_VG_MODEL vgModel,const unsigned long Edge, const unsigned long Point_i);
+
+    /*!
+     * \brief Method returns jBAY redistribution constant
+     * \param[in] vgModel - Type of VG model
+     * \param[in] Edge - Index of the edge
+     * \param[in] Point_i - Index of the point
+     * \param[in] Volume - Volume of the cells
+     * \return Interpolation coefficient
+     */
+    const su2double Get_redistributionConstant(const ENUM_VG_MODEL vgModel,const unsigned long Edge, const unsigned long Point_i, const su2double Volume);
+  };
+
+  vector<Vortex_Generator*> VGs{nullptr}; /*!< \brief Vector storing all VGs data structures */
+
+
+public:
+  CSourceBAYModel(unsigned short val_ndim, unsigned short val_nVar, const CConfig* config);
+
+  ~CSourceBAYModel() {
+    for (auto* iVg : VGs) {
+      delete iVg;
+    }
+    VGs.clear();
+  }
+
+  /*!
+   * \brief Source term integration for a VG Bay model source.
+   * \param[in] config - Definition of the particular problem.
+   * \return Lightweight const-view of residual and Jacobian.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override;
+
+  void ReduceVariables(void);
+
+  /*!
+   * \brief Source term initialization
+   */
+  void UpdateSource(const CConfig* config) override;
+
+};
diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp
index a6925c28e94..0e016d855a6 100644
--- a/SU2_CFD/include/output/CFlowOutput.hpp
+++ b/SU2_CFD/include/output/CFlowOutput.hpp
@@ -335,6 +335,14 @@ class CFlowOutput : public CFVMOutput{
    */
   void WriteForcesBreakdown(const CConfig *config, const CSolver *flow_solver) const;
 
+//Added by max
+  /*!
+   * \brief Write restart file for the vortex generator models
+   * \param[in] config - Definition of the particular problem per zone.
+   */
+  void WriteVGModelRestartConfig(const CConfig* config) const;
+
+//End added by max
   /*!
    * \brief Set the time averaged output fields.
    */
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
index 5e0e34b7a3e..85b85385f1b 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
@@ -2432,4 +2432,14 @@ class CFVMFlowSolverBase : public CSolver {
    */
   inline const su2activevector* GetEdgeMassFluxes() const final { return &EdgeMassFluxes; }
 
+  //Added by Max
+  /*!
+   * \brief Preprocess source terms.
+   * \param[in] geometry - Geometrical definition of the particular problem
+   * \param[in] numeric - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void PreprocessSources(CGeometry* geometry,CNumerics** numerics,CConfig* config);
+//end added by max
+
 };
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index c3051311c30..671d5173155 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -3049,3 +3049,37 @@ void CFVMFlowSolverBase<V, FlowRegime>::ComputeAxisymmetricAuxGradients(CGeometr
     SetAuxVar_Gradient_LS(geometry, config);
   }
 }
+
+//Added by max
+template <class V, ENUM_REGIME FlowRegime>
+void CFVMFlowSolverBase<V, FlowRegime>::PreprocessSources(CGeometry* geometry, CNumerics** numerics, CConfig* config){
+  
+  const bool bay = config->GetVGModel() != ENUM_VG_MODEL::NONE;
+  
+  /*Preprocess BAY model*/
+  if (bay) {
+    CNumerics* second_numerics = numerics[SOURCE_SECOND_TERM+omp_get_thread_num() * MAX_TERMS];
+    unsigned long iPoint, jPoint;
+    AD::StartNoSharedReading();
+
+    /* Loop over the edges */
+    
+    for (auto color : EdgeColoring) {
+      SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize))
+      for (auto k = 0ul; k < color.size; ++k) {
+        auto iEdge = color.indices[k];
+        iPoint = geometry->edges->GetNode(iEdge, 0);
+        jPoint = geometry->edges->GetNode(iEdge, 1);
+        second_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint));
+        second_numerics->SetNormal(geometry->edges->GetNormal(iEdge));
+        second_numerics->SetIndex(geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetGlobalIndex(jPoint));
+        second_numerics->SetVolume(geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetVolume(jPoint));
+        second_numerics->SetEdge(iEdge);
+        second_numerics->UpdateSource(config);
+      }
+      END_SU2_OMP_FOR
+    }
+    AD::EndNoSharedReading();
+  }
+}
+//end added by max
\ No newline at end of file
diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
index 6f96628ddc1..3a2b6668cbd 100644
--- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp
+++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
@@ -421,4 +421,5 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
    * \param[in] config - The particular config.
    */
   void ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) final;
+
 };
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index 2b114a7714b..b8375130687 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -171,6 +171,8 @@ class CSolver {
    */
   void SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config);
 
+   
+
 private:
 
   /*!
@@ -1578,6 +1580,12 @@ class CSolver {
                                       CNumerics *numerics,
                                       CConfig *config,
                                       unsigned short iMesh) { }
+  //Added by max
+  void ReadVGConfigFile(CConfig* config);
+
+  void InitializeVGVariables(const unsigned short nVgs_file, std::vector<std::string>& lines_configVg, CConfig* config);
+
+  // End added by max
 
   /*!
    * \brief A virtual member.
@@ -4326,6 +4334,9 @@ class CSolver {
     }
     END_SU2_OMP_FOR
   }
+    // Added by max
+  virtual void PreprocessSources(CGeometry* geometry,CNumerics** numerics,CConfig* config){};
+  // End added by max
 
 protected:
   /*!
@@ -4426,4 +4437,5 @@ class CSolver {
     }
   }
 
+
 };
diff --git a/SU2_CFD/include/variables/CFlowVariable.hpp b/SU2_CFD/include/variables/CFlowVariable.hpp
index 5a28c16ea6b..4eec3ece26c 100644
--- a/SU2_CFD/include/variables/CFlowVariable.hpp
+++ b/SU2_CFD/include/variables/CFlowVariable.hpp
@@ -51,6 +51,9 @@ class CFlowVariable : public CVariable {
    * ---*/
   MatrixType Vorticity; /*!< \brief Vorticity of the flow field. */
   VectorType StrainMag; /*!< \brief Magnitude of rate of strain tensor. */
+   //Added by max
+  VectorType VG_Locations; /*!< \brief Cell locations where the Vg model is applied */
+  //end Added by max
 
   /*!
    * \brief Constructor of the class.
@@ -267,4 +270,18 @@ class CFlowVariable : public CVariable {
    * \return Vector of magnitudes.
    */
   inline su2activevector& GetStrainMag() { return StrainMag; }
+  
+  /*!
+   * \brief Get status of the vortex generator model per node
+   * \param[in] iPoint - Node index
+   * \return 0 if the vg model is not active in the cell, 1 if the model is active
+   */
+  inline su2double Get_VGLocations(unsigned long iPoint) const  override {return VG_Locations(iPoint);}
+
+  /*!
+   * \brief Set status of the vortex generator model per node
+   * \param[in] iPoint - Node index
+   * \param[in] bool_vgLoc - status of the vg model
+   */
+  inline void Set_VGLocations(unsigned long iPoint, su2double bool_vgLoc) override  {VG_Locations(iPoint)=bool_vgLoc;}
 };
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index 358824ae8d7..cac61bbec5e 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -2327,6 +2327,22 @@ class CVariable {
    */
   inline virtual unsigned long GetNewtonSolverIterations(unsigned long iPoint) const { return 0; }
 
+  //Added by max DEBUG REMOVE
+   /*!
+   * \brief Get status of the vortex generator model per node
+   * \param[in] iPoint - Node index
+   * \return 0 if the vg model is not active in the cell, 1 if the model is active
+   */
+  virtual inline su2double Get_VGLocations(unsigned long iPoint) const  {return 0;}
+
+  /*!
+   * \brief Set status of the vortex generator model per node
+   * \param[in] iPoint - Node index
+   * \param[in] bool_vgLoc - status of the vg model
+   */
+  virtual inline void Set_VGLocations(unsigned long iPoint, su2double bool_vgLoc)  {}
+ //end added by max
+
   /*!
    * \brief LUT premixed flamelet: virtual functions for the speciesflameletvariable LUT
    */
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 7160cc93e6f..355ef0e4349 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -1039,6 +1039,10 @@ void CDriver::InitializeSolver(CConfig* config, CGeometry** geometry, CSolver **
 
   PreprocessInlet(solver, geometry, config);
 
+  //Added by max
+  if(config->GetVGModel()!=ENUM_VG_MODEL::NONE) solver[MESH_0][FLOW_SOL]->ReadVGConfigFile(config);
+  //End added by max
+
 }
 
 void CDriver::PreprocessInlet(CSolver ***solver, CGeometry **geometry, CConfig *config) const {
@@ -1889,6 +1893,12 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver
         /*--- Vorticity Confinement term as a second source term to allow the use of VC along with other source terms ---*/
         numerics[iMGlevel][FLOW_SOL][source_second_term] = new CSourceVorticityConfinement(nDim, nVar_Flow, config);
       }
+      //Added by Max
+      else if (config->GetVGModel()!=ENUM_VG_MODEL::NONE){
+        numerics[iMGlevel][FLOW_SOL][source_second_term] = new CSourceBAYModel(nDim, nVar_Flow, config);
+        solver[iMGlevel][FLOW_SOL]->PreprocessSources(geometry[iMGlevel],numerics[iMGlevel][FLOW_SOL],config);
+      }
+      //End added by Max
       else
         numerics[iMGlevel][FLOW_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Flow, config);
 
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index 9ed1305f3ce..a10747b8fb8 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -837,3 +837,237 @@ CNumerics::ResidualType<> CSourceRadiation::ComputeResidual(const CConfig *confi
 
   return ResidualType<>(residual, jacobian, nullptr);
 }
+
+//Added by Max
+
+CSourceBAYModel::CSourceBAYModel(unsigned short val_ndim, unsigned short val_nVar, const CConfig* config)
+    : CSourceBase_Flow(val_ndim, val_nVar, config) {
+
+  if (val_ndim != 3) SU2_MPI::Error("BAY model can be used only in 3D", CURRENT_FUNCTION);
+
+  implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+
+  nVgs = config->Get_nVGs();
+
+  VGs.resize(nVgs);
+
+  calibrationConstant = config->GetVGConstant();
+
+  for (unsigned short iVG = 0; iVG < nVgs; iVG++) {
+    VGs[iVG] = new Vortex_Generator(config, iVG);
+  };
+};
+
+CNumerics::ResidualType<> CSourceBAYModel::ComputeResidual(const CConfig* config) {
+
+/*Calculate total volume across ranks*/
+#if HAVE_MPI
+  if (first_computation) {
+    ReduceVariables();
+  }
+#endif
+
+  /*Zero the jacobian and residual vector*/
+  residual[0] = 0.0;
+  for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+    residual[iVar] = 0.0;
+    for (unsigned short jVar = 0; jVar < nVar; jVar++) {
+      jacobian[iVar][jVar] = 0.0;
+    }
+  }
+
+  /*Iterate over VGs and check if the edge is part of the VG domain*/
+
+  for (auto* iVG : VGs) {
+    if (iVG->PointInVG(config->GetVGModel(), Point_i, Edge)){
+      const auto S = iVG->Svg;
+      const auto b = iVG->b;
+      const auto n = iVG->n;
+      const auto Vtot = iVG->Vtot;
+      const auto t = iVG->t;
+
+      const auto rho = V_i[nDim+2];
+
+      const auto interpolation_coeff = iVG->Get_interpolationCoefficient(config->GetVGModel(), Edge, Point_i);
+      const auto redistribution_const=iVG->Get_redistributionConstant(config->GetVGModel(), Edge, Point_i, Volume);
+
+      /*Compute velocity at the vg surface*/
+
+      const su2double u[3]{interpolation_coeff * V_i[1] + (1 - interpolation_coeff) * V_j[1],
+                           interpolation_coeff * V_i[2] + (1 - interpolation_coeff) * V_j[2],
+                           interpolation_coeff * V_i[3] + (1 - interpolation_coeff) * V_j[3]};
+
+      /*Compute Residual*/
+
+      su2double momentumResidual[3];
+
+      su2double un = GeometryToolbox::DotProduct(nDim, u, n);
+      su2double ut = GeometryToolbox::DotProduct(nDim, u, t);
+
+      su2double k = (calibrationConstant * S * Volume / Vtot) * rho * un * ut / GeometryToolbox::Norm(nDim, u);
+      GeometryToolbox::CrossProduct(u, b, momentumResidual);
+      if (pow(GeometryToolbox::Norm(nDim, u), 3) > EPS && fabs(k) > EPS) {
+        residual[1] = redistribution_const * k * momentumResidual[0];
+        residual[2] = redistribution_const * k * momentumResidual[1];
+        residual[3] = redistribution_const * k * momentumResidual[2];
+
+        if (implicit) {
+          /*Compute Jacobian*/
+
+          /* Cross product contribution */
+          jacobian[1][2] = b[2];
+          jacobian[1][3] = -b[1];
+
+          jacobian[2][1] = -b[2];
+          jacobian[2][3] = b[0];
+
+          jacobian[3][1] = b[1];
+          jacobian[3][2] = -b[0];
+
+            /* Dot product contribution */
+            for (unsigned short iVar = 1; iVar < nDim + 1; iVar++) {
+            for (unsigned short jVar = 1; jVar < nDim + 1; jVar++) {
+              jacobian[iVar][jVar] += n[jVar - 1] * momentumResidual[iVar - 1] * ut +
+                          un * momentumResidual[iVar - 1] * t[jVar - 1] +
+                          un * momentumResidual[iVar - 1] * ut * (-u[jVar - 1] / GeometryToolbox::SquaredNorm(nDim, u));
+              jacobian[iVar][jVar] *= ((redistribution_const * calibrationConstant * S * Volume / Vtot) /
+                           GeometryToolbox::Norm(nDim, u) / rho);
+            }
+            }
+        }
+      }
+    }
+  }
+
+  return ResidualType<>(residual, jacobian, nullptr);
+}
+
+CSourceBAYModel::Vortex_Generator::Vortex_Generator(const CConfig* config, unsigned short iVG) {
+
+  b=config->Get_bVG(iVG);
+  n=config->Get_nVG(iVG);
+  t=config->Get_tVG(iVG);
+
+  coords_vg=config->GetVGcoord(iVG);
+  nPoints = config->Get_nPointsVg();
+
+  Svg=config->Get_Svg(iVG);
+
+};
+
+CSourceBAYModel::Vortex_Generator::~Vortex_Generator(){
+
+  for(auto itr = EdgesBay.begin(); itr != EdgesBay.end(); itr++) {
+    delete itr->second;
+  }
+}
+
+const bool CSourceBAYModel::Vortex_Generator::PointInVG(const ENUM_VG_MODEL vgModel,const unsigned long& Point_i,const unsigned long& Edge){
+  switch (vgModel) {
+          case ENUM_VG_MODEL::BAY:
+            return pointsBAY.find(Point_i)!= pointsBAY.end();
+          case ENUM_VG_MODEL::JBAY:
+            return EdgesBay.find(Edge)!= EdgesBay.end();
+          default:
+            return false;
+            break; //avoid compiler warning
+        }
+};
+
+const su2double CSourceBAYModel::Vortex_Generator::Get_interpolationCoefficient(const ENUM_VG_MODEL vgModel,const unsigned long Edge, const unsigned long Point_i){
+      if(vgModel==ENUM_VG_MODEL::JBAY){
+          const auto edgeInfo = *(EdgesBay.find(Edge)->second);
+          if (Point_i == edgeInfo.iPoint)
+            return edgeInfo.iDistance / (edgeInfo.iDistance + edgeInfo.jDistance);
+          else
+            return edgeInfo.jDistance / (edgeInfo.iDistance + edgeInfo.jDistance);
+      }
+      else return 1.0;
+};
+
+const su2double CSourceBAYModel::Vortex_Generator::Get_redistributionConstant(const ENUM_VG_MODEL vgModel,const unsigned long Edge, const unsigned long Point_i, const su2double Volume){
+  if(vgModel==ENUM_VG_MODEL::JBAY){
+          auto edgeInfo = *(EdgesBay.find(Edge)->second);
+          if (Point_i == edgeInfo.iPoint)
+            return edgeInfo.cells_volume / (Volume + (edgeInfo.iDistance / edgeInfo.jDistance) * (edgeInfo.cells_volume - Volume));
+          else
+            return edgeInfo.cells_volume / (Volume + (edgeInfo.jDistance / edgeInfo.iDistance) * (edgeInfo.cells_volume - Volume));
+      }
+      else return 1.0;
+}
+
+void CSourceBAYModel::UpdateSource(const CConfig* config) {
+
+  su2double iDistance;
+
+  /*Check if edge intersects a VG*/
+
+  for (auto* iVG : VGs) {
+
+    if (iVG->EdgeIntersectsVG(Coord_i,Coord_j,Normal,iDistance)) {
+
+        const su2double pointDistance = GeometryToolbox::Distance(3, Coord_i, Coord_j);
+        su2double distanceOld;
+        unsigned long jEdge;
+
+        auto alreadySelected = iVG->Check_edge_map(Point_i,jEdge,distanceOld)||iVG->Check_edge_map(Point_j,jEdge,distanceOld);
+
+        if (alreadySelected && distanceOld > pointDistance) {
+          iVG->pointsBAY.erase(iVG->EdgesBay[jEdge]->iPoint);
+          iVG->pointsBAY.erase(iVG->EdgesBay[jEdge]->jPoint);
+          iVG->addVGcellVolume(-iVG->EdgesBay[jEdge]->cells_volume);  // remove old volume from total volume
+          delete iVG->EdgesBay.find(jEdge)->second;
+          iVG->EdgesBay.erase(jEdge);
+          alreadySelected=false;
+        }
+
+        if(!alreadySelected){
+          iVG->pointsBAY.insert(make_pair(Point_i, Edge));
+          iVG->pointsBAY.insert(make_pair(Point_j, Edge));
+          Edge_info_VGModel* edge_info = new Edge_info_VGModel{Point_i, Point_j, iDistance, pointDistance - iDistance, Volume};
+          iVG->EdgesBay.insert(make_pair(Edge, edge_info));
+          iVG->addVGcellVolume(Volume);
+        }
+      }
+    }
+  }
+
+
+bool CSourceBAYModel::Vortex_Generator::Check_edge_map(const unsigned long Point, unsigned long &jEdge,su2double &distanceOld) {
+  if (this->pointsBAY.find(Point) != this->pointsBAY.end()) {
+    jEdge = this->pointsBAY[Point];
+    auto edgeOld = this->EdgesBay.find(jEdge)->second;
+    distanceOld = edgeOld->iDistance + edgeOld->jDistance;
+    return true;
+  } else
+    return false;
+}
+
+bool CSourceBAYModel::Vortex_Generator::EdgeIntersectsVG(const su2double* Coord_i,
+                                                         const su2double* Coord_j, const su2double* Normal, su2double& distanceToVg) {
+
+  su2double vg_edge_intersection[3];
+  if (GeometryToolbox::IntersectEdge(3, coords_vg[0], n, Coord_i, Coord_j)) {
+    distanceToVg = GeometryToolbox::LinePlaneIntersection<su2double, 3>(
+        Coord_i, Normal, coords_vg[0], n, vg_edge_intersection);
+
+    if (GeometryToolbox::PointInConvexPolygon(3, coords_vg, vg_edge_intersection, 4)) return true;
+  }
+  return false;
+}
+
+void CSourceBAYModel::ReduceVariables(void){
+  for (auto iVG : VGs) {
+      su2double V_tot_glob;
+      su2double V_loc{iVG->Vtot};
+
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+      SU2_MPI::Allreduce(&V_loc, &V_tot_glob, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
+
+      iVG->Vtot = V_tot_glob;
+      first_computation = false;
+    }
+}
+
+//End added by Max
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 5ffee0d23f5..7c295ea0f77 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1489,6 +1489,9 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
     AddVolumeOutput("TURB_DELTA_TIME", "Turb_Delta_Time", "TIMESTEP", "Value of the local timestep for the turbulence variables");
     AddVolumeOutput("TURB_CFL", "Turb_CFL", "TIMESTEP", "Value of the local CFL for the turbulence variables");
   }
+   //Added by max DEBUG remove
+  AddVolumeOutput("VG_LOC", "VG_Locations", "VG_MODEL", "Debug variable for vg locations");
+  //end Added by max
 }
 
 void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry,
@@ -1513,6 +1516,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     }
     SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint)));
   }
+  
 
   const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE);
 
@@ -1568,7 +1572,12 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint));
     SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
   }
+  //Added by max DEBUG REMOVE
+  if(config->GetVGModel()!=ENUM_VG_MODEL::NONE){
+    SetVolumeOutputValue("VG_LOC",iPoint, Node_Flow->Get_VGLocations(iPoint));
+  }
 
+//End added by Max
   switch (config->GetKind_Species_Model()) {
 
     case SPECIES_MODEL::SPECIES_TRANSPORT: {
@@ -2355,6 +2364,9 @@ void CFlowOutput::WriteAdditionalFiles(CConfig *config, CGeometry *geometry, CSo
     WriteForcesBreakdown(config, solver_container[FLOW_SOL]);
   }
 
+  if(config->GetVGModel()!=ENUM_VG_MODEL::NONE && config->GetTime_Domain()){
+    WriteVGModelRestartConfig(config);
+  }
 }
 
 void CFlowOutput::WriteMetaData(const CConfig *config){
@@ -3915,6 +3927,45 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo
   // clang-format on
 }
 
+//Added by max
+
+void CFlowOutput::WriteVGModelRestartConfig(const CConfig* config) const {
+  if(rank!=MASTER_NODE) return;
+  
+  auto filename = config->GetUnsteady_FileName(config->GetVGFileName(),curTimeIter,".cfg");
+  ofstream file;
+  file.open(filename);
+  
+  auto nVgs = config->Get_nVGs();
+  su2double l,h1,h2,beta;
+  su2double **vg_coord,*uhat,*unhat;
+  unsigned short iDim;
+
+  file.precision(5);
+
+  for(unsigned short iVG=0; iVG<nVgs;iVG++){
+
+    vg_coord = config->GetVGcoord(iVG);
+    beta =config->Get_betaVg(iVG);
+    uhat = config->Get_uhatVg(iVG);
+    unhat=config->Get_bVG(iVG);
+
+    l = GeometryToolbox::Distance(nDim,vg_coord[0],vg_coord[1]);
+    h1 = GeometryToolbox::Distance(nDim,vg_coord[0],vg_coord[3]);
+    h2 = GeometryToolbox::Distance(nDim,vg_coord[1],vg_coord[2]);
+
+    file<<l<<" "<<h1<<" "<<h2<<" "<<beta<<" ";
+    
+    for(iDim=0;iDim<nDim;iDim++) file<<vg_coord[0][iDim]<<" ";
+    for(iDim=0;iDim<nDim;iDim++) file<<unhat[iDim]<<" ";
+    for(iDim=0;iDim<nDim;iDim++) file<<uhat[iDim]<<" ";
+    file<<'\n';
+  }
+
+  file.close();
+
+}
+//End added by max
 bool CFlowOutput::WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile){
 
   bool writeRestart = false;
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index 16cdfe59a4b..120a63b3cdd 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -2040,6 +2040,8 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
   const bool ideal_gas        = (config->GetKind_FluidModel() == STANDARD_AIR) ||
                                 (config->GetKind_FluidModel() == IDEAL_GAS);
   const bool rans             = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
+  const bool bay = config->GetVGModel() == ENUM_VG_MODEL::BAY;
+  const bool jbay = config->GetVGModel() == ENUM_VG_MODEL::JBAY;
 
   /*--- Pick one numerics object per thread. ---*/
   CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS];
@@ -2239,6 +2241,76 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
     END_SU2_OMP_FOR
 
   }
+  // Added by Max
+
+  if (jbay) {
+    unsigned long jPoint;
+    CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num() * MAX_TERMS];
+    AD::StartNoSharedReading();
+
+    /* Loop over the edges for the jBAY model */
+    for (auto color : EdgeColoring) {
+      SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize))
+      for (auto k = 0ul; k < color.size; ++k) {
+        auto iEdge = color.indices[k];
+
+        /* Get the points defining the edge */
+        iPoint = geometry->edges->GetNode(iEdge, 0);
+        jPoint = geometry->edges->GetNode(iEdge, 1);
+
+        /* Set edge index */
+        second_numerics->SetEdge(iEdge);
+
+        /* Set variables at the edge extremes*/
+        second_numerics->SetDensity(nodes->GetDensity(iPoint), nodes->GetDensity(iPoint));
+        second_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint));
+
+        for (unsigned short kPoint = 0; kPoint < geometry->edges->GetnNodes(); kPoint++) {
+          const auto Point = geometry->edges->GetNode(iEdge, kPoint);
+          /* Set cell volume */
+          second_numerics->SetVolume(geometry->nodes->GetVolume(Point));
+          /* Set point global index */
+          second_numerics->SetIndex(geometry->nodes->GetGlobalIndex(Point), geometry->nodes->GetGlobalIndex(Point));
+          /* Compute residual */
+          auto residual = second_numerics->ComputeResidual(config);
+          /* Set VG_Loc flag*/
+          if (residual.residual[1] != 0.0 || residual.residual[2] != 0.0 || residual.residual[3] != 0.0)
+            nodes->Set_VGLocations(Point, 1.0);
+          LinSysRes.AddBlock(Point, residual);
+          if (implicit) Jacobian.AddBlock2Diag(Point, residual.jacobian_i);
+        }
+      }
+      END_SU2_OMP_FOR
+    }
+    AD::EndNoSharedReading();
+  }
+
+  if (bay) {
+    AD::StartNoSharedReading();
+    CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num() * MAX_TERMS];
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    /* Loop over the points */
+    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+      /* Set index */
+      second_numerics->SetIndex(geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetGlobalIndex(iPoint));
+      /* Set density*/
+      second_numerics->SetDensity(nodes->GetDensity(iPoint), nodes->GetDensity(iPoint));
+      /* Set primitive */
+      second_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));
+      /* Set Volume */
+      second_numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+      /* Compute residual */
+      auto residual = second_numerics->ComputeResidual(config);
+      /* Set VG_Loc flag*/
+      if (residual.residual[1] != 0.0 || residual.residual[2] != 0.0 || residual.residual[3] != 0.0)
+        nodes->Set_VGLocations(iPoint, 1.0);
+      LinSysRes.AddBlock(iPoint, residual);
+      if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+    }
+    END_SU2_OMP_FOR
+    AD::EndNoSharedReading();
+  }
+  // End added by max
 
   /*--- Check if a verification solution is to be computed. ---*/
 
@@ -9529,4 +9601,4 @@ void CEulerSolver::GatherInOutAverageValues(CConfig *config, CGeometry *geometry
       NuOut[markerTP -1][iSpan]                  = nuOut;
     }
   }
-}
+}
\ No newline at end of file
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 5b81394bd9a..b6ecddb0a0b 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -1392,7 +1392,8 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
   const bool energy         = config->GetEnergy_Equation();
   const bool streamwise_periodic             = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE);
   const bool streamwise_periodic_temperature = config->GetStreamwise_Periodic_Temperature();
-
+  const bool bay = config->GetVGModel() == ENUM_VG_MODEL::BAY;
+  const bool jbay = config->GetVGModel() == ENUM_VG_MODEL::JBAY;
   AD::StartNoSharedReading();
 
   if (body_force) {
@@ -1769,6 +1770,76 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
       }// for iMarker
     }// if !streamwise_periodic_temperature
   }// if streamwise_periodic
+  // Added by Max
+
+  if (jbay) {
+    unsigned long jPoint;
+    CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num() * MAX_TERMS];
+    AD::StartNoSharedReading();
+
+    /* Loop over the edges for the jBAY model */
+    for (auto color : EdgeColoring) {
+      SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize))
+      for (auto k = 0ul; k < color.size; ++k) {
+        auto iEdge = color.indices[k];
+
+        /* Get the points defining the edge */
+        iPoint = geometry->edges->GetNode(iEdge, 0);
+        jPoint = geometry->edges->GetNode(iEdge, 1);
+
+        /* Set edge index */
+        second_numerics->SetEdge(iEdge);
+
+        /* Set variables at the edge extremes*/
+        second_numerics->SetDensity(nodes->GetDensity(iPoint), nodes->GetDensity(iPoint));
+        second_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint));
+
+        for (unsigned short kPoint = 0; kPoint < geometry->edges->GetnNodes(); kPoint++) {
+          const auto Point = geometry->edges->GetNode(iEdge, kPoint);
+          /* Set cell volume */
+          second_numerics->SetVolume(geometry->nodes->GetVolume(Point));
+          /* Set point global index */
+          second_numerics->SetIndex(geometry->nodes->GetGlobalIndex(Point), geometry->nodes->GetGlobalIndex(Point));
+          /* Compute residual */
+          auto residual = second_numerics->ComputeResidual(config);
+          /* Set VG_Loc flag*/
+          if (residual.residual[1] != 0.0 || residual.residual[2] != 0.0 || residual.residual[3] != 0.0)
+            nodes->Set_VGLocations(Point, 1.0);
+          LinSysRes.AddBlock(Point, residual);
+          if (implicit) Jacobian.AddBlock2Diag(Point, residual.jacobian_i);
+        }
+      }
+      END_SU2_OMP_FOR
+    }
+    AD::EndNoSharedReading();
+  }
+
+  if (bay) {
+    AD::StartNoSharedReading();
+    CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num() * MAX_TERMS];
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    /* Loop over the points */
+    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+      /* Set index */
+      second_numerics->SetIndex(geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetGlobalIndex(iPoint));
+      /* Set density*/
+      second_numerics->SetDensity(nodes->GetDensity(iPoint), nodes->GetDensity(iPoint));
+      /* Set primitive */
+      second_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));
+      /* Set Volume */
+      second_numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+      /* Compute residual */
+      auto residual = second_numerics->ComputeResidual(config);
+      /* Set VG_Loc flag*/
+      if (residual.residual[1] != 0.0 || residual.residual[2] != 0.0 || residual.residual[3] != 0.0)
+        nodes->Set_VGLocations(iPoint, 1.0);
+      LinSysRes.AddBlock(iPoint, residual);
+      if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+    }
+    END_SU2_OMP_FOR
+    AD::EndNoSharedReading();
+  }
+  // End added by max
 
   /*--- Check if a verification solution is to be computed. ---*/
 
@@ -3199,4 +3270,4 @@ void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, con
   if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
     adj_sol[0] = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop);
   }
-}
+}
\ No newline at end of file
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 03c79b46fc3..ede6fef29c8 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -4312,3 +4312,139 @@ void CSolver::SavelibROM(CGeometry *geometry, CConfig *config, bool converged) {
 #endif
 
 }
+
+//Added by max
+void CSolver::ReadVGConfigFile(CConfig* config) {
+  string filename, line;
+  unsigned short nVgs = 0;
+  vector<string> lines_vgConfig;
+  bool unsteady_restart = config->GetTime_Domain() && config->GetRestart();
+  bool ts_2nd = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND;
+
+  /*--- Get the vortex generator filename. ---*/
+
+  if (unsteady_restart) {
+
+    int restart_iter = config->GetRestart_Iter()-1-ts_2nd;
+    filename = config->GetUnsteady_FileName(config->GetVGFileName(), restart_iter, ".cfg");
+
+  } else {
+    filename = config->GetVGFileName();
+  }
+
+  /*--- Read the configuation file ---*/
+
+  ifstream file;
+  file.open(filename);
+
+  if (!file.is_open()) {
+    SU2_MPI::Error("Unable to open the vortex generator configuration file "+filename, CURRENT_FUNCTION);
+  }
+
+  while (std::getline(file, line)) {
+    if (line.size()<13 || line.front() == '#') continue;
+    nVgs++;
+    lines_vgConfig.push_back(line);
+  };
+  file.close();
+
+  InitializeVGVariables(nVgs, lines_vgConfig, config);
+}
+
+void CSolver::InitializeVGVariables(const unsigned short nVgs, std::vector<std::string> &lines_vgConfig, CConfig * config)
+{
+  su2double l,h1,h2,beta,p1[3],un[3],u[3],uc[3];
+
+  /*--- Allocate and compute the required varibales for the vortex generator model  ---*/
+
+  su2double **vg_b = new su2double*[nVgs];
+  su2double **vg_n = new su2double*[nVgs];
+  su2double **vg_t = new su2double*[nVgs];
+  su2double **vg_uhat = new su2double*[nVgs];
+  su2double ***vg_coord = new su2double**[nVgs];
+  su2double *Svg = new su2double[nVgs];
+  su2double *betaVg = new su2double[nVgs];
+
+  for (unsigned short iVG = 0; iVG < nVgs; iVG++) {
+
+    /* --- Parse variables ---*/
+    istringstream vg_line{lines_vgConfig[iVG]};
+
+    vg_line >> l;
+    vg_line >> h1;
+    vg_line >> h2;
+    vg_line >> beta;
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      vg_line >> p1[iDim];
+    }
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      vg_line >> un[iDim];
+    }
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      vg_line >> u[iDim];
+    }
+
+    /*--- Allocate Variable ---*/
+
+    vg_b[iVG] = new su2double[3];
+    vg_uhat[iVG] = new su2double[3];
+    vg_n[iVG] = new su2double[3];
+    vg_t[iVG] = new su2double[3];
+
+    /*--- Compute the VG variables ---*/
+
+    betaVg[iVG] = beta;
+    beta *= PI_NUMBER / 180.0;
+
+    GeometryToolbox::CrossProduct(un, u, uc);
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      vg_t[iVG][iDim] = u[iDim] * cos(beta) + uc[iDim] * sin(beta);
+      // vg_n[iVG][iDim] = -u[iDim] * sin(beta) - uc[iDim] * cos(beta);
+      vg_b[iVG][iDim] = un[iDim];
+      vg_uhat[iVG][iDim] = u[iDim];
+    }
+
+    GeometryToolbox::CrossProduct(vg_t[iVG],vg_b[iVG],vg_n[iVG]);
+
+    GeometryToolbox::NormalizeVector(nDim, vg_t[iVG]);
+    GeometryToolbox::NormalizeVector(nDim, vg_b[iVG]);
+    GeometryToolbox::NormalizeVector(nDim, vg_n[iVG]);
+
+    /*--- Set variables in CConfig  ---*/
+
+    vg_coord[iVG] = new su2double*[4];
+
+    vg_coord[iVG][0] = new su2double[3]{p1[0], p1[1], p1[2]};
+
+    vg_coord[iVG][1] = new su2double[3]{p1[0] + vg_t[iVG][0] * l,
+                                        p1[1] + vg_t[iVG][1] * l,
+                                        p1[2] + vg_t[iVG][2] * l};
+
+    su2double* p2=vg_coord[iVG][1];
+
+    vg_coord[iVG][2] = new su2double[3]{p2[0] + vg_b[iVG][0] * h2,
+                                        p2[1] + vg_b[iVG][1] * h2,
+                                        p2[2] + vg_b[iVG][2] * h2};
+
+    vg_coord[iVG][3] = new su2double[3]{p1[0] + vg_b[iVG][0] * h1,
+                                        p1[1] + vg_b[iVG][1] * h1,
+                                        p1[2] + vg_b[iVG][2] * h1};
+
+    Svg[iVG] = 0.5 * (h1 + h2) * l;
+  };
+
+  /*--- Set the variables in CConfig  ---*/
+
+  config->Set_nVG(vg_n);
+  config->Set_bVG(vg_b);
+  config->Set_tVG(vg_t);
+  config->SetVGCoord(vg_coord);
+  config->Set_Svg(Svg);
+  config->Set_nVGs(nVgs);
+  config->Set_uhatVg(vg_uhat);
+  config->Set_betaVg(betaVg);
+}//End added by max
\ No newline at end of file
diff --git a/SU2_CFD/src/variables/CFlowVariable.cpp b/SU2_CFD/src/variables/CFlowVariable.cpp
index c094bc4e043..154556b6ec6 100644
--- a/SU2_CFD/src/variables/CFlowVariable.cpp
+++ b/SU2_CFD/src/variables/CFlowVariable.cpp
@@ -97,6 +97,11 @@ CFlowVariable::CFlowVariable(unsigned long npoint, unsigned long ndim, unsigned
   if (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE) {
     HB_Source.resize(nPoint, nVar) = su2double(0.0);
   }
+
+  //Added by max
+  /* Initialize Vg locations */
+  if(config->GetVGModel()!=ENUM_VG_MODEL::NONE) VG_Locations.resize(nPoint)=su2double(0.0);
+  // End added by max
 }
 
 void CFlowVariable::SetSolution_New() {
diff --git a/config_template.cfg b/config_template.cfg
index 74832ef98db..cbe14db8331 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -995,7 +995,22 @@ STREAMWISE_PERIODIC_TEMPERATURE= NO
 % Upon convergence, the area averaged inlet temperature will be INC_TEMPERATURE_INIT.
 % Defaults to 0.0.
 STREAMWISE_PERIODIC_OUTLET_HEAT= 0.0
-
+%
+% -------------------- BAY Source term vortex generator model -----------------%
+%
+% Chose vortex generator model (default=NONE, BAY, JBAY)
+VG_MODEL = NONE
+%
+% Specify BAY model constant. Default value is 10.
+VG_CONST = 10
+%
+% Specify vortex generator definition file. 
+% It should include the following colums to define the vortex generator:
+% L(lenght) h1(height leading edge) h2(height trailing edge) beta(VG angle in degrees) 
+% p1x p2y p1z(leaidng edge x,y,z coordinates) unx uny unz(surfce normal vector x,y,z components) 
+% ux uy uz(surface tangential vector x,y,z components)
+VG_CONFIG = vg.cfg
+%
 % -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
 %
 % Euler wall boundary marker(s) (NONE = no marker)