diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index d1f59abbd79..e3137e8869f 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1113,6 +1113,7 @@ enum class SA_OPTIONS {
   ROT,      /*!< \brief Rotation correction. */
   BC,       /*!< \brief Bas-Cakmakcioclu transition. */
   EXP,      /*!< \brief Allow experimental combinations of options (according to TMR). */
+  PGOMGA,   /*!< \brief Use PGOMGA term. */
 };
 static const MapType<std::string, SA_OPTIONS> SA_Options_Map = {
   MakePair("NONE", SA_OPTIONS::NONE)
@@ -1124,6 +1125,7 @@ static const MapType<std::string, SA_OPTIONS> SA_Options_Map = {
   MakePair("ROTATION", SA_OPTIONS::ROT)
   MakePair("BCM", SA_OPTIONS::BC)
   MakePair("EXPERIMENTAL", SA_OPTIONS::EXP)
+  MakePair("PGOMGA", SA_OPTIONS::PGOMGA)
 };
 
 /*!
@@ -1136,6 +1138,7 @@ struct SA_ParsedOptions {
   bool comp = false;                      /*!< \brief Use compressibility correction. */
   bool rot = false;                       /*!< \brief Use rotation correction. */
   bool bc = false;                        /*!< \brief BC transition. */
+  bool pgomga = false;                    /*!< \brief Use PGOMGA term. */
 };
 
 /*!
@@ -1173,6 +1176,7 @@ inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned sh
   SAParsedOptions.comp = IsPresent(SA_OPTIONS::COMP);
   SAParsedOptions.rot = IsPresent(SA_OPTIONS::ROT);
   SAParsedOptions.bc = IsPresent(SA_OPTIONS::BC);
+  SAParsedOptions.pgomga = IsPresent(SA_OPTIONS::PGOMGA);
 
   /*--- Validate user settings when not in experimental mode. ---*/
   if (!IsPresent(SA_OPTIONS::EXP)) {
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index 450112df8f7..773e00d5096 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -76,6 +76,7 @@ class CNumerics {
   su2double
   Eddy_Viscosity_i,  /*!< \brief Eddy viscosity at point i. */
   Eddy_Viscosity_j;  /*!< \brief Eddy viscosity at point j. */
+
   su2double
   turb_ke_i,  /*!< \brief Turbulent kinetic energy at point i. */
   turb_ke_j;  /*!< \brief Turbulent kinetic energy at point j. */
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index ce75b38f115..406306b06bd 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -50,9 +50,12 @@ struct CSAVariables {
   const su2double cr1 = 0.5;
   const su2double CRot = 1.0;
   const su2double c2 = 0.7, c3 = 0.9;
+  const su2double cpw1=0.94,cpw2=2.66,cpw3=4.8;//pgomega by Xiao He 2022
+  su2double beta_PGO = 0.0;//pgomega by Xiao He 2022
 
   /*--- List of auxiliary functions ---*/
   su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2;
+  
 
   /*--- List of helpers ---*/
   su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad;
@@ -129,7 +132,6 @@ class CSourceBase_TurbSA : public CNumerics {
     Jacobian_i = &Jacobian_Buffer;
   }
 
-
   /*!
    * \brief Residual for source term integration.
    * \param[in] config - Definition of the particular problem.
@@ -142,6 +144,7 @@ class CSourceBase_TurbSA : public CNumerics {
     AD::StartPreacc();
     AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, ScalarVar_i[0], Volume, dist_i, roughness_i);
     AD::SetPreaccIn(Vorticity_i, 3);
+    AD::SetPreaccIn(Coord_i, 3);
     AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim);
     AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim);
 
@@ -152,7 +155,43 @@ class CSourceBase_TurbSA : public CNumerics {
     Jacobian_i[0] = 0.0;
 
     /*--- Evaluate Omega with a rotational correction term. ---*/
-
+    // pgomega by Xiao He 2022
+      if (options.pgomga) {
+        
+        /*!\brief  PG-OMEGA*/
+        su2double dpds = 0, dpomg = 0, omg = 0, velocityMag = 0;
+        su2double Velocity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double Vorticity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double RotationalVelocity[3] = {0.0, 0.0, 0.0};
+        su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};
+        for (int iDim = 0; iDim < nDim; iDim++) {
+          RotationalVelocity[iDim] = config->GetRotation_Rate(iDim) / config->GetOmega_Ref();
+        }
+        
+        GeometryToolbox::CrossProduct( RotationalVelocity,Coord_i, RotationVelocityCrossR);
+        
+        for (int iDim = 0; iDim < nDim; iDim++) {
+          Velocity_Rel[iDim] = V_i[idx.Velocity() + iDim] - RotationVelocityCrossR[iDim];
+          Vorticity_Rel[iDim] = Vorticity_i[iDim] - RotationalVelocity[iDim];
+         
+        }
+
+        velocityMag = max(GeometryToolbox::Norm(3, Velocity_Rel), 0.001);
+        omg = max(GeometryToolbox::Norm(3, Vorticity_Rel), 0.001);
+        
+        for (int iDim = 0; iDim < nDim; iDim++) dpds += PrimVar_Grad_i[idx.Pressure()][iDim] * Velocity_Rel[iDim];
+
+        if (dpds >= 0.0) {
+          su2double RefReynold = 1.0e6;
+          dpomg = GeometryToolbox::DotProduct(3, PrimVar_Grad_i[idx.Pressure()], Vorticity_Rel) * laminar_viscosity *
+                  RefReynold / (omg * pow(velocityMag, 3) * pow(density, 2));  // reynold should be modified
+          //cout<<"dpomg="<<dpomg<<endl;
+          dpomg=abs(dpomg);
+          var.beta_PGO = var.cpw1 * tanh(var.cpw2 * pow(dpomg, var.cpw3));
+        }
+
+        //var.Omega = omg * (1.0 + var.beta_PGO);
+      }
     Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var);
 
     /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/
@@ -231,12 +270,11 @@ class CSourceBase_TurbSA : public CNumerics {
         var.intermittency = intermittency_eff_i;
         var.interDestrFactor = 1;
 
-      } else if (transition_LM){
-
+      } else if (transition_LM) {
         var.intermittency = intermittency_eff_i;
-        //var.intermittency = 1.0;
-        // Is wrong the reference from NASA?
-        // Original max(min(gamma, 0.5), 1.0) always gives 1 as result.
+        // var.intermittency = 1.0;
+        //  Is wrong the reference from NASA?
+        //  Original max(min(gamma, 0.5), 1.0) always gives 1 as result.
         var.interDestrFactor = min(max(intermittency_i, 0.5), 1.0);
 
       } else {
@@ -244,6 +282,8 @@ class CSourceBase_TurbSA : public CNumerics {
         var.intermittency = 1.0;
         var.interDestrFactor = 1.0;
       }
+      
+      
 
       /*--- Compute production, destruction and cross production and jacobian ---*/
       su2double Production = 0.0, Destruction = 0.0, CrossProduction = 0.0;
@@ -283,7 +323,7 @@ struct Omega {
 struct Bsl {
   template <class MatrixType>
   static void get(const su2double* vorticity, unsigned short, const MatrixType&, CSAVariables& var) {
-    var.Omega = GeometryToolbox::Norm(3, vorticity);
+    var.Omega = GeometryToolbox::Norm(3, vorticity)*(1.0 + var.beta_PGO);
   }
 };
 
diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp
index f7f7c4ea512..7d82cd4afab 100644
--- a/SU2_CFD/include/variables/CTurbSAVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp
@@ -42,6 +42,7 @@ class CTurbSAVariable final : public CTurbVariable {
   VectorType DES_LengthScale;
   VectorType Vortex_Tilting;
 
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -87,4 +88,5 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
 
+
 };
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 5ffee0d23f5..b9c47924383 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1413,6 +1413,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) {
 
   if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) {
     AddVolumeOutput("EDDY_VISCOSITY", "Eddy_Viscosity", "PRIMITIVE", "Turbulent eddy viscosity");
+
   }
 
 }
@@ -1544,6 +1545,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     SetVolumeOutputValue("EDDY_VISCOSITY", iPoint, Node_Flow->GetEddyViscosity(iPoint));
     SetVolumeOutputValue("TURB_DELTA_TIME", iPoint, Node_Turb->GetDelta_Time(iPoint));
     SetVolumeOutputValue("TURB_CFL", iPoint, Node_Turb->GetLocalCFL(iPoint));
+
   }
 
   if (config->GetSAParsedOptions().bc) {
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 40e19bd024c..856555a67cf 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -344,7 +344,9 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
 
     numerics->SetScalarVar(nodes->GetSolution(iPoint), nullptr);
     numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nullptr);
-
+    /*--- Set Coord ---*/
+    if(config->GetSAParsedOptions().pgomga)
+    numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr);
     /*--- Set volume ---*/
 
     numerics->SetVolume(geometry->nodes->GetVolume(iPoint));