diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 8a51da7a7286..3ad89a515039 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -4662,8 +4662,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Streamwise Periodic Flow currently only implemented for incompressible flow.", CURRENT_FUNCTION);
     if (Kind_Solver == INC_EULER)
       SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
-    if (nMarker_PerBound != 2)
-      SU2_MPI::Error("Streamwise Periodic Flow currently only implemented for one Periodic Marker pair. Combining Streamwise and Spanwise periodicity not possible in the moment.", CURRENT_FUNCTION);
+    if (nMarker_PerBound == 0)
+      SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
     if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
       SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
     if (DiscreteAdjoint && Kind_Streamwise_Periodic == ENUM_STREAMWISE_PERIODIC::MASSFLOW)
diff --git a/SU2_CFD/include/interfaces/CInterface.hpp b/SU2_CFD/include/interfaces/CInterface.hpp
index 0bd9e4f70fe2..ea4b487cbfb4 100644
--- a/SU2_CFD/include/interfaces/CInterface.hpp
+++ b/SU2_CFD/include/interfaces/CInterface.hpp
@@ -70,6 +70,7 @@ class CInterface {
   unsigned short nSpanMaxAllZones = 0;
 
   unsigned short nVar = 0;
+  static constexpr size_t MAXNDIM = 3;  /*!< \brief Max number of space dimensions, used in some static arrays. */
 
 public:
   /*!
diff --git a/SU2_CFD/include/output/CAdjFlowOutput.hpp b/SU2_CFD/include/output/CAdjFlowCompOutput.hpp
similarity index 100%
rename from SU2_CFD/include/output/CAdjFlowOutput.hpp
rename to SU2_CFD/include/output/CAdjFlowCompOutput.hpp
diff --git a/SU2_CFD/include/output/CAdjHeatOutput.hpp b/SU2_CFD/include/output/CAdjHeatOutput.hpp
index 5b5586591ded..2102be9ad849 100644
--- a/SU2_CFD/include/output/CAdjHeatOutput.hpp
+++ b/SU2_CFD/include/output/CAdjHeatOutput.hpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_fea_discadj.hpp
+ * \file CAdjHeatOutput.hpp
  * \brief Headers of the adjoint heat output.
  * \author T. Albring
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/include/output/CFVMOutput.hpp b/SU2_CFD/include/output/CFVMOutput.hpp
new file mode 100644
index 000000000000..7e7da9e7a380
--- /dev/null
+++ b/SU2_CFD/include/output/CFVMOutput.hpp
@@ -0,0 +1,69 @@
+/*!
+ * \file CFVMOutput.hpp
+ * \brief  Headers of the Finite Volume Method output.
+ * \author T. Kattmann
+ * \version 7.1.1 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2021, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "COutput.hpp"
+
+class CFVMOutput : public COutput{
+protected:
+  /*!
+   * \brief Constructor of the class
+   */
+  CFVMOutput(CConfig *config, unsigned short nDim, bool femOutput);
+
+  /*!
+   * \brief Destructor of the class.
+   */
+  ~CFVMOutput(void) = default;
+
+  /*!
+   * \brief Add Coordinates to output.
+   */
+  void AddCoordinates();
+
+  /*!
+   * \brief Load the coordinates.
+   */
+  template<class T>
+  inline void LoadCoordinates(const T& Coord, const unsigned long iPoint) {
+    SetVolumeOutputValue("COORD-X", iPoint, Coord[0]);
+    SetVolumeOutputValue("COORD-Y", iPoint, Coord[1]);
+    if (nDim == 3)
+      SetVolumeOutputValue("COORD-Z", iPoint, Coord[2]);
+  }
+
+  /*!
+   * \brief Add common FVM outputs.
+   */
+  void AddCommonFVMOutputs(const CConfig* config);
+
+  /*!
+   * \brief Load common FVM outputs.
+   */
+  void LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint);
+};
diff --git a/SU2_CFD/include/output/CFlowIncOutput.hpp b/SU2_CFD/include/output/CFlowIncOutput.hpp
index af106a330932..9e1618c697ee 100644
--- a/SU2_CFD/include/output/CFlowIncOutput.hpp
+++ b/SU2_CFD/include/output/CFlowIncOutput.hpp
@@ -1,5 +1,5 @@
-/*!
- * \file CFlowIncCompOutput.hpp
+/*!
+ * \file CFlowIncOutput.hpp
  * \brief  Headers of the incompressible flow output.
  * \author T. Albring, R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp
index 8e15d607ae1d..b5665420442c 100644
--- a/SU2_CFD/include/output/CFlowOutput.hpp
+++ b/SU2_CFD/include/output/CFlowOutput.hpp
@@ -27,10 +27,10 @@
 
 #pragma once
 
-#include "COutput.hpp"
+#include "CFVMOutput.hpp"
 #include "../variables/CVariable.hpp"
 
-class CFlowOutput : public COutput{
+class CFlowOutput : public CFVMOutput{
 public:
   /*!
    * \brief Constructor of the class
@@ -142,14 +142,4 @@ class CFlowOutput : public COutput{
    */
   void LoadTimeAveragedData(unsigned long iPoint, CVariable *node_flow);
 
-  /*!
-   * \brief Add common FVM outputs.
-   */
-  void AddCommonFVMOutputs(const CConfig* config);
-
-  /*!
-   * \brief Load common FVM outputs.
-   */
-  void LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint);
-
 };
diff --git a/SU2_CFD/include/output/CHeatOutput.hpp b/SU2_CFD/include/output/CHeatOutput.hpp
index cb9b2b1b0f4e..06f37d1a3344 100644
--- a/SU2_CFD/include/output/CHeatOutput.hpp
+++ b/SU2_CFD/include/output/CHeatOutput.hpp
@@ -1,6 +1,6 @@
 /*!
  * \file CHeatOutput.hpp
- * \brief  Headers of the heat output.
+ * \brief Headers of the heat output.
  * \author R. Sanchez, T. Albring.
  * \version 7.1.1 "Blackbird"
  *
@@ -27,7 +27,7 @@
 
 #pragma once
 
-#include "COutput.hpp"
+#include "CFVMOutput.hpp"
 
 
 /*! \class CHeatOutput
@@ -35,7 +35,7 @@
  *  \author R. Sanchez, T. Albring.
  *  \date June 5, 2018.
  */
-class CHeatOutput final: public COutput {
+class CHeatOutput final: public CFVMOutput {
 public:
 
   /*!
diff --git a/SU2_CFD/include/output/CMeshOutput.hpp b/SU2_CFD/include/output/CMeshOutput.hpp
index 5c736090c4cc..34c01ebc4357 100644
--- a/SU2_CFD/include/output/CMeshOutput.hpp
+++ b/SU2_CFD/include/output/CMeshOutput.hpp
@@ -1,6 +1,6 @@
 /*!
  * \file CMeshOutput.hpp
- * \brief  Headers of the mesh output.
+ * \brief Headers of the mesh output.
  * \author R. Sanchez, T. Albring.
  * \version 7.1.1 "Blackbird"
  *
diff --git a/SU2_CFD/include/output/CMultizoneOutput.hpp b/SU2_CFD/include/output/CMultizoneOutput.hpp
index 1be5ea3e753c..28ff8d1b5557 100644
--- a/SU2_CFD/include/output/CMultizoneOutput.hpp
+++ b/SU2_CFD/include/output/CMultizoneOutput.hpp
@@ -1,5 +1,5 @@
 /*!
- * \file CDriverOutput.hpp
+ * \file CMultizoneOutput.hpp
  * \brief Headers of the main subroutines for screen and history output in multizone problems.
  * \author R. Sanchez, T. Albring
  * \version 7.1.1 "Blackbird"
@@ -44,7 +44,7 @@
 using namespace std;
 
 /*!
- * \class CDriverOutput
+ * \class CMultizoneOutput
  * \brief Class for writing the multizone output.
  * \author R. Sanchez, T. Albring.
  */
@@ -53,7 +53,7 @@ class CMultizoneOutput final: public COutput {
 protected:
   unsigned short nZone; //!< Number of zones
 
-  string bgs_res_name; //!< Block-Gauss seidel residual name
+  string bgs_res_name; //!< Block-Gauss Seidel residual name
   bool write_zone;     //!< Boolean indicating whether the individual zones write to screen
 
 public:
diff --git a/SU2_CFD/include/output/COutputFactory.hpp b/SU2_CFD/include/output/COutputFactory.hpp
index 3e659e979546..0aa2dd25e457 100644
--- a/SU2_CFD/include/output/COutputFactory.hpp
+++ b/SU2_CFD/include/output/COutputFactory.hpp
@@ -1,5 +1,5 @@
 /*!
- * \file COutput.hpp
+ * \file COutputFactory.hpp
  * \brief Headers of the output class.
  * \author T.Albring
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/include/output/COutputLegacy.hpp b/SU2_CFD/include/output/COutputLegacy.hpp
index fbf334a66cab..9f34d8e45292 100644
--- a/SU2_CFD/include/output/COutputLegacy.hpp
+++ b/SU2_CFD/include/output/COutputLegacy.hpp
@@ -1,7 +1,7 @@
 /*!
- * \file output.hpp
+ * \file COutputLegacy.hpp
  * \brief Headers of the main subroutines for generating the file outputs.
- *        The subroutines and functions are in the <i>output_structure.cpp</i> file.
+ *        The subroutines and functions are in the <i>output_structure_legacy.cpp</i> file.
  * \author F. Palacios, T. Economon, M. Colonno
  * \version 7.1.1 "Blackbird"
  *
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index 1898006dbac1..49fac6b79b7b 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -78,7 +78,7 @@ class CVariable {
   MatrixType Solution_Min;   /*!< \brief Min solution for limiter computation. */
 
   MatrixType AuxVar;             /*!< \brief Auxiliary variable for gradient computation. */
-  CVectorOfMatrix Grad_AuxVar;   /*!< \brief Gradient of the auxiliary variables  of the problem. */
+  CVectorOfMatrix Grad_AuxVar;   /*!< \brief Gradient of the auxiliary variables of the problem. */
 
   VectorType Max_Lambda_Inv;   /*!< \brief Maximun inviscid eingenvalue. */
   VectorType Max_Lambda_Visc;  /*!< \brief Maximun viscous eingenvalue. */
diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am
index 5f7975302eb8..6b7f7854ac9c 100644
--- a/SU2_CFD/obj/Makefile.am
+++ b/SU2_CFD/obj/Makefile.am
@@ -136,6 +136,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \
   ../src/output/output_physics.cpp \
   ../src/output/CMeshOutput.cpp \
   ../src/output/CElasticityOutput.cpp \
+  ../src/output/CFVMOutput.cpp \
   ../src/output/CFlowOutput.cpp \
   ../src/output/CFlowCompOutput.cpp \
   ../src/output/CFlowCompFEMOutput.cpp \
diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
index 885d155fb2b7..657785783e85 100644
--- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
@@ -345,18 +345,6 @@ void CDiscAdjSinglezoneDriver::SetObjFunction(){
 
     solver[FLOW_SOL]->SetTotal_ComboObj(0.0);
 
-//    if (config->GetnMarker_Analyze() != 0)
-//      output->SpecialOutput_AnalyzeSurface(solver[FLOW_SOL], geometry, config, false);
-
-//    if ((config->GetnMarker_Analyze() != 0) && compressible)
-//      output->SpecialOutput_Distortion(solver[FLOW_SOL], geometry, config, false);
-
-//    if (config->GetnMarker_NearFieldBound() != 0)
-//      output->SpecialOutput_SonicBoom(solver[FLOW_SOL], geometry, config, false);
-
-//    if (config->GetPlot_Section_Forces())
-//      output->SpecialOutput_SpanLoad(solver[FLOW_SOL], geometry, config, false);
-
     /*--- Surface based obj. function ---*/
 
     solver[FLOW_SOL]->Evaluate_ObjFunc(config);
diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp
index 314f21daace5..ea5af8b59457 100644
--- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp
@@ -139,14 +139,11 @@ CMultizoneDriver::~CMultizoneDriver(void) {
 void CMultizoneDriver::StartSolver() {
 
   /*--- Find out the minimum of all references times and then set each zone to this (same) value.
-   * (To ensure that all zones run synchronously in time, be it a dimensional or non-dimensionalized one.) ---*/
+        To ensure that all zones run synchronously in time, be it a dimensional or non-dimensionalized one. ---*/
 
-  su2double Time_Ref = config_container[ZONE_0]->GetTime_Ref();
-
-  for (iZone = 1; iZone < nZone; iZone++) {
-    if (config_container[iZone]->GetTime_Ref() < Time_Ref)
-      Time_Ref = config_container[iZone]->GetTime_Ref();
-  }
+  su2double Time_Ref = std::numeric_limits<su2double>::max();
+  for (iZone = 0; iZone < nZone; iZone++)
+    Time_Ref = min(Time_Ref, config_container[iZone]->GetTime_Ref());
 
   for (iZone = 0; iZone < nZone; iZone++) {
 
diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
index 50baa8b386de..9996f85c410f 100644
--- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
+++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
@@ -40,90 +40,62 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
                                                 const CConfig *donor_config, unsigned long Marker_Donor,
                                                 unsigned long Vertex_Donor, unsigned long Point_Donor) {
 
+  /*--- Compute distance of donor point to PointNormal for T-gradient/heatflux computation ---*/
   const auto nDim = donor_geometry->GetnDim();
-
-  su2double Twall, Tnormal, dTdn, rho_cp_solid, Prandtl_Lam, laminar_viscosity,
-      thermal_diffusivity, thermal_conductivity=0.0, thermal_conductivityND,
-      heat_flux_density=0.0, conductivity_over_dist=0.0;
-
-  /*--- Check whether the current zone is a solid zone or a fluid zone ---*/
-
-  const bool compressible_flow = (donor_config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE);
-  const bool incompressible_flow = (donor_config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) && donor_config->GetEnergy_Equation();
-  const bool heat_equation = (donor_config->GetKind_Solver() == HEAT_EQUATION) ||
-                             (donor_config->GetKind_Solver() == DISC_ADJ_HEAT);
-
   const auto Coord = donor_geometry->nodes->GetCoord(Point_Donor);
   const auto PointNormal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal_Neighbor();
   const auto Coord_Normal = donor_geometry->nodes->GetCoord(PointNormal);
 
-  Twall = 0.0; Tnormal = 0.0; dTdn = 0.0;
-
-  su2double Edge_Vector[3] = {0.0};
+  su2double Edge_Vector[MAXNDIM] = {0.0};
   GeometryToolbox::Distance(nDim, Coord_Normal, Coord, Edge_Vector);
-  su2double dist = GeometryToolbox::Norm(nDim, Edge_Vector);
+  const su2double dist = GeometryToolbox::Norm(nDim, Edge_Vector);
 
   /*--- Retrieve temperature solution and its gradient ---*/
 
-  if (compressible_flow) {
-
-    Twall   = donor_solution->GetNodes()->GetPrimitive(Point_Donor,0);
-    Tnormal = donor_solution->GetNodes()->GetPrimitive(PointNormal,0);
-
-    dTdn = (Twall - Tnormal)/dist;
-  }
-  else if (incompressible_flow) {
-
-    Twall   = donor_solution->GetNodes()->GetTemperature(Point_Donor);
-    Tnormal = donor_solution->GetNodes()->GetTemperature(PointNormal);
+  const su2double Twall   = donor_solution->GetNodes()->GetTemperature(Point_Donor);
+  const su2double Tnormal = donor_solution->GetNodes()->GetTemperature(PointNormal);
 
-    dTdn = (Twall - Tnormal)/dist;
-  }
-  else if (heat_equation) {
-    Twall   = donor_solution->GetNodes()->GetSolution(Point_Donor,0);
-    Tnormal = donor_solution->GetNodes()->GetSolution(PointNormal,0);
-
-    // TODO: Check if these improve accuracy, if needed at all
-    //    const auto Normal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal();
-    //    su2double Area = GeometryToolbox::Norm(nDim, Normal);
-    //    for (iDim = 0; iDim < nDim; iDim++) {
-    //      dTdn += (Twall - Tnormal)/dist * (Edge_Vector[iDim]/dist) * (Normal[iDim]/Area);
-    //    }
-
-    dTdn = (Twall - Tnormal)/dist;
-  }
-  else {
+  const su2double dTdn = (Twall - Tnormal)/dist;
+  // TODO: Check if these improve accuracy, if needed at all
+  //    const auto Normal = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNormal();
+  //    su2double Area = GeometryToolbox::Norm(nDim, Normal);
+  //    for (iDim = 0; iDim < nDim; iDim++) {
+  //      dTdn += (Twall - Tnormal)/dist * (Edge_Vector[iDim]/dist) * (Normal[iDim]/Area);
+  //    }
 
-    SU2_MPI::Error("Transfer of conjugate heat variables failed (non-supported donor solver).", CURRENT_FUNCTION);
-  }
+  /*--- Calculate the heat flux density (temperature gradient times thermal conductivity) and
+        thermal conductivity divided by distance. ---*/
+  su2double thermal_conductivity = 0.0;
+  su2double heat_flux_density = 0.0;
+  su2double conductivity_over_dist = 0.0;
 
-  /*--- Calculate the heat flux density (temperature gradient times thermal conductivity) ---*/
+  const bool compressible_flow = (donor_config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE);
+  const bool incompressible_flow = (donor_config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) && donor_config->GetEnergy_Equation();
 
   if (compressible_flow) {
 
-    su2double Gamma         = donor_config->GetGamma();
-    su2double Gas_Constant  = donor_config->GetGas_ConstantND();
-    su2double Cp            = (Gamma / (Gamma - 1.0)) * Gas_Constant;
+    const su2double Gamma = donor_config->GetGamma();
+    const su2double Gas_Constant = donor_config->GetGas_ConstantND();
+    const su2double Cp = (Gamma / (Gamma - 1.0)) * Gas_Constant;
 
-    Prandtl_Lam             = donor_config->GetPrandtl_Lam();
-    laminar_viscosity       = donor_solution->GetNodes()->GetLaminarViscosity(Point_Donor); // TDE check for consistency
-    Cp                      = (Gamma / (Gamma - 1.0)) * Gas_Constant;
+    const su2double Prandtl_Lam = donor_config->GetPrandtl_Lam();
+    const su2double laminar_viscosity = donor_solution->GetNodes()->GetLaminarViscosity(Point_Donor); // TDE check for consistency
 
-    thermal_conductivityND  = Cp*(laminar_viscosity/Prandtl_Lam);
-    heat_flux_density       = thermal_conductivityND*dTdn;
+    const su2double thermal_conductivityND = Cp*(laminar_viscosity/Prandtl_Lam);
+    heat_flux_density = thermal_conductivityND*dTdn;
 
     if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
         (donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
-      thermal_conductivity    = thermal_conductivityND*donor_config->GetViscosity_Ref();
-      conductivity_over_dist  = thermal_conductivity/dist;
+      thermal_conductivity   = thermal_conductivityND*donor_config->GetViscosity_Ref();
+      conductivity_over_dist = thermal_conductivity/dist;
     }
   }
   else if (incompressible_flow) {
 
     const auto iPoint = donor_geometry->vertex[Marker_Donor][Vertex_Donor]->GetNode();
 
-    thermal_conductivityND  = donor_solution->GetNodes()->GetThermalConductivity(iPoint);
+    const su2double thermal_conductivityND  = donor_solution->GetNodes()->GetThermalConductivity(iPoint);
     heat_flux_density       = thermal_conductivityND*dTdn;
 
     if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
@@ -140,24 +112,25 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
                                  *donor_config->GetViscosity_Ref();
           break;
 
-        default:
+        case CONDUCTIVITYMODEL::POLYNOMIAL:
+          SU2_MPI::Error("Polynomial Conductivity model not implemented for CHT interface.", CURRENT_FUNCTION);
           break;
       }
 
       conductivity_over_dist  = thermal_conductivity/dist;
     }
   }
-  else if (heat_equation) {
+  else if (donor_config->GetHeatProblem()) {
 
     /*--- Heat solver stand-alone case ---*/
 
-    thermal_diffusivity     = donor_config->GetThermalDiffusivity_Solid();
-    heat_flux_density       = thermal_diffusivity*dTdn;
+    const su2double thermal_diffusivity = donor_config->GetThermalDiffusivity_Solid();
+    heat_flux_density = thermal_diffusivity*dTdn;
 
     if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
         (donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
-      rho_cp_solid            = donor_config->GetSpecific_Heat_Cp()*donor_config->GetDensity_Solid();
+      const su2double rho_cp_solid = donor_config->GetSpecific_Heat_Cp()*donor_config->GetDensity_Solid();
       conductivity_over_dist  = thermal_diffusivity*rho_cp_solid/dist;
     }
   }
diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build
index ee288818ef4e..736100a1b812 100644
--- a/SU2_CFD/src/meson.build
+++ b/SU2_CFD/src/meson.build
@@ -21,6 +21,7 @@ su2_cfd_src += files(['output/COutputFactory.cpp',
                       'output/CElasticityOutput.cpp',
                       'output/CFlowCompOutput.cpp',
                       'output/CFlowCompFEMOutput.cpp',
+                      'output/CFVMOutput.cpp',
                       'output/CFlowOutput.cpp',
                       'output/CFlowIncOutput.cpp',
                       'output/CHeatOutput.cpp',
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index b7e583757d47..da7512418d23 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -713,7 +713,7 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
 
     residual[nDim+1] = Volume * scalar_factor * dot_product;
 
-    /*--- If a RANS turbulence model ias used an additional source term, based on the eddy viscosity gradient is added. ---*/
+    /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
     if(turbulent) {
 
       /*--- Compute a scalar factor ---*/
diff --git a/SU2_CFD/src/output/CAdjFlowCompOutput.cpp b/SU2_CFD/src/output/CAdjFlowCompOutput.cpp
index 6d5db2133b20..2cf99b95e396 100644
--- a/SU2_CFD/src/output/CAdjFlowCompOutput.cpp
+++ b/SU2_CFD/src/output/CAdjFlowCompOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_adj_flow_comp.cpp
+ * \file CAdjFlowCompOutput.cpp
  * \brief Main subroutines for flow discrete adjoint output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
@@ -26,7 +26,7 @@
  */
 
 
-#include "../../include/output/CAdjFlowOutput.hpp"
+#include "../../include/output/CAdjFlowCompOutput.hpp"
 
 #include "../../../Common/include/geometry/CGeometry.hpp"
 #include "../../include/solvers/CSolver.hpp"
diff --git a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp
index f51f0e6c7959..8bcec5adeb9c 100644
--- a/SU2_CFD/src/output/CAdjFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CAdjFlowIncOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_adj_flow_inc.cpp
+ * \file CAdjFlowIncOutput.cpp
  * \brief Main subroutines for flow discrete adjoint output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CAdjHeatOutput.cpp b/SU2_CFD/src/output/CAdjHeatOutput.cpp
index 3c3cbd0cabdf..0570f97ad680 100644
--- a/SU2_CFD/src/output/CAdjHeatOutput.cpp
+++ b/SU2_CFD/src/output/CAdjHeatOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_adj_heat.cpp
+ * \file CAdjHeatOutput.cpp
  * \brief Main subroutines for flow discrete adjoint output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CBaselineOutput.cpp b/SU2_CFD/src/output/CBaselineOutput.cpp
index a37004c092f8..253ac5d432bc 100644
--- a/SU2_CFD/src/output/CBaselineOutput.cpp
+++ b/SU2_CFD/src/output/CBaselineOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_baseline.cpp
+ * \file CBaselineOutput.cpp
  * \brief Main subroutines for flow discrete adjoint output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp
index 55d2f9f84a08..c48dab833dba 100644
--- a/SU2_CFD/src/output/CElasticityOutput.cpp
+++ b/SU2_CFD/src/output/CElasticityOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_elasticity.cpp
+ * \file CElasticityOutput.cpp
  * \brief Main subroutines for FEA output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CFVMOutput.cpp b/SU2_CFD/src/output/CFVMOutput.cpp
new file mode 100644
index 000000000000..364f1f357eba
--- /dev/null
+++ b/SU2_CFD/src/output/CFVMOutput.cpp
@@ -0,0 +1,79 @@
+/*!
+ * \file CFVMOutput.cpp
+ * \brief Main subroutines for Finite Volume Method output
+ * \author T. Kattmann
+ * \version 7.1.1 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2021, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "../../include/output/CFVMOutput.hpp"
+#include "../../../Common/include/geometry/CGeometry.hpp"
+
+
+CFVMOutput::CFVMOutput(CConfig *config, unsigned short nDim, bool fem_output) : COutput (config, nDim, fem_output){ }
+
+void CFVMOutput::AddCoordinates() {
+
+  // Grid coordinates
+  AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
+  AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
+  if (nDim == 3)
+    AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
+}
+
+
+void CFVMOutput::AddCommonFVMOutputs(const CConfig *config) {
+
+  // Mesh quality metrics
+  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
+  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
+  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
+
+  AddVolumeOutput("RANK", "rank", "MPI", "Rank of the MPI-partition");
+
+  for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
+    stringstream key, name;
+    key << "MG_" << iMesh;
+    name << "Coarse_Grid_" << iMesh;
+    AddVolumeOutput(key.str(), name.str(), "MULTIGRID", "Coarse mesh");
+  }
+}
+
+void CFVMOutput::LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint) {
+
+  // Mesh quality metrics, computed in CPhysicalGeometry::ComputeMeshQualityStatistics.
+  if (config->GetWrt_MeshQuality()) {
+    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
+    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
+    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
+  }
+
+  SetVolumeOutputValue("RANK", iPoint, rank);
+
+  if (config->GetWrt_MultiGrid()) {
+    for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
+      stringstream key;
+      key << "MG_" << iMesh;
+      SetVolumeOutputValue(key.str(), iPoint, geometry->CoarseGridColor(iPoint,iMesh-1));
+    }
+  }
+}
diff --git a/SU2_CFD/src/output/CFlowCompFEMOutput.cpp b/SU2_CFD/src/output/CFlowCompFEMOutput.cpp
index 91cade8a86c9..5eb60cbe660d 100644
--- a/SU2_CFD/src/output/CFlowCompFEMOutput.cpp
+++ b/SU2_CFD/src/output/CFlowCompFEMOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_flow_comp_fem.cpp
+ * \file CFlowCompFEMOutput.cpp
  * \brief Main subroutines for compressible flow output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp
index 5f053a80275f..0126222cf4ec 100644
--- a/SU2_CFD/src/output/CFlowCompOutput.cpp
+++ b/SU2_CFD/src/output/CFlowCompOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_flow_comp.cpp
+ * \file CFlowCompOutput.cpp
  * \brief Main subroutines for compressible flow output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
@@ -291,10 +291,7 @@ void CFlowCompOutput::SetHistoryOutputFields(CConfig *config){
 void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){
 
   // Grid coordinates
-  AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
-  AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
-  if (nDim == 3)
-    AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
+  AddCoordinates();
 
   // Solution variables
   AddVolumeOutput("DENSITY",    "Density",    "SOLUTION", "Density");
@@ -433,17 +430,13 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
 
   CVariable* Node_Flow = solver[FLOW_SOL]->GetNodes();
   CVariable* Node_Turb = nullptr;
+  const auto Node_Geo  = geometry->nodes;
 
   if (config->GetKind_Turb_Model() != NONE){
     Node_Turb = solver[TURB_SOL]->GetNodes();
   }
 
-  CPoint*    Node_Geo  = geometry->nodes;
-
-  SetVolumeOutputValue("COORD-X", iPoint,  Node_Geo->GetCoord(iPoint, 0));
-  SetVolumeOutputValue("COORD-Y", iPoint,  Node_Geo->GetCoord(iPoint, 1));
-  if (nDim == 3)
-    SetVolumeOutputValue("COORD-Z", iPoint, Node_Geo->GetCoord(iPoint, 2));
+  LoadCoordinates(Node_Geo->GetCoord(iPoint), iPoint);
 
   SetVolumeOutputValue("DENSITY",    iPoint, Node_Flow->GetSolution(iPoint, 0));
   SetVolumeOutputValue("MOMENTUM-X", iPoint, Node_Flow->GetSolution(iPoint, 1));
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
index 229d61f866c5..57398f448e82 100644
--- a/SU2_CFD/src/output/CFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_flow_inc.cpp
+ * \file CFlowIncOutput.cpp
  * \brief Main subroutines for incompressible flow output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
@@ -373,10 +373,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
 void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
 
   // Grid coordinates
-  AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
-  AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
-  if (nDim == 3)
-    AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
+  AddCoordinates();
 
   // SOLUTION variables
   AddVolumeOutput("PRESSURE",   "Pressure",   "SOLUTION", "Pressure");
@@ -521,6 +518,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   CVariable* Node_Heat = nullptr;
   CVariable* Node_Turb = nullptr;
   CVariable* Node_Rad = nullptr;
+  const auto Node_Geo = geometry->nodes;
 
   if (config->GetKind_Turb_Model() != NONE){
     Node_Turb = solver[TURB_SOL]->GetNodes();
@@ -529,12 +527,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
     Node_Heat = solver[HEAT_SOL]->GetNodes();
   }
 
-  CPoint*    Node_Geo  = geometry->nodes;
-
-  SetVolumeOutputValue("COORD-X", iPoint,  Node_Geo->GetCoord(iPoint, 0));
-  SetVolumeOutputValue("COORD-Y", iPoint,  Node_Geo->GetCoord(iPoint, 1));
-  if (nDim == 3)
-    SetVolumeOutputValue("COORD-Z", iPoint, Node_Geo->GetCoord(iPoint, 2));
+  LoadCoordinates(Node_Geo->GetCoord(iPoint), iPoint);
 
   SetVolumeOutputValue("PRESSURE",   iPoint, Node_Flow->GetSolution(iPoint, 0));
   SetVolumeOutputValue("VELOCITY-X", iPoint, Node_Flow->GetSolution(iPoint, 1));
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 815b400b13a4..27d5d78fec48 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_flow.cpp
+ * \file CFlowOutput.cpp
  * \brief Main subroutines for compressible flow output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
@@ -30,11 +30,10 @@
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../../include/solvers/CSolver.hpp"
 
-CFlowOutput::CFlowOutput(CConfig *config, unsigned short nDim, bool fem_output) : COutput (config, nDim, fem_output){
+CFlowOutput::CFlowOutput(CConfig *config, unsigned short nDim, bool fem_output) : CFVMOutput (config, nDim, fem_output){
 
 }
 
-
 CFlowOutput::~CFlowOutput(void){}
 
 void CFlowOutput::AddAnalyzeSurfaceOutput(CConfig *config){
@@ -2735,41 +2734,6 @@ bool CFlowOutput::WriteVolume_Output(CConfig *config, unsigned long Iter, bool f
   return false || force_writing;
 }
 
-void CFlowOutput::AddCommonFVMOutputs(const CConfig *config) {
-
-  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
-  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
-  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
-
-  AddVolumeOutput("RANK", "rank", "MPI", "Rank of the MPI-partition");
-
-  for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
-    stringstream key, name;
-    key << "MG_" << iMesh;
-    name << "Coarse_Grid_" << iMesh;
-    AddVolumeOutput(key.str(), name.str(), "MULTIGRID", "Coarse mesh");
-  }
-}
-
-void CFlowOutput::LoadCommonFVMOutputs(const CConfig* config, const CGeometry* geometry, unsigned long iPoint) {
-
-  if (config->GetWrt_MeshQuality()) {
-    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
-    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
-    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
-  }
-
-  SetVolumeOutputValue("RANK", iPoint, rank);
-
-  if (config->GetWrt_MultiGrid()) {
-    for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); ++iMesh) {
-      stringstream key;
-      key << "MG_" << iMesh;
-      SetVolumeOutputValue(key.str(), iPoint, geometry->CoarseGridColor(iPoint,iMesh-1));
-    }
-  }
-}
-
 void CFlowOutput::SetTimeAveragedFields(){
   AddVolumeOutput("MEAN_DENSITY", "MeanDensity", "TIME_AVERAGE", "Mean density");
   AddVolumeOutput("MEAN_VELOCITY-X", "MeanVelocity_x", "TIME_AVERAGE", "Mean velocity x-component");
diff --git a/SU2_CFD/src/output/CHeatOutput.cpp b/SU2_CFD/src/output/CHeatOutput.cpp
index 831c2afa98b1..0e277da57877 100644
--- a/SU2_CFD/src/output/CHeatOutput.cpp
+++ b/SU2_CFD/src/output/CHeatOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_heat.cpp
+ * \file CHeatOutput.cpp
  * \brief Main subroutines for the heat solver output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
@@ -30,7 +30,7 @@
 #include "../../../Common/include/geometry/CGeometry.hpp"
 #include "../../include/solvers/CSolver.hpp"
 
-CHeatOutput::CHeatOutput(CConfig *config, unsigned short nDim) : COutput(config, nDim, false) {
+CHeatOutput::CHeatOutput(CConfig *config, unsigned short nDim) : CFVMOutput(config, nDim, false) {
 
   multiZone = config->GetMultizone_Problem();
 
@@ -117,10 +117,7 @@ void CHeatOutput::SetHistoryOutputFields(CConfig *config){
 void CHeatOutput::SetVolumeOutputFields(CConfig *config){
 
   // Grid coordinates
-  AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
-  AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
-  if (nDim == 3)
-    AddVolumeOutput("COORD-Z", "z", "COORDINATES","z-component of the coordinate vector");
+  AddCoordinates();
 
   // SOLUTION
   AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature");
@@ -131,26 +128,17 @@ void CHeatOutput::SetVolumeOutputFields(CConfig *config){
   // Residuals
   AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature");
 
-  // Mesh quality metrics, computed in CPhysicalGeometry::ComputeMeshQualityStatistics.
-  AddVolumeOutput("ORTHOGONALITY", "Orthogonality", "MESH_QUALITY", "Orthogonality Angle (deg.)");
-  AddVolumeOutput("ASPECT_RATIO",  "Aspect_Ratio",  "MESH_QUALITY", "CV Face Area Aspect Ratio");
-  AddVolumeOutput("VOLUME_RATIO",  "Volume_Ratio",  "MESH_QUALITY", "CV Sub-Volume Ratio");
-
-  // MPI-Rank
-  AddVolumeOutput("RANK", "rank", "MPI", "Rank of the MPI-partition");
+  AddCommonFVMOutputs(config);
 }
 
 
 void CHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
 
   CVariable* Node_Heat = solver[HEAT_SOL]->GetNodes();
-  CPoint*    Node_Geo  = geometry->nodes;
+  const auto Node_Geo  = geometry->nodes;
 
   // Grid coordinates
-  SetVolumeOutputValue("COORD-X", iPoint,  Node_Geo->GetCoord(iPoint, 0));
-  SetVolumeOutputValue("COORD-Y", iPoint,  Node_Geo->GetCoord(iPoint, 1));
-  if (nDim == 3)
-    SetVolumeOutputValue("COORD-Z", iPoint, Node_Geo->GetCoord(iPoint, 2));
+  LoadCoordinates(Node_Geo->GetCoord(iPoint), iPoint);
 
   // SOLUTION
   SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0));
@@ -158,15 +146,7 @@ void CHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver *
   // Residuals
   SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[HEAT_SOL]->LinSysRes(iPoint, 0));
 
-  // Mesh quality metrics
-  if (config->GetWrt_MeshQuality()) {
-    SetVolumeOutputValue("ORTHOGONALITY", iPoint, geometry->Orthogonality[iPoint]);
-    SetVolumeOutputValue("ASPECT_RATIO",  iPoint, geometry->Aspect_Ratio[iPoint]);
-    SetVolumeOutputValue("VOLUME_RATIO",  iPoint, geometry->Volume_Ratio[iPoint]);
-  }
-
-  // MPI-Rank
-  SetVolumeOutputValue("RANK", iPoint, rank);
+  LoadCommonFVMOutputs(config, geometry, iPoint);
 }
 
 void CHeatOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex){
diff --git a/SU2_CFD/src/output/CMeshOutput.cpp b/SU2_CFD/src/output/CMeshOutput.cpp
index ea2969ecbd29..5b27d3e7bbee 100644
--- a/SU2_CFD/src/output/CMeshOutput.cpp
+++ b/SU2_CFD/src/output/CMeshOutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_mesh.cpp
+ * \file CMeshOutput.cpp
  * \brief Main subroutines for the heat solver output
  * \author R. Sanchez
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/CNEMOCompOutput.cpp b/SU2_CFD/src/output/CNEMOCompOutput.cpp
index 291127f4e08b..779b72aa86af 100644
--- a/SU2_CFD/src/output/CNEMOCompOutput.cpp
+++ b/SU2_CFD/src/output/CNEMOCompOutput.cpp
@@ -296,10 +296,7 @@ void CNEMOCompOutput::SetVolumeOutputFields(CConfig *config){
   unsigned short nSpecies = config->GetnSpecies();
 
   // Grid coordinates
-  AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
-  AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
-  if (nDim == 3)
-    AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
+  AddCoordinates();
 
   // Solution variables
   for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
@@ -439,17 +436,13 @@ void CNEMOCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
   CVariable* Node_Flow = solver[FLOW_SOL]->GetNodes();
   CVariable* Node_Turb = NULL;
   unsigned short nSpecies = config->GetnSpecies();
+  const auto Node_Geo = geometry->nodes;
 
   if (config->GetKind_Turb_Model() != NONE){
     Node_Turb = solver[TURB_SOL]->GetNodes();
   }
 
-  auto*    Node_Geo  = geometry->nodes;
-
-  SetVolumeOutputValue("COORD-X", iPoint,  Node_Geo->GetCoord(iPoint, 0));
-  SetVolumeOutputValue("COORD-Y", iPoint,  Node_Geo->GetCoord(iPoint, 1));
-  if (nDim == 3)
-    SetVolumeOutputValue("COORD-Z", iPoint, Node_Geo->GetCoord(iPoint, 2));
+  LoadCoordinates(Node_Geo->GetCoord(iPoint), iPoint);
 
   for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     SetVolumeOutputValue("DENSITY_" + std::to_string(iSpecies),   iPoint, Node_Flow->GetSolution(iPoint, iSpecies));
diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp
index efc29fe036eb..fe736fbccb99 100644
--- a/SU2_CFD/src/output/COutput.cpp
+++ b/SU2_CFD/src/output/COutput.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_structure.cpp
+ * \file COutput.cpp
  * \brief Main subroutines for output solver information
  * \author F. Palacios, T. Economon
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/output/COutputFactory.cpp b/SU2_CFD/src/output/COutputFactory.cpp
index 2feb4cd3ac3d..5814c85767b6 100644
--- a/SU2_CFD/src/output/COutputFactory.cpp
+++ b/SU2_CFD/src/output/COutputFactory.cpp
@@ -33,7 +33,7 @@
 #include "../../include/output/CAdjElasticityOutput.hpp"
 #include "../../include/output/CFlowCompOutput.hpp"
 #include "../../include/output/CNEMOCompOutput.hpp"
-#include "../../include/output/CAdjFlowOutput.hpp"
+#include "../../include/output/CAdjFlowCompOutput.hpp"
 #include "../../include/output/CFlowCompFEMOutput.hpp"
 #include "../../include/output/CFlowIncOutput.hpp"
 #include "../../include/output/CAdjFlowIncOutput.hpp"
diff --git a/SU2_CFD/src/output/output_structure_legacy.cpp b/SU2_CFD/src/output/output_structure_legacy.cpp
index acfa5d7f1b5f..fc65caa3aeb6 100644
--- a/SU2_CFD/src/output/output_structure_legacy.cpp
+++ b/SU2_CFD/src/output/output_structure_legacy.cpp
@@ -1,5 +1,5 @@
 /*!
- * \file output_structure.cpp
+ * \file output_structure_legacy.cpp
  * \brief Main subroutines for output solver information
  * \author F. Palacios, T. Economon
  * \version 7.1.1 "Blackbird"
diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp
index f8a8b8620333..e31233ae6e46 100644
--- a/SU2_CFD/src/solvers/CHeatSolver.cpp
+++ b/SU2_CFD/src/solvers/CHeatSolver.cpp
@@ -1264,9 +1264,6 @@ void CHeatSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
 
 void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) {
 
-  unsigned long total_index;
-  su2double Delta, Vol, *local_Res_TruncError;
-
   /*--- Set maximum residual to zero ---*/
 
   SetResToZero();
@@ -1277,29 +1274,20 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
 
     /*--- Read the residual ---*/
 
-    local_Res_TruncError = nodes->GetResTruncError(iPoint);
-
-    /*--- Read the volume ---*/
-
-    Vol = geometry->nodes->GetVolume(iPoint);
+    su2double* local_Res_TruncError = nodes->GetResTruncError(iPoint);
 
     /*--- Modify matrix diagonal to assure diagonal dominance ---*/
 
     if (nodes->GetDelta_Time(iPoint) != 0.0) {
 
-      if(flow) {
-        Delta = Vol / nodes->GetDelta_Time(iPoint);
-        Jacobian.AddVal2Diag(iPoint, Delta);
-      }
-      else if (heat_equation) {
-        Delta = Vol / nodes->GetDelta_Time(iPoint);
-        Jacobian.AddVal2Diag(iPoint, Delta);
-      }
+      // Identical for flow and heat
+      const su2double Delta = geometry->nodes->GetVolume(iPoint) / nodes->GetDelta_Time(iPoint);
+      Jacobian.AddVal2Diag(iPoint, Delta);
 
     } else {
       Jacobian.SetVal2Diag(iPoint, 1.0);
       for (auto iVar = 0u; iVar < nVar; iVar++) {
-        total_index = iPoint*nVar + iVar;
+        const auto total_index = iPoint*nVar + iVar;
         LinSysRes[total_index] = 0.0;
         local_Res_TruncError[iVar] = 0.0;
       }
@@ -1308,7 +1296,7 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
     /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/
 
     for (auto iVar = 0; iVar < nVar; iVar++) {
-      total_index = iPoint*nVar+iVar;
+      const auto total_index = iPoint*nVar+iVar;
       LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]);
       LinSysSol[total_index] = 0.0;
       Residual_RMS[iVar] += LinSysRes[total_index]*LinSysRes[total_index];
@@ -1320,7 +1308,7 @@ void CHeatSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
 
   for (auto iPoint = nPointDomain; iPoint < nPoint; iPoint++) {
     for (auto iVar = 0u; iVar < nVar; iVar++) {
-      total_index = iPoint*nVar + iVar;
+      const auto total_index = iPoint*nVar + iVar;
       LinSysRes[total_index] = 0.0;
       LinSysSol[total_index] = 0.0;
     }
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index f1d6d54963e9..6a564e5f1118 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -131,9 +131,11 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
   for (auto iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
 
-    /*--- Only "outlet"/donor periodic marker ---*/
+    /*--- Only "outlet"/donor periodic marker. GetnMarker_Periodic counts from 1: First all the inlet marker from
+          the periodic marker pairs and then the all the outlets. I.e. in order to get the outlet of the first pair
+          we need to divide the number of periodic markers by 2 and add 1 (because count starts at 1). ---*/
     if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY &&
-        config->GetMarker_All_PerBound(iMarker) == 2) {
+        config->GetMarker_All_PerBound(iMarker) == (config->GetnMarker_Periodic()/2 + 1)) {
 
       for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {
 
@@ -145,7 +147,7 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
           const auto AreaNormal = geometry->vertex[iMarker][iVertex]->GetNormal();
 
-          auto FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
+          const su2double FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
 
           /*--- m_dot = dot_prod(n*v) * A * rho, with n beeing unit normal. ---*/
           MassFlow_Local += nodes->GetProjVel(iPoint, AreaNormal) * nodes->GetDensity(iPoint);
@@ -154,7 +156,7 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
           Average_Density_Local += FaceArea * nodes->GetDensity(iPoint);
 
-          /*--- Due to periodicty, temperatures are equal one the inlet(1) and outlet(2) ---*/
+          /*--- Due to periodicity, temperatures are equal one the inlet(1) and outlet(2) ---*/
           Temperature_Local += FaceArea * nodes->GetTemperature(iPoint);
 
         } // if domain
@@ -230,7 +232,7 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
         /*--- Identify the boundary by string name and retrive heatflux from config ---*/
         const auto Marker_StringTag = config->GetMarker_All_TagBound(iMarker);
-        const auto Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
+        const su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
 
         for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {
 
@@ -240,7 +242,7 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 
           const auto AreaNormal = geometry->vertex[iMarker][iVertex]->GetNormal();
 
-          auto FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
+          const su2double FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
 
           HeatFlow_Local += FaceArea * (-1.0) * Wall_HeatFlux/config->GetHeat_Flux_Ref();;
         } // loop Vertices
@@ -263,7 +265,7 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config,
   const auto InnerIter = config->GetInnerIter();
 
   /*--- Reference node on inlet periodic marker to compute relative distance along periodic translation vector. ---*/
-  const su2double* ReferenceNode = geometry->GetStreamwise_Periodic_RefNode();
+  const auto ReferenceNode = geometry->GetStreamwise_Periodic_RefNode();
 
   /*--- Compute square of the distance between the 2 periodic surfaces. ---*/
   const su2double norm2_translation = GeometryToolbox::SquaredNorm(nDim, config->GetPeriodic_Translation(0));
@@ -445,12 +447,12 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
     }
     else { // ISOTHERMAL
 
-      auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+      const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
 
       /*--- Get coordinates of i & nearest normal and compute distance ---*/
 
-      auto Coord_i = geometry->nodes->GetCoord(iPoint);
-      auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
+      const auto Coord_i = geometry->nodes->GetCoord(iPoint);
+      const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
       su2double Edge_Vector[MAXNDIM];
       GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector);
       su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector);
@@ -555,12 +557,12 @@ void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **sol
 
       /*--- Compute closest normal neighbor ---*/
 
-      auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+      const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
 
       /*--- Get coordinates of i & nearest normal and compute distance ---*/
 
-      auto Coord_i = geometry->nodes->GetCoord(iPoint);
-      auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
+      const auto Coord_i = geometry->nodes->GetCoord(iPoint);
+      const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
       su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_j, Coord_i);
 
       /*--- Compute wall temperature from both temperatures ---*/
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
index 6821416134a5..a32ece5723b8 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
@@ -44,7 +44,7 @@ SCREEN_WRT_FREQ_OUTER= 50
 %
 HISTORY_OUTPUT= ( ITER, BGS_RES[0], RMS_RES[0], BGS_RES[1], RMS_RES[1],\
  FLOW_COEFF[0], HEAT[0], AERO_COEFF[0], HEAT[1],\
- LINSOL[0], LINSOL[1])
+ LINSOL[0], LINSOL[1], TAVG_HEAT[1])
 %
 OUTPUT_FILES= (RESTART, PARAVIEW)
 OUTPUT_WRT_FREQ= 1
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py
index 003f742ce090..9751520d2be4 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py
@@ -80,7 +80,7 @@
 adjoint.addExpected("solution_adj_avtp_1_00053.dat")
 
 # gradient projection
-dot = adjont = ExternalRun("DOT",dot_ad_command,True)
+dot = ExternalRun("DOT",dot_ad_command,True)
 dot.setMaxTries(max_tries)
 dot.addConfig("chtMaster.cfg")
 dot.addData("fluid.cfg") # zonal cfg's can be symlinked
@@ -106,6 +106,7 @@
 
 # Driver --------------------------------------------------------------- #
 
+# The input variable is the constraint tolerance which is not used for our purpose of finite differences
 driver = ExteriorPenaltyDriver(0.005)
 driver.addObjective("min", tavgT)
 
@@ -120,6 +121,7 @@
 
 # Primal simulation for each deformed DV
 for iLoop in range(0, nDV, 1):
+    print("Computing deformed primal ", iLoop, "/", nDV-1)
     x = driver.getInitial()
     x[iLoop] = 1e-4 # DV_VALUE
     driver.fun(x)
@@ -127,10 +129,12 @@
 
 # Undeformed/initial primal last in order to have the correct solution in
 # the WorkindDirectory for the following adjoint
+print("Computing baseline primal")
 x = driver.getInitial()
 driver.fun(x) # baseline evaluation
 
 # Compute discrete adjoint gradient
+print("Computing discrete adjoint gradient")
 driver.grad(x)
 
 his.close()
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
index 8e59310c1a94..bcbe8a7ba2d0 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
@@ -10,17 +10,30 @@
 FDvals = pd.read_csv("doe.his")
 
 # additional values
+nDV = 18
 FDstep = 1e-4
 FDstring = '  tavgT'
 DAstring = 'AVG_TEMPERATURE gradient '
 
 # create FD gradient
-FDgrad = (FDvals[FDstring].iloc[:18] - FDvals[FDstring].iloc[18]) / FDstep
+FDgrad = (FDvals[FDstring].iloc[:nDV] - FDvals[FDstring].iloc[nDV]) / FDstep
 
 # absolute difference
 absoluteDiff = DAgrad[DAstring] - FDgrad
-print("DAgrad - FDgrad\n", absoluteDiff)
+print("absolute diff = DAgrad - FDgrad")
 
 # relative difference in percent
-relDiffPercent = abs(DAgrad[DAstring] - FDgrad)/DAgrad[DAstring] * 100
-print("(DAgrad - FDgrad) / DAgrad * 100\n", relDiffPercent)
+relDiffPercent = (DAgrad[DAstring] - FDgrad)/abs(DAgrad[DAstring]) * 100
+print("relative diff = (DAgrad - FDgrad) / abs(DAgrad) * 100")
+
+
+print('')
+print('+-----------+-------------------+-------------------+-------------------+-------------------+')
+print('| DV number |       DA gradient |       FD gradient |     absolute diff | relative diff [%] |')
+print('+-----------+-------------------+-------------------+-------------------+-------------------+')
+
+for i_dv in range(0, nDV,1):
+    print('|{0:10d} |{1:18.10f} |{2:18.10f} |{3:18.10f} |{4:18.10f} |'.format(i_dv, DAgrad[DAstring].iloc[i_dv], FDgrad[i_dv], absoluteDiff[i_dv], relDiffPercent[i_dv]))
+
+print('+-----------+-------------------+-------------------+-------------------+-------------------+')
+print('')
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 54728bb237a0..77faa7d3ecf9 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1252,19 +1252,6 @@ def main():
     cht_incompressible.tol       = 0.00001
     test_list.append(cht_incompressible)
 
-    # CHT incompressible unsteady
-    cht_incompressible_unsteady           = TestCase('cht_incompressible_unsteady')
-    cht_incompressible_unsteady.cfg_dir   = "coupled_cht/incomp_2d_unsteady"
-    cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
-    cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-1.305471, -0.080372, -0.080376, -0.080372] #last 4 columns
-    cht_incompressible_unsteady.su2_exec  = "mpirun -n 2 SU2_CFD"
-    cht_incompressible_unsteady.timeout   = 1600
-    cht_incompressible_unsteady.multizone = True
-    cht_incompressible_unsteady.unsteady  = True
-    cht_incompressible_unsteady.tol       = 0.00001
-    test_list.append(cht_incompressible_unsteady)
-
     # CHT compressible
     cht_compressible           = TestCase('cht_compressible')
     cht_compressible.cfg_dir   = "coupled_cht/comp_2d"
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 6a715ac54781..040d9fba0207 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -1416,19 +1416,6 @@ def main():
     cht_incompressible.tol       = 0.00001
     test_list.append(cht_incompressible)
 
-    # CHT incompressible unsteady
-    cht_incompressible_unsteady           = TestCase('cht_incompressible_unsteady')
-    cht_incompressible_unsteady.cfg_dir   = "coupled_cht/incomp_2d_unsteady"
-    cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
-    cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-1.303588, -0.080377, -0.080380, -0.080377] #last 4 columns
-    cht_incompressible_unsteady.su2_exec  = "SU2_CFD"
-    cht_incompressible_unsteady.timeout   = 1600
-    cht_incompressible_unsteady.multizone = True
-    cht_incompressible_unsteady.unsteady  = True
-    cht_incompressible_unsteady.tol       = 0.00001
-    test_list.append(cht_incompressible_unsteady)
-
      # CHT compressible
     cht_incompressible           = TestCase('cht_compressible')
     cht_incompressible.cfg_dir   = "coupled_cht/comp_2d"
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
index 8ae969d93c79..01bf850b4e2d 100644
--- a/TestCases/tutorials.py
+++ b/TestCases/tutorials.py
@@ -42,6 +42,20 @@ def main():
     ### RUN TUTORIAL CASES             ###
     ######################################
 
+
+    # CHT incompressible unsteady
+    cht_incompressible_unsteady           = TestCase('cht_incompressible_unsteady')
+    cht_incompressible_unsteady.cfg_dir   = "../Tutorials/multiphysics/unsteady_cht/"
+    cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
+    cht_incompressible_unsteady.test_iter = 2
+    cht_incompressible_unsteady.test_vals = [-2.659390, -2.533160, -0.080399, -0.080399, -0.080399, -12.421450,  0.0000e+00, 0.0000e+00, 2.3824e+02] #last 9 columns
+    cht_incompressible_unsteady.su2_exec  = "mpirun -n 2 SU2_CFD"
+    cht_incompressible_unsteady.timeout   = 1600
+    cht_incompressible_unsteady.multizone = True
+    cht_incompressible_unsteady.unsteady  = True
+    cht_incompressible_unsteady.tol       = 0.00001
+    test_list.append(cht_incompressible_unsteady)
+
     ### Incompressible Flow
 
     # 2D pin case massflow periodic with heatflux BC and prescribed extracted outlet heat