diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index d53da0506cbf..95ed111e94a5 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -102,6 +102,7 @@ const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */
 
 const su2double ONE2 = 0.5;         /*!< \brief One divided by two. */
+const su2double ONE3 = 1.0 / 3.0;   /*!< \brief One divided by three. */
 const su2double TWO3 = 2.0 / 3.0;   /*!< \brief Two divided by three. */
 const su2double FOUR3 = 4.0 / 3.0;  /*!< \brief Four divided by three. */
 
@@ -2213,7 +2214,7 @@ enum MPI_QUANTITIES {
   SOLUTION_FEA_OLD     = 26,  /*!< \brief FEA solution old communication. */
   MESH_DISPLACEMENTS   = 27,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      = 28,  /*!< \brief Solution at time n. */
-  SOLUTION_TIME_N1     = 29,   /*!< \brief Solution at time n-1. */
+  SOLUTION_TIME_N1     = 29,  /*!< \brief Solution at time n-1. */
   PRIMITIVE            = 30   /*!< \brief Primitive solution communication. */
 };
 
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index 368fb4cd384b..81157b088ed2 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -181,8 +181,8 @@ class CNumerics {
   **TurbPsi_Grad_i,  /*!< \brief Gradient of adjoint turbulent variables at point i. */
   **TurbPsi_Grad_j;  /*!< \brief Gradient of adjoint turbulent variables at point j. */
   su2double
-  *AuxVar_Grad_i,    /*!< \brief Gradient of an auxiliary variable at point i. */
-  *AuxVar_Grad_j;    /*!< \brief Gradient of an auxiliary variable at point i. */
+  **AuxVar_Grad_i,    /*!< \brief Gradient of an auxiliary variable at point i. */
+  **AuxVar_Grad_j;    /*!< \brief Gradient of an auxiliary variable at point i. */
   const su2double *RadVar_Source;  /*!< \brief Source term from the radiative heat transfer equation. */
   su2double
   *Coord_i,      /*!< \brief Cartesians coordinates of point i. */
@@ -562,12 +562,13 @@ class CNumerics {
 
   /*!
    * \brief Set the gradient of the auxiliary variables.
-   * \param[in] val_auxvargrad_i - Gradient of the auxiliary variable at point i.
-   * \param[in] val_auxvargrad_j - Gradient of the auxiliary variable at point j.
+   * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i.
+   * \param[in] val_auxvar_grad_j - Gradient of the auxiliary variable at point j.
    */
-  inline void SetAuxVarGrad(su2double *val_auxvargrad_i, su2double *val_auxvargrad_j) {
-    AuxVar_Grad_i = val_auxvargrad_i;
-    AuxVar_Grad_j = val_auxvargrad_j;
+  inline void SetAuxVarGrad(su2double **val_auxvar_grad_i,
+                            su2double **val_auxvar_grad_j) {
+    AuxVar_Grad_i = val_auxvar_grad_i;
+    AuxVar_Grad_j = val_auxvar_grad_j;
   }
 
   /*!
diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp
index f9edd3bd7f0e..ea9df90c9119 100644
--- a/SU2_CFD/include/numerics/flow/flow_sources.hpp
+++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp
@@ -63,7 +63,16 @@ class CSourceBase_Flow : public CNumerics {
  * \ingroup SourceDiscr
  * \author F. Palacios
  */
-class CSourceAxisymmetric_Flow final : public CSourceBase_Flow {
+class CSourceAxisymmetric_Flow : public CSourceBase_Flow {
+protected:
+    bool implicit, viscous, rans;
+    su2double yinv{0.0};
+    
+  /*!
+  * \brief Diffusion residual of the axisymmetric source term.
+  */
+  void ResidualDiffusion();
+    
 public:
   /*!
    * \brief Constructor of the class.
@@ -74,12 +83,31 @@ class CSourceAxisymmetric_Flow final : public CSourceBase_Flow {
   CSourceAxisymmetric_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
 
   /*!
-   * \brief Residual of the rotational frame source term.
+   * \brief Residual of the axisymmetric source term.
    * \param[in] config - Definition of the particular problem.
    * \return Lightweight const-view of residual and Jacobian.
    */
   ResidualType<> ComputeResidual(const CConfig* config) override;
+  
+};
 
+/*!
+ * \class CSourceGeneralAxisymmetric_Flow
+ * \brief Class for source term for solving axisymmetric problems for a general (non ideal) fluid.
+ * \ingroup SourceDiscr
+ * \author F. Dittmann
+ */
+class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow {
+public:
+  
+  using CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow;
+  /*!
+   * \brief Residual of the general axisymmetric source term.
+   * \param[in] config - Definition of the particular problem.
+   * \return Lightweight const-view of residual and Jacobian.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override;
+  
 };
 
 /*!
diff --git a/SU2_CFD/include/solvers/CAdjEulerSolver.hpp b/SU2_CFD/include/solvers/CAdjEulerSolver.hpp
index d7884f7ee1d4..2d17d2f2e914 100644
--- a/SU2_CFD/include/solvers/CAdjEulerSolver.hpp
+++ b/SU2_CFD/include/solvers/CAdjEulerSolver.hpp
@@ -78,6 +78,13 @@ class CAdjEulerSolver : public CSolver {
    */
   inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }
 
+  /*!
+   * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config);
+
 public:
 
   /*!
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index e5b92958e152..74df30d7107a 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -585,13 +585,6 @@ class CSolver {
    */
   void SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config);
 
-  /*!
-   * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config);
-
   /*!
    * \brief Add External to Solution vector.
    */
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index c04937ec29ce..8d22ca3d9657 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -77,8 +77,8 @@ class CVariable {
   MatrixType Solution_Max;   /*!< \brief Max solution for limiter computation. */
   MatrixType Solution_Min;   /*!< \brief Min solution for limiter computation. */
 
-  VectorType AuxVar;       /*!< \brief Auxiliar variable for gradient computation. */
-  MatrixType Grad_AuxVar;  /*!< \brief Gradient of the auxiliar variable. */
+  MatrixType AuxVar;             /*!< \brief Auxiliar variable for gradient computation. */
+  CVectorOfMatrix Grad_AuxVar;   /*!< \brief Gradient of the auxilliary variables  of the problem. */
 
   VectorType Max_Lambda_Inv;   /*!< \brief Maximun inviscid eingenvalue. */
   VectorType Max_Lambda_Visc;  /*!< \brief Maximun viscous eingenvalue. */
@@ -105,6 +105,7 @@ class CVariable {
   unsigned long nPrimVarGrad = 0;    /*!< \brief Number of primitives for which a gradient is computed. */
   unsigned long nSecondaryVar = 0;     /*!< \brief Number of secondary variables. */
   unsigned long nSecondaryVarGrad = 0;   /*!< \brief Number of secondaries for which a gradient is computed. */
+  unsigned long nAuxVar = 0; /*!< \brief Number of auxiliary variables. */
 
   /*--- Only allow default construction by derived classes. ---*/
   CVariable() = default;
@@ -137,6 +138,11 @@ class CVariable {
    */
   virtual ~CVariable() = default;
 
+  /*!
+   * \brief Get the number of auxiliary variables.
+   */
+  inline unsigned long GetnAuxVar() const { return nAuxVar; }
+
   /*!
    * \brief Set the value of the solution, all variables.
    * \param[in] iPoint - Point index.
@@ -564,65 +570,62 @@ class CVariable {
   inline su2double GetLocalCFL(unsigned long iPoint) const { return LocalCFL(iPoint); }
 
   /*!
-   * \brief Set auxiliar variables, we are looking for the gradient of that variable.
-   * \param[in] iPoint - Point index.
-   * \param[in] val_auxvar - Value of the auxiliar variable.
+   * \brief Get the entire Aux matrix of the problem.
+   * \return Reference to the aux var  matrix.
    */
-  inline void SetAuxVar(unsigned long iPoint, su2double val_auxvar) { AuxVar(iPoint) = val_auxvar; }
+  inline const MatrixType& GetAuxVar(void) const { return AuxVar; }
 
   /*!
-   * \brief Get the value of the auxiliary variable.
-   * \param[in] iPoint - Point index.
-   * \return Value of the auxiliary variable.
+   * \brief Get the Aux var value at Point i, variable j.
    */
-  inline su2double GetAuxVar(unsigned long iPoint) const { return AuxVar(iPoint); }
+  inline su2double GetAuxVar(unsigned long iPoint, unsigned long iVar = 0) const { return AuxVar(iPoint,iVar); }
 
   /*!
-   * \brief Get the auxiliary variable.
-   * \return 2D view of the auxiliary variable.
+   * \brief Set auxiliary variables.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Varriable indexs
+   * \param[in] val_auxvar - Value of the auxiliar variable.
    */
-  inline C2DDummyLastView<const VectorType> GetAuxVar(void) const {
-    return C2DDummyLastView<const VectorType>(AuxVar);
+  inline void SetAuxVar(unsigned long iPoint, unsigned long iVar, const su2double auxvar) {
+    AuxVar(iPoint,iVar) = auxvar;
   }
 
   /*!
-   * \brief Set the value of the auxiliary variable gradient.
+   * \brief Set value of auxillary gradients.
    * \param[in] iPoint - Point index.
+   * \param[in] iVar - Index of the variable.
    * \param[in] iDim - Index of the dimension.
-   * \param[in] val_gradient - Value of the gradient for the index <i>iDim</i>.
+   * \param[in] value - Value of the gradient.
    */
-  inline void SetAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_gradient) { Grad_AuxVar(iPoint,iDim) = val_gradient; }
+  inline void SetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) {
+    Grad_AuxVar(iPoint,iVar,iDim) = value;
+  }
 
   /*!
-   * \brief Add a value to the auxiliary variable gradient.
-   * \param[in] iPoint - Point index.
-   * \param[in] iDim - Index of the dimension.
-   * \param[in] val_value - Value of the gradient to be added for the index <i>iDim</i>.
+   * \brief Get the gradient of the auxilary variables.
+   * \return Reference to gradient.
    */
-  inline void AddAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_value) { Grad_AuxVar(iPoint,iDim) += val_value;}
+  inline CVectorOfMatrix& GetAuxVarGradient(void) { return Grad_AuxVar; }
 
   /*!
-   * \brief Get the gradient of the auxiliary variable.
+   * \brief Get the value of the auxilliary gradient.
    * \param[in] iPoint - Point index.
-   * \return Value of the gradient of the auxiliary variable.
-   */
-  inline su2double *GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; }
-
-  /*!
-   * \brief Get the gradient of the auxiliary variable.
-   * \return 3D view of the gradient of the auxiliary variable.
+   * \param[in] iVar - Index of the variable.
+   * \param[in] iDim - Index of the dimension.
+   * \return Value of the solution gradient.
    */
-  inline C3DDummyMiddleView<MatrixType> GetAuxVarGradient() {
-    return C3DDummyMiddleView<MatrixType>(Grad_AuxVar);
+  inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const {
+    return Grad_AuxVar(iPoint,iVar,iDim);
   }
 
   /*!
-   * \brief Get the gradient of the auxiliary variable.
+   * \brief Get the value of the auxilliary gradient.
    * \param[in] iPoint - Point index.
-   * \param[in] iDim - Index of the dimension.
-   * \return Value of the gradient of the auxiliary variable for the dimension <i>iDim</i>.
+   * \return Value of the solution gradient.
    */
-  inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iDim) const { return Grad_AuxVar(iPoint,iDim); }
+  inline su2double** GetAuxVarGradient(unsigned long iPoint) {
+    return Grad_AuxVar[iPoint];
+  }
 
   /*!
    * \brief Add a value to the truncation error.
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 7443285b235e..55fdf052db5b 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -1864,8 +1864,10 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol
       else if (config->GetAxisymmetric() == YES) {
         if (incompressible)
           numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config);
-        else
+        else if (ideal_gas)
           numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config);
+        else
+          numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGeneralAxisymmetric_Flow(nDim, nVar_Flow, config);
       }
       else if (config->GetGravityForce() == YES) {
         numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGravity(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 8d21f94676b0..b3a48b19b2ce 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -1,4 +1,4 @@
-/*!
+/*!
  * \file flow_sources.cpp
  * \brief Implementation of numerics classes for integration
  *        of source terms in fluid flow problems.
@@ -50,16 +50,18 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi
 
   Gamma = config->GetGamma();
   Gamma_Minus_One = Gamma - 1.0;
+  
+  implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
+  viscous = config->GetViscous();
+  rans = (config->GetKind_Turb_Model() != NONE);
 
 }
 
 CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfig* config) {
 
-  su2double yinv, Pressure_i, Enthalpy_i, Velocity_i, sq_vel;
+  su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel;
   unsigned short iDim, iVar, jVar;
 
-  bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-
   if (Coord_i[1] > EPS) {
 
     yinv = 1.0/Coord_i[1];
@@ -77,6 +79,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
     residual[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0];
     residual[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]);
     residual[3] = yinv*Volume*Enthalpy_i*U_i[2];
+    
+    /*--- Inviscid component of the source term. ---*/
 
     if (implicit) {
       jacobian[0][0] = 0.0;
@@ -104,7 +108,117 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
           jacobian[iVar][jVar] *= yinv*Volume;
 
     }
+    
+    /*--- Add the viscous terms if necessary. ---*/
+
+    if (viscous) ResidualDiffusion();
+
+  }
+
+  else {
+
+    for (iVar=0; iVar < nVar; iVar++)
+      residual[iVar] = 0.0;
+
+    if (implicit) {
+      for (iVar=0; iVar < nVar; iVar++) {
+        for (jVar=0; jVar < nVar; jVar++)
+          jacobian[iVar][jVar] = 0.0;
+      }
+    }
+
+  }
+
+  return ResidualType<>(residual, jacobian, nullptr);
+}
+
+void CSourceAxisymmetric_Flow::ResidualDiffusion(){
+  
+  if (!rans){ turb_ke_i = 0.0; }
+  
+  su2double laminar_viscosity_i    = V_i[nDim+5];
+  su2double eddy_viscosity_i       = V_i[nDim+6];
+  su2double thermal_conductivity_i = V_i[nDim+7];
+  su2double heat_capacity_cp_i     = V_i[nDim+8];
+  
+  su2double total_viscosity_i = laminar_viscosity_i + eddy_viscosity_i;
+  su2double total_conductivity_i = thermal_conductivity_i + heat_capacity_cp_i*eddy_viscosity_i/Prandtl_Turb;
+  
+  su2double u = U_i[1]/U_i[0];
+  su2double v = U_i[2]/U_i[0];
+  
+  residual[0] -= 0.0;
+  residual[1] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][0]) 
+                                                                -TWO3*AuxVar_Grad_i[0][0]);
+  residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv)
+                                                                -TWO3*AuxVar_Grad_i[0][1]);
+  residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1])
+                                                 +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0]
+                                                 -v*yinv+U_i[0]*turb_ke_i))
+                                                 -total_conductivity_i*PrimVar_Grad_i[0][1])
+                                                 -TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1]));
+}
+  
+
+CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const CConfig* config) {
+  unsigned short iVar, jVar;
+
+  if (Coord_i[1] > EPS) {
+
+    yinv = 1.0/Coord_i[1];
+    
+    su2double Density_i = U_i[0];
+    su2double Velocity1_i = U_i[1]/U_i[0];
+    su2double Velocity2_i = U_i[2]/U_i[0];
+    su2double Energy_i = U_i[3]/U_i[0];
+    
+    su2double Pressure_i = V_j[3];
+    su2double Enthalpy_i = Energy_i + Pressure_i/Density_i;
+    
+    /*--- Inviscid component of the source term. ---*/
+
+    residual[0] = yinv*Volume*U_i[2];
+    residual[1] = yinv*Volume*U_i[1]*Velocity2_i;
+    residual[2] = yinv*Volume*U_i[2]*Velocity2_i;
+    residual[3] = yinv*Volume*U_i[2]*Enthalpy_i;
+
+    if (implicit) {
+      
+      su2double dPdrho_e_i = S_i[0];
+      su2double dPde_rho_i = S_i[1];
+    
+      jacobian[0][0] = 0.0;
+      jacobian[0][1] = 0.0;
+      jacobian[0][2] = 1.0;
+      jacobian[0][3] = 0.0;
+
+      jacobian[1][0] = -Velocity1_i*Velocity2_i;
+      jacobian[1][1] = Velocity2_i;
+      jacobian[1][2] = Velocity1_i;
+      jacobian[1][3] = 0.0;
+
+      jacobian[2][0] = -Velocity2_i*Velocity2_i;
+      jacobian[2][1] = 0.0;
+      jacobian[2][2] = 2*Velocity2_i;
+      jacobian[2][3] = 0.0;
+
+      jacobian[3][0] = Velocity2_i*(dPdrho_e_i + dPde_rho_i/Density_i*(Velocity1_i*Velocity1_i
+                                                                        + Velocity2_i*Velocity2_i
+                                                                        - Energy_i) - Enthalpy_i);
+      jacobian[3][1] = -Velocity1_i*Velocity2_i/Density_i *dPde_rho_i;
+      jacobian[3][2] = Enthalpy_i - Velocity2_i*Velocity2_i/Density_i *dPde_rho_i;
+      jacobian[3][3] = Velocity2_i + Velocity2_i/Density_i *dPde_rho_i;
+
+      for (iVar=0; iVar < nVar; iVar++)
+        for (jVar=0; jVar < nVar; jVar++)
+          jacobian[iVar][jVar] *= yinv*Volume;
+
+    }
+
+    /*--- Add the viscous terms if necessary. ---*/
 
+    if (viscous) ResidualDiffusion();
+    
   }
 
   else {
@@ -207,10 +321,10 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo
       /*--- Viscous terms. ---*/
 
       residual[0] -= 0.0;
-      residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0]);
+      residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0][0]);
       residual[2] -= Volume*(yinv*2.0*total_viscosity*PrimVar_Grad_i[2][1] -
-                                 yinv*yinv*2.0*total_viscosity*Velocity_i[1] -
-                                 TWO3*AuxVar_Grad_i[1]);
+                             yinv* yinv*2.0*total_viscosity*Velocity_i[1] -
+                             TWO3*AuxVar_Grad_i[0][1]);
       residual[3] -= Volume*yinv*Thermal_Conductivity_i*PrimVar_Grad_i[nDim+1][1];
 
     }
diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
index 828ad71e3069..3da1fca0a123 100644
--- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
@@ -2385,7 +2385,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver
           conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1];
           for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1];
 
-          nodes->SetAuxVar(iPoint,conspsi);
+          nodes->SetAuxVar(iPoint,0,conspsi);
 
           /*--- Also load the auxiliary variable for first neighbors ---*/
 
@@ -2396,7 +2396,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver
             Enthalpy = solver_container[FLOW_SOL]->GetNodes()->GetEnthalpy(Neigh);
             conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1];
             for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1];
-            nodes->SetAuxVar(Neigh,conspsi);
+            nodes->SetAuxVar(Neigh,0,conspsi);
           }
         }
       }
@@ -2423,7 +2423,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver
           Area = sqrt(Area);
 
           PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint);
-          ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint);
+          ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint)[0];
           ConsPsi = nodes->GetAuxVar(iPoint);
 
           d_press = 0.0; grad_v = 0.0; v_gradconspsi = 0.0;
@@ -4884,3 +4884,115 @@ void CAdjEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf
   Restart_Vars = nullptr; Restart_Data = nullptr;
 
 }
+
+void CAdjEulerSolver::SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config) {
+
+  unsigned short iDim, jDim, iNeigh, iMarker;
+  unsigned short nDim = geometry->GetnDim();
+  unsigned long iPoint, jPoint, iVertex;
+  su2double *Coord_i, *Coord_j, AuxVar_i, AuxVar_j;
+  su2double **Smatrix, *Cvector;
+
+  Smatrix = new su2double* [nDim];
+  Cvector = new su2double [nDim];
+  for (iDim = 0; iDim < nDim; iDim++)
+    Smatrix[iDim] = new su2double [nDim];
+
+
+  /*--- Loop over boundary markers to select those for Euler or NS walls ---*/
+  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+
+    if (config->GetSolid_Wall(iMarker)) {
+
+      /*--- Loop over points on the surface (Least-Squares approximation) ---*/
+      for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+        iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+        if (geometry->nodes->GetDomain(iPoint)) {
+          Coord_i = geometry->nodes->GetCoord(iPoint);
+          AuxVar_i = nodes->GetAuxVar(iPoint);
+
+          /*--- Inizialization of variables ---*/
+          for (iDim = 0; iDim < nDim; iDim++)
+            Cvector[iDim] = 0.0;
+          su2double r11 = 0.0, r12 = 0.0, r13 = 0.0, r22 = 0.0, r23 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0;
+
+          for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) {
+            jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+            Coord_j = geometry->nodes->GetCoord(jPoint);
+            AuxVar_j = nodes->GetAuxVar(jPoint);
+
+            su2double weight = 0;
+            for (iDim = 0; iDim < nDim; iDim++)
+              weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]);
+
+            /*--- Sumations for entries of upper triangular matrix R ---*/
+            r11 += (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight;
+            r12 += (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight;
+            r22 += (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight;
+            if (nDim == 3) {
+              r13 += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight;
+              r23_a += (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight;
+              r23_b += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight;
+              r33 += (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight;
+            }
+
+            /*--- Entries of c:= transpose(A)*b ---*/
+            for (iDim = 0; iDim < nDim; iDim++)
+              Cvector[iDim] += (Coord_j[iDim]-Coord_i[iDim])*(AuxVar_j-AuxVar_i)/weight;
+          }
+
+          /*--- Entries of upper triangular matrix R ---*/
+          r11 = sqrt(r11);
+          r12 = r12/r11;
+          r22 = sqrt(r22-r12*r12);
+          if (nDim == 3) {
+            r13 = r13/r11;
+            r23 = r23_a/r22 - r23_b*r12/(r11*r22);
+            r33 = sqrt(r33-r23*r23-r13*r13);
+          }
+          /*--- S matrix := inv(R)*traspose(inv(R)) ---*/
+          if (nDim == 2) {
+            su2double detR2 = (r11*r22)*(r11*r22);
+            Smatrix[0][0] = (r12*r12+r22*r22)/detR2;
+            Smatrix[0][1] = -r11*r12/detR2;
+            Smatrix[1][0] = Smatrix[0][1];
+            Smatrix[1][1] = r11*r11/detR2;
+          }
+          else {
+            su2double detR2 = (r11*r22*r33)*(r11*r22*r33);
+            su2double z11, z12, z13, z22, z23, z33; // aux vars
+            z11 = r22*r33;
+            z12 = -r12*r33;
+            z13 = r12*r23-r13*r22;
+            z22 = r11*r33;
+            z23 = -r11*r23;
+            z33 = r11*r22;
+            Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2;
+            Smatrix[0][1] = (z12*z22+z13*z23)/detR2;
+            Smatrix[0][2] = (z13*z33)/detR2;
+            Smatrix[1][0] = Smatrix[0][1];
+            Smatrix[1][1] = (z22*z22+z23*z23)/detR2;
+            Smatrix[1][2] = (z23*z33)/detR2;
+            Smatrix[2][0] = Smatrix[0][2];
+            Smatrix[2][1] = Smatrix[1][2];
+            Smatrix[2][2] = (z33*z33)/detR2;
+          }
+          /*--- Computation of the gradient: S*c ---*/
+          su2double product;
+          for (iDim = 0; iDim < nDim; iDim++) {
+            product = 0.0;
+            for (jDim = 0; jDim < nDim; jDim++)
+              product += Smatrix[iDim][jDim]*Cvector[jDim];
+            nodes->SetAuxVarGradient(iPoint, 0, iDim, product);
+          }
+        }
+      } /*--- End of loop over surface points ---*/
+    }
+  }
+
+  /*--- Memory deallocation ---*/
+  for (iDim = 0; iDim < nDim; iDim++)
+    delete [] Smatrix[iDim];
+  delete [] Cvector;
+  delete [] Smatrix;
+}
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index e5e3c966957e..0dfa329c4ed1 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -3039,13 +3039,17 @@ void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned
 void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,
                                    CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
 
-  const bool implicit         = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const bool implicit         = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+  const bool viscous          = config->GetViscous();
   const bool rotating_frame   = config->GetRotating_Frame();
   const bool axisymmetric     = config->GetAxisymmetric();
   const bool gravity          = (config->GetGravityForce() == YES);
   const bool harmonic_balance = (config->GetTime_Marching() == HARMONIC_BALANCE);
   const bool windgust         = config->GetWind_Gust();
   const bool body_force       = config->GetBody_Force();
+  const bool ideal_gas        = (config->GetKind_FluidModel() == STANDARD_AIR) ||
+                                (config->GetKind_FluidModel() == IDEAL_GAS);
+  const bool rans             = (config->GetKind_Turb_Model() != NONE);
 
   /*--- Pick one numerics object per thread. ---*/
   CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS];
@@ -3107,6 +3111,34 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
   if (axisymmetric) {
 
+    /*--- For viscous problems, we need an additional gradient. ---*/
+    if (viscous) {
+
+      for (iPoint = 0; iPoint < nPoint; iPoint++) {
+
+        su2double yCoord          = geometry->nodes->GetCoord(iPoint, 1);
+        su2double yVelocity       = nodes->GetVelocity(iPoint,1);
+        su2double xVelocity       = nodes->GetVelocity(iPoint,0);
+        su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint);
+
+        if (yCoord > EPS){
+          su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord;
+          nodes->SetAuxVar(iPoint, 0, nu_v_on_y);
+          nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity);
+          nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity);
+        }
+      }
+
+      /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/
+      if (config->GetKind_Gradient_Method() == GREEN_GAUSS) {
+        SetAuxVar_Gradient_GG(geometry, config);
+      }
+      if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) {
+        SetAuxVar_Gradient_LS(geometry, config);
+      }
+
+    }
+
     /*--- loop over points ---*/
     SU2_OMP_FOR_DYN(omp_chunk_size)
     for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
@@ -3120,6 +3152,27 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
       /*--- Set y coordinate ---*/
       numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
 
+      /*--- Set primitive variables for viscous terms and/or generalised source ---*/
+      if (!ideal_gas || viscous) numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));
+
+      /*--- Set secondary variables for generalised source ---*/
+      if (!ideal_gas) numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
+
+      if (viscous) {
+
+        /*--- Set gradient of primitive variables ---*/
+        numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
+
+        /*--- Set gradient of auxillary variables ---*/
+        numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
+
+        /*--- Set turbulence kinetic energy ---*/
+        if (rans){
+          CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
+          numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0));
+        }
+      }
+
       /*--- Compute Source term Residual ---*/
       auto residual = numerics->ComputeResidual(config);
 
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index 8fce143bfd0d..cdc7f9ff6a1d 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -3475,7 +3475,7 @@ void CFEASolver::FilterElementDensities(CGeometry *geometry, const CConfig *conf
         sum += w * element_properties[iElem]->GetPhysicalDensity();
         vol += w;
       }
-      nodes->SetAuxVar(iPoint, sum/vol);
+      nodes->SetAuxVar(iPoint, 0, sum/vol);
     }
   }
 
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 3de75bc87501..991c10eaeb08 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -1572,7 +1572,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
 
         /*--- Set the auxilairy variable for this node. ---*/
 
-        nodes->SetAuxVar(iPoint, AuxVar);
+        nodes->SetAuxVar(iPoint, 0, AuxVar);
 
       }
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 6743086dca2f..ede7553dd01a 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -1647,7 +1647,7 @@ void CSolver::GetCommCountAndType(const CConfig* config,
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       break;
     case AUXVAR_GRADIENT:
-      COUNT_PER_POINT  = nDim;
+      COUNT_PER_POINT  = nDim*base_nodes->GetnAuxVar();
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       break;
     case MESH_DISPLACEMENTS:
@@ -1778,8 +1778,11 @@ void CSolver::InitiateComms(CGeometry *geometry,
               bufDSend[buf_offset+iVar] = base_nodes->GetLimiter_Primitive(iPoint, iVar);
             break;
           case AUXVAR_GRADIENT:
-            for (iDim = 0; iDim < nDim; iDim++)
-              bufDSend[buf_offset+iDim] = base_nodes->GetAuxVarGradient(iPoint, iDim);
+            for (iVar = 0; iVar < base_nodes->GetnAuxVar(); iVar++){
+              for (iDim = 0; iDim < nDim; iDim++){
+                bufDSend[buf_offset+iVar*nDim+iDim] = base_nodes->GetAuxVarGradient(iPoint, iVar, iDim);
+              }
+            }
             break;
           case SOLUTION_FEA:
             for (iVar = 0; iVar < nVar; iVar++) {
@@ -1953,8 +1956,11 @@ void CSolver::CompleteComms(CGeometry *geometry,
               base_nodes->SetLimiter_Primitive(iPoint, iVar, bufDRecv[buf_offset+iVar]);
             break;
           case AUXVAR_GRADIENT:
-            for (iDim = 0; iDim < nDim; iDim++)
-              base_nodes->SetAuxVarGradient(iPoint, iDim, bufDRecv[buf_offset+iDim]);
+            for (iVar = 0; iVar < base_nodes->GetnAuxVar(); iVar++){
+              for (iDim = 0; iDim < nDim; iDim++){
+                base_nodes->SetAuxVarGradient(iPoint, iVar, iDim, bufDRecv[buf_offset+iVar*nDim+iDim]);
+              }
+            }
             break;
           case SOLUTION_FEA:
             for (iVar = 0; iVar < nVar; iVar++) {
@@ -2500,7 +2506,7 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config)
   auto gradient = base_nodes->GetAuxVarGradient();
 
   computeGradientsGreenGauss(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry,
-                             *config, solution, 0, 1, gradient);
+                             *config, solution, 0, base_nodes->GetnAuxVar(), gradient);
 }
 
 void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) {
@@ -2511,7 +2517,7 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config)
   auto& rmatrix  = base_nodes->GetRmatrix();
 
   computeGradientsLeastSquares(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config,
-                               weighted, solution, 0, 1, gradient, rmatrix);
+                               weighted, solution, 0, base_nodes->GetnAuxVar(), gradient, rmatrix);
 }
 
 void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config, bool reconstruction) {
@@ -2592,118 +2598,6 @@ void CSolver::SetGridVel_Gradient(CGeometry *geometry, const CConfig *config) {
                                true, gridVel, 0, nDim, gridVelGrad, rmatrix);
 }
 
-void CSolver::SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config) {
-
-  unsigned short iDim, jDim, iNeigh, iMarker;
-  unsigned short nDim = geometry->GetnDim();
-  unsigned long iPoint, jPoint, iVertex;
-  su2double *Coord_i, *Coord_j, AuxVar_i, AuxVar_j;
-  su2double **Smatrix, *Cvector;
-
-  Smatrix = new su2double* [nDim];
-  Cvector = new su2double [nDim];
-  for (iDim = 0; iDim < nDim; iDim++)
-    Smatrix[iDim] = new su2double [nDim];
-
-
-  /*--- Loop over boundary markers to select those for Euler or NS walls ---*/
-  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-
-    if (config->GetSolid_Wall(iMarker)) {
-
-      /*--- Loop over points on the surface (Least-Squares approximation) ---*/
-      for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
-        iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
-        if (geometry->nodes->GetDomain(iPoint)) {
-          Coord_i = geometry->nodes->GetCoord(iPoint);
-          AuxVar_i = base_nodes->GetAuxVar(iPoint);
-
-          /*--- Inizialization of variables ---*/
-          for (iDim = 0; iDim < nDim; iDim++)
-            Cvector[iDim] = 0.0;
-          su2double r11 = 0.0, r12 = 0.0, r13 = 0.0, r22 = 0.0, r23 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0;
-
-          for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) {
-            jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
-            Coord_j = geometry->nodes->GetCoord(jPoint);
-            AuxVar_j = base_nodes->GetAuxVar(jPoint);
-
-            su2double weight = 0;
-            for (iDim = 0; iDim < nDim; iDim++)
-              weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]);
-
-            /*--- Sumations for entries of upper triangular matrix R ---*/
-            r11 += (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight;
-            r12 += (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight;
-            r22 += (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight;
-            if (nDim == 3) {
-              r13 += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight;
-              r23_a += (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight;
-              r23_b += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight;
-              r33 += (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight;
-            }
-
-            /*--- Entries of c:= transpose(A)*b ---*/
-            for (iDim = 0; iDim < nDim; iDim++)
-              Cvector[iDim] += (Coord_j[iDim]-Coord_i[iDim])*(AuxVar_j-AuxVar_i)/weight;
-          }
-
-          /*--- Entries of upper triangular matrix R ---*/
-          r11 = sqrt(r11);
-          r12 = r12/r11;
-          r22 = sqrt(r22-r12*r12);
-          if (nDim == 3) {
-            r13 = r13/r11;
-            r23 = r23_a/r22 - r23_b*r12/(r11*r22);
-            r33 = sqrt(r33-r23*r23-r13*r13);
-          }
-          /*--- S matrix := inv(R)*traspose(inv(R)) ---*/
-          if (nDim == 2) {
-            su2double detR2 = (r11*r22)*(r11*r22);
-            Smatrix[0][0] = (r12*r12+r22*r22)/detR2;
-            Smatrix[0][1] = -r11*r12/detR2;
-            Smatrix[1][0] = Smatrix[0][1];
-            Smatrix[1][1] = r11*r11/detR2;
-          }
-          else {
-            su2double detR2 = (r11*r22*r33)*(r11*r22*r33);
-            su2double z11, z12, z13, z22, z23, z33; // aux vars
-            z11 = r22*r33;
-            z12 = -r12*r33;
-            z13 = r12*r23-r13*r22;
-            z22 = r11*r33;
-            z23 = -r11*r23;
-            z33 = r11*r22;
-            Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2;
-            Smatrix[0][1] = (z12*z22+z13*z23)/detR2;
-            Smatrix[0][2] = (z13*z33)/detR2;
-            Smatrix[1][0] = Smatrix[0][1];
-            Smatrix[1][1] = (z22*z22+z23*z23)/detR2;
-            Smatrix[1][2] = (z23*z33)/detR2;
-            Smatrix[2][0] = Smatrix[0][2];
-            Smatrix[2][1] = Smatrix[1][2];
-            Smatrix[2][2] = (z33*z33)/detR2;
-          }
-          /*--- Computation of the gradient: S*c ---*/
-          su2double product;
-          for (iDim = 0; iDim < nDim; iDim++) {
-            product = 0.0;
-            for (jDim = 0; jDim < nDim; jDim++)
-              product += Smatrix[iDim][jDim]*Cvector[jDim];
-            base_nodes->SetAuxVarGradient(iPoint, iDim, product);
-          }
-        }
-      } /*--- End of loop over surface points ---*/
-    }
-  }
-
-  /*--- Memory deallocation ---*/
-  for (iDim = 0; iDim < nDim; iDim++)
-    delete [] Smatrix[iDim];
-  delete [] Cvector;
-  delete [] Smatrix;
-}
-
 void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) {
 
   auto kindLimiter = static_cast<ENUM_LIMITER>(config->GetKind_SlopeLimit());
diff --git a/SU2_CFD/src/variables/CAdjEulerVariable.cpp b/SU2_CFD/src/variables/CAdjEulerVariable.cpp
index 09458a37942a..4081cd0b75f1 100644
--- a/SU2_CFD/src/variables/CAdjEulerVariable.cpp
+++ b/SU2_CFD/src/variables/CAdjEulerVariable.cpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -52,7 +52,7 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2
   if (config->GetReconstructionGradientRequired()) {
     Gradient_Aux.resize(nPoint,nVar,nDim,0.0);
   }
-  
+
   if (config->GetLeastSquaresRequired()) {
     Rmatrix.resize(nPoint,nDim,nDim,0.0);
   }
@@ -84,8 +84,9 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2
   }
 
   /*--- Allocate auxiliar vector for sensitivity computation ---*/
-  AuxVar.resize(nPoint);
-  Grad_AuxVar.resize(nPoint,nDim);
+  nAuxVar = 1;
+  AuxVar.resize(nPoint, nAuxVar) = su2double(0.0);
+  Grad_AuxVar.resize(nPoint, nAuxVar, nDim);
 
   /*--- Allocate and initializate projection vector for wall boundary condition ---*/
   ForceProj_Vector.resize(nPoint,nDim) = su2double(0.0);
@@ -101,11 +102,11 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2
     Set_BGSSolution_k();
 
   Sensor.resize(nPoint);
-  
+
   /* Non-physical point (first-order) initialization. */
   Non_Physical.resize(nPoint) = false;
   Non_Physical_Counter.resize(nPoint) = 0;
-  
+
 }
 
 bool CAdjEulerVariable::SetPrimVar(unsigned long iPoint, su2double SharpEdge_Distance, bool check, CConfig *config) {
diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp
index f791136ebe5a..43ade694b9d6 100644
--- a/SU2_CFD/src/variables/CEulerVariable.cpp
+++ b/SU2_CFD/src/variables/CEulerVariable.cpp
@@ -121,6 +121,12 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2
     Gradient_Aux.resize(nPoint,nPrimVarGrad,nDim,0.0);
   }
 
+  if (config->GetAxisymmetric()){
+    nAuxVar = 3;
+    Grad_AuxVar.resize(nPoint,nAuxVar,nDim,0.0);
+    AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);
+  }
+
   if (config->GetLeastSquaresRequired()) {
     Rmatrix.resize(nPoint,nDim,nDim,0.0);
   }
diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp
index 7e6d0c3c7525..376db841bc69 100644
--- a/SU2_CFD/src/variables/CFEAVariable.cpp
+++ b/SU2_CFD/src/variables/CFEAVariable.cpp
@@ -88,7 +88,10 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig
 
   if (multizone) Set_BGSSolution_k();
 
-  if (config->GetTopology_Optimization()) AuxVar.resize(nPoint);
+  if (config->GetTopology_Optimization()) {
+    nAuxVar = 1;
+    AuxVar.resize(nPoint);
+  }
 }
 
 void CFEAVariable::SetSolution_Vel_time_n() { Solution_Vel_time_n = Solution_Vel; }
diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp
index 72b11117b91c..47fee15f56f7 100644
--- a/SU2_CFD/src/variables/CIncEulerVariable.cpp
+++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp
@@ -34,8 +34,6 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
 
   bool dual_time    = (config->GetTime_Marching() == DT_STEPPING_1ST) ||
                       (config->GetTime_Marching() == DT_STEPPING_2ND);
-  bool viscous      = config->GetViscous();
-  bool axisymmetric = config->GetAxisymmetric();
 
   /*--- Allocate and initialize the primitive variables and gradients ---*/
 
@@ -105,10 +103,6 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
     Rmatrix.resize(nPoint,nDim,nDim,0.0);
   }
 
-  /*--- If axisymmetric and viscous, we need an auxiliary gradient. ---*/
-
-  if (axisymmetric && viscous) Grad_AuxVar.resize(nPoint,nDim);
-
   if (config->GetMultizone_Problem())
     Set_BGSSolution_k();
 
diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp
index e08933790ecb..06a5fdbd67d8 100644
--- a/SU2_CFD/src/variables/CIncNSVariable.cpp
+++ b/SU2_CFD/src/variables/CIncNSVariable.cpp
@@ -31,10 +31,17 @@
 CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature,
                                unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) :
                                CIncEulerVariable(pressure, velocity, temperature, npoint, ndim, nvar, config) {
+
   Vorticity.resize(nPoint,3);
   StrainMag.resize(nPoint);
   DES_LengthScale.resize(nPoint) = su2double(0.0);
   Max_Lambda_Visc.resize(nPoint);
+
+  if (config->GetAxisymmetric()) {
+    nAuxVar = 1;
+    AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);
+    Grad_AuxVar.resize(nPoint,nAuxVar,nDim);
+  }
 }
 
 bool CIncNSVariable::SetVorticity_StrainMag() {