diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 3719219bfee8..d3e6e3f1aa60 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -570,7 +570,7 @@ class CConfig {
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
-  su2double Relaxation_Factor_AdjFlow;  /*!< \brief Relaxation coefficient of the linear solver adjoint mean flow. */
+  su2double Relaxation_Factor_Adjoint;  /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
   su2double Relaxation_Factor_CHT;      /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
   su2double AdjTurb_Linear_Error;       /*!< \brief Min error of the turbulent adjoint linear solver for the implicit formulation. */
   su2double EntropyFix_Coeff;           /*!< \brief Entropy fix coefficient. */
@@ -3998,10 +3998,9 @@ class CConfig {
   su2double GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; }
 
   /*!
-   * \brief Get the relaxation coefficient of the linear solver for the implicit formulation.
-   * \return relaxation coefficient of the linear solver for the implicit formulation.
+   * \brief Get the relaxation factor for solution updates of adjoint solvers.
    */
-  su2double GetRelaxation_Factor_AdjFlow(void) const { return Relaxation_Factor_AdjFlow; }
+  su2double GetRelaxation_Factor_Adjoint(void) const { return Relaxation_Factor_Adjoint; }
 
   /*!
    * \brief Get the relaxation coefficient of the CHT coupling.
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index ec8979084c72..b1d0e74c5fbd 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1670,8 +1670,8 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0);
   /* DESCRIPTION: Custom number of threads used for additive domain decomposition for ILU and LU_SGS (0 is "auto"). */
   addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0);
-  /* DESCRIPTION: Relaxation of the flow equations solver for the implicit formulation */
-  addDoubleOption("RELAXATION_FACTOR_ADJFLOW", Relaxation_Factor_AdjFlow, 1.0);
+  /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */
+  addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
   /* DESCRIPTION: Relaxation of the CHT coupling */
   addDoubleOption("RELAXATION_FACTOR_CHT", Relaxation_Factor_CHT, 1.0);
   /* DESCRIPTION: Roe coefficient */
@@ -2832,20 +2832,9 @@ void CConfig::SetConfig_Parsing(char case_filename[MAX_STRING_SIZE]) {
           newString.append(": invalid option name");
           newString.append(". Check current SU2 options in config_template.cfg.");
           newString.append("\n");
-          if (!option_name.compare("EXT_ITER")) newString.append("Option EXT_ITER is deprecated as of v7.0. Please use TIME_ITER, OUTER_ITER or ITER \n"
-                                                                 "to specify the number of time iterations, outer multizone iterations or iterations, respectively.");
-          if (!option_name.compare("UNST_TIMESTEP")) newString.append("UNST_TIMESTEP is now TIME_STEP.\n");
-          if (!option_name.compare("UNST_TIME")) newString.append("UNST_TIME is now MAX_TIME.\n");
-          if (!option_name.compare("UNST_INT_ITER")) newString.append("UNST_INT_ITER is now INNER_ITER.\n");
-          if (!option_name.compare("RESIDUAL_MINVAL")) newString.append("RESIDUAL_MINVAL is now CONV_RESIDUAL_MINVAL.\n");
-          if (!option_name.compare("STARTCONV_ITER")) newString.append("STARTCONV_ITER is now CONV_STARTITER.\n");
-          if (!option_name.compare("CAUCHY_ELEMS")) newString.append("CAUCHY_ELEMS is now CONV_CAUCHY_ELEMS.\n");
-          if (!option_name.compare("CAUCHY_EPS")) newString.append("CAUCHY_EPS is now CONV_CAUCHY_EPS.\n");
-          if (!option_name.compare("OUTPUT_FORMAT")) newString.append("OUTPUT_FORMAT is now TABULAR_FORMAT.\n");
-          if (!option_name.compare("PHYSICAL_PROBLEM")) newString.append("PHYSICAL_PROBLEM is now SOLVER.\n");
-          if (!option_name.compare("REGIME_TYPE")) newString.append("REGIME_TYPE has been removed.\n "
-                                                                    "If you want use the incompressible solver, \n"
-                                                                    "use INC_EULER, INC_NAVIER_STOKES or INC_RANS as value of the SOLVER option.");
+          if (!option_name.compare("RELAXATION_FACTOR_ADJFLOW"))
+            newString.append("Option RELAXATION_FACTOR_ADJFLOW is now RELAXATION_FACTOR_ADJOINT, "
+                             "and it also applies to discrete adjoint problems\n.");
           errorString.append(newString);
           err_count++;
         continue;
diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
index 1041ad8941a7..56d775016af4 100644
--- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
+++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
@@ -180,6 +180,12 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
    */
   void Add_External_To_Solution(unsigned short iZone);
 
+  /*!
+   * \brief Puts Solution into SolutionOld.
+   * \param[in] iZone - Zone index.
+   */
+  void Set_SolutionOld_To_Solution(unsigned short iZone);
+
   /*!
    * \brief Extract contribution of iZone to jZone with BGS relaxation.
    * \param[in] iZone - Source zone (the one that was initialized).
diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp
index 2935d2c96e92..2528c88f3bf3 100644
--- a/SU2_CFD/include/solvers/CFEASolver.hpp
+++ b/SU2_CFD/include/solvers/CFEASolver.hpp
@@ -55,8 +55,8 @@ class CFEASolver : public CSolver {
   su2double Conv_Check[3];          /*!< \brief Current values for convergence check: UTOL, RTOL, ETOL. */
   su2double FSI_Conv[2];            /*!< \brief Values to check the convergence of the FSI problem. */
 
-  unsigned long idxIncrement;       /*!< \brief Index of the current load increment */
-  su2double loadIncrement;          /*!< \brief Coefficient that determines the amount of load which is applied. */
+  unsigned long idxIncrement = 0;   /*!< \brief Index of the current load increment */
+  su2double loadIncrement = 1.0;    /*!< \brief Coefficient that determines the amount of load which is applied. */
 
   su2double WAitken_Dyn;            /*!< \brief Aitken's dynamic coefficient. */
   su2double WAitken_Dyn_tn1;        /*!< \brief Aitken's dynamic coefficient in the previous iteration. */
diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
index 1403f674602e..941b6852fc89 100644
--- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
@@ -256,7 +256,7 @@ void CDiscAdjMultizoneDriver::Run() {
           Add_External_To_Solution(iZone);
         }
         else {
-          /*--- If we restarted, Solution already has all contribution,
+          /*--- If we restarted, Solution already has all contributions,
            *    we run only one inner iter to compute the cross terms. ---*/
           eval_transfer = true;
         }
@@ -274,7 +274,12 @@ void CDiscAdjMultizoneDriver::Run() {
                                                     solver_container, numerics_container, config_container,
                                                     surface_movement, grid_movement, FFDBox, iZone, INST_0);
 
-        /*--- Print out the convergence data to screen and history file ---*/
+        /*--- This is done explicitly here for multizone cases, only in inner iterations and not when
+         *    extracting cross terms so that the adjoint residuals in each zone still make sense. ---*/
+
+        Set_SolutionOld_To_Solution(iZone);
+
+        /*--- Print out the convergence data to screen and history file. ---*/
 
         bool converged = iteration_container[iZone][INST_0]->Monitor(output_container[iZone], integration_container,
                                                     geometry_container, solver_container, numerics_container,
@@ -293,6 +298,7 @@ void CDiscAdjMultizoneDriver::Run() {
           /*--- Extracting adjoints for solvers in jZone w.r.t. to the output of all solvers in iZone,
            *    that is, for the cases iZone != jZone we are evaluating cross derivatives between zones. ---*/
 
+          config_container[jZone]->SetInnerIter(0);
           iteration_container[jZone][INST_0]->Iterate(output_container[jZone], integration_container, geometry_container,
                                                       solver_container, numerics_container, config_container,
                                                       surface_movement, grid_movement, FFDBox, jZone, INST_0);
@@ -898,6 +904,15 @@ void CDiscAdjMultizoneDriver::Add_External_To_Solution(unsigned short iZone) {
   }
 }
 
+void CDiscAdjMultizoneDriver::Set_SolutionOld_To_Solution(unsigned short iZone) {
+
+  for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {
+    auto solver = solver_container[iZone][INST_0][MESH_0][iSol];
+    if (solver && solver->GetAdjoint())
+      solver->GetNodes()->Set_OldSolution();
+  }
+}
+
 void CDiscAdjMultizoneDriver::Update_Cross_Term(unsigned short iZone, unsigned short jZone) {
 
   for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {
diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp
index c6e8fede6aa4..2cfd3fedba80 100644
--- a/SU2_CFD/src/iteration_structure.cpp
+++ b/SU2_CFD/src/iteration_structure.cpp
@@ -473,11 +473,12 @@ void CFluidIteration::Iterate(COutput *output,
                                                                      RUNTIME_RADIATION_SYS, val_iZone, val_iInst);
   }
 
-  /*--- Adapt the CFL number using an exponential progression
-   with under-relaxation approach. ---*/
+  /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/
 
   if (config[val_iZone]->GetCFL_Adapt() == YES) {
-    solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone]);
+    SU2_OMP_PARALLEL
+    solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst],
+                                                                   solver[val_iZone][val_iInst], config[val_iZone]);
   }
 
   /*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/
diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
index 2c0e8e523317..8e4cbcd5ceb7 100644
--- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp
@@ -2301,7 +2301,7 @@ void CAdjEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **sol
 
   for (iPoint = 0; iPoint < nPointDomain; iPoint++)
     for (iVar = 0; iVar < nVar; iVar++) {
-      nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_AdjFlow()*LinSysSol[iPoint*nVar+iVar]);
+      nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_Adjoint()*LinSysSol[iPoint*nVar+iVar]);
     }
 
   /*--- MPI solution ---*/
diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp
index 6bf7e05d468a..2cadc23514a3 100644
--- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp
@@ -667,8 +667,6 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co
     }
   }
 
-  if(multizone) nodes->Set_OldSolution();
-
   SetResidual_RMS(geometry, config);
 
 }
diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
index 6523ea691858..bff05069cc2e 100644
--- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
+++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
@@ -476,26 +476,26 @@ void CDiscAdjSolver::SetAdj_ObjFunc(CGeometry *geometry, CConfig *config) {
 
 void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config){
 
-  bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
-  bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
-  bool multizone = config->GetMultizone_Problem();
+  const bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
+  const bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
+  const bool multizone = config->GetMultizone_Problem();
 
-  unsigned short iVar;
-  unsigned long iPoint;
-  su2double residual;
+  const su2double relax = (config->GetInnerIter()==0)? 1.0 : config->GetRelaxation_Factor_Adjoint();
 
   /*--- Set Residuals to zero ---*/
 
-  for (iVar = 0; iVar < nVar; iVar++) {
-      SetRes_RMS(iVar,0.0);
-      SetRes_Max(iVar,0.0,0);
+  for (auto iVar = 0u; iVar < nVar; iVar++) {
+    SetRes_RMS(iVar,0.0);
+    SetRes_Max(iVar,0.0,0);
   }
 
-  /*--- Set the old solution ---*/
+  /*--- Set the old solution and compute residuals. ---*/
 
   if(!multizone) nodes->Set_OldSolution();
 
-  for (iPoint = 0; iPoint < nPoint; iPoint++) {
+  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
+
+    const su2double isdomain = (iPoint < nPointDomain)? 1.0 : 0.0;
 
     /*--- Extract the adjoint solution ---*/
 
@@ -506,13 +506,22 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
       direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);
     }
 
-    /*--- Store the adjoint solution ---*/
+    /*--- Relax and store the adjoint solution, compute the residuals. ---*/
 
-    nodes->SetSolution(iPoint,Solution);
+    for (auto iVar = 0u; iVar < nVar; iVar++) {
+      su2double residual = relax*(Solution[iVar]-nodes->GetSolution_Old(iPoint,iVar));
+      nodes->AddSolution(iPoint, iVar, residual);
+
+      residual *= isdomain;
+      AddRes_RMS(iVar,pow(residual,2));
+      AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
+    }
   }
 
+  SetResidual_RMS(geometry, config);
+
   if (time_n_needed) {
-    for (iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
 
       /*--- Extract the adjoint solution at time n ---*/
 
@@ -523,8 +532,9 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
       nodes->Set_Solution_time_n(iPoint,Solution);
     }
   }
+
   if (time_n1_needed) {
-    for (iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
 
       /*--- Extract the adjoint solution at time n-1 ---*/
 
@@ -536,21 +546,6 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
     }
   }
 
-  /*--- Set the residuals ---*/
-
-  for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
-    for (iVar = 0; iVar < nVar; iVar++) {
-      residual = nodes->GetSolution(iPoint,iVar) - nodes->GetSolution_Old(iPoint,iVar);
-
-      AddRes_RMS(iVar,residual*residual);
-      AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
-    }
-  }
-
-  if(multizone) nodes->Set_OldSolution();
-
-  SetResidual_RMS(geometry, config);
-
 }
 
 void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) {
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index 1c1864f5f667..962a748acb00 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -48,8 +48,6 @@ CFEASolver::CFEASolver(bool mesh_deform_mode) : CSolver(mesh_deform_mode) {
   Total_CFEA = 0.0;
   WAitken_Dyn = 0.0;
   WAitken_Dyn_tn1 = 0.0;
-  idxIncrement = 0;
-  loadIncrement = 1.0;
 
   element_container = new CElement** [MAX_TERMS]();
   for (unsigned short iTerm = 0; iTerm < MAX_TERMS; iTerm++)
@@ -124,8 +122,6 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() {
   Total_CFEA        = 0.0;
   WAitken_Dyn       = 0.0;
   WAitken_Dyn_tn1   = 0.0;
-  idxIncrement      = 0;
-  loadIncrement     = 0.0;
 
   SetFSI_ConvValue(0,0.0);
   SetFSI_ConvValue(1,0.0);
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index f6c1ee273446..0b8e47a4c90a 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -2320,8 +2320,6 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
                              CSolver   ***solver_container,
                              CConfig   *config) {
 
-  /// TODO: Add OpenMP stuff to this method.
-
   /* Adapt the CFL number on all multigrid levels using an
    exponential progression with under-relaxation approach. */
 
@@ -2330,6 +2328,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
   const su2double CFLFactorIncrease = config->GetCFL_AdaptParam(1);
   const su2double CFLMin            = config->GetCFL_AdaptParam(2);
   const su2double CFLMax            = config->GetCFL_AdaptParam(3);
+  const bool fullComms              = (config->GetComm_Level() == COMM_FULL);
 
   for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
 
@@ -2365,6 +2364,9 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
     /* Check that we are meeting our nonlinear residual reduction target
      over time so that we do not get stuck in limit cycles. */
 
+    SU2_OMP_MASTER
+    { /* Only the master thread updates the shared variables. */
+
     Old_Func = New_Func;
     unsigned short Res_Count = 100;
     if (NonLinRes_Series.size() == 0) NonLinRes_Series.resize(Res_Count,0.0);
@@ -2408,17 +2410,30 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
      Reset the array so that we delay the next decrease for some iterations. */
 
     if (fabs(NonLinRes_Value) < 0.1*New_Func) {
-      reduceCFL = true;
       NonLinRes_Counter = 0;
       for (unsigned short iCounter = 0; iCounter < Res_Count; iCounter++)
         NonLinRes_Series[iCounter] = New_Func;
     }
 
+    } /* End SU2_OMP_MASTER, now all threads update the CFL number. */
+    SU2_OMP_BARRIER
+
+    if (fabs(NonLinRes_Value) < 0.1*New_Func) {
+      reduceCFL = true;
+    }
+
     /* Loop over all points on this grid and apply CFL adaption. */
 
-    su2double myCFLMin = 1e30;
-    su2double myCFLMax = 0.0;
-    su2double myCFLSum = 0.0;
+    su2double myCFLMin = 1e30, myCFLMax = 0.0, myCFLSum = 0.0;
+
+    SU2_OMP_MASTER
+    if ((iMesh == MESH_0) && fullComms) {
+      Min_CFL_Local = 1e30;
+      Max_CFL_Local = 0.0;
+      Avg_CFL_Local = 0.0;
+    }
+
+    SU2_OMP_FOR_STAT(roundUpDiv(geometry[iMesh]->GetnPointDomain(),omp_get_max_threads()))
     for (unsigned long iPoint = 0; iPoint < geometry[iMesh]->GetnPointDomain(); iPoint++) {
 
       /* Get the current local flow CFL number at this point. */
@@ -2476,20 +2491,37 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
         solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL);
       }
 
-      /* Store min and max CFL for reporting on fine grid. */
+      /* Store min and max CFL for reporting on the fine grid. */
 
-      myCFLMin = min(CFL,myCFLMin);
-      myCFLMax = max(CFL,myCFLMax);
-      myCFLSum += CFL;
+      if ((iMesh == MESH_0) && fullComms) {
+        myCFLMin = min(CFL,myCFLMin);
+        myCFLMax = max(CFL,myCFLMax);
+        myCFLSum += CFL;
+      }
 
     }
 
     /* Reduce the min/max/avg local CFL numbers. */
 
-    SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
-    SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
-    SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
+    if ((iMesh == MESH_0) && fullComms) {
+      SU2_OMP_CRITICAL
+      { /* OpenMP reduction. */
+        Min_CFL_Local = min(Min_CFL_Local,myCFLMin);
+        Max_CFL_Local = max(Max_CFL_Local,myCFLMax);
+        Avg_CFL_Local += myCFLSum;
+      }
+      SU2_OMP_BARRIER
+
+      SU2_OMP_MASTER
+      { /* MPI reduction. */
+        myCFLMin = Min_CFL_Local; myCFLMax = Max_CFL_Local; myCFLSum = Avg_CFL_Local;
+        SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+        SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+        SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
+      }
+      SU2_OMP_BARRIER
+    }
 
   }
 
diff --git a/TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg b/TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg
index 5f163c014380..b2db44316a48 100644
--- a/TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg
+++ b/TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg
@@ -174,7 +174,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/cont_adj_euler/wedge/inv_wedge_ROE_multiobj.cfg b/TestCases/cont_adj_euler/wedge/inv_wedge_ROE_multiobj.cfg
index 3df8aaf6ed09..48a5ed2b9f11 100644
--- a/TestCases/cont_adj_euler/wedge/inv_wedge_ROE_multiobj.cfg
+++ b/TestCases/cont_adj_euler/wedge/inv_wedge_ROE_multiobj.cfg
@@ -182,7 +182,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/deformation/naca4412/def_NACA4412.cfg b/TestCases/deformation/naca4412/def_NACA4412.cfg
index 858c07a59d6f..2aa7b22dbb7b 100644
--- a/TestCases/deformation/naca4412/def_NACA4412.cfg
+++ b/TestCases/deformation/naca4412/def_NACA4412.cfg
@@ -208,7 +208,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.15
diff --git a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg
index 38c0323a7cf8..3ad4de8b5ecf 100644
--- a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg
+++ b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg
@@ -351,7 +351,7 @@ CONV_NUM_METHOD_ADJFLOW= JST
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.96
diff --git a/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg b/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg
index feb7314f4015..f2cd3ccfa097 100644
--- a/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg
+++ b/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg
@@ -240,7 +240,7 @@ CONV_NUM_METHOD_ADJFLOW= JST
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 0.95
+RELAXATION_FACTOR_ADJOINT= 0.95
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/disc_adj_fsi/Airfoil_2d/configFEA.cfg b/TestCases/disc_adj_fsi/Airfoil_2d/configFEA.cfg
index 65d6aa959100..caed22ee0683 100755
--- a/TestCases/disc_adj_fsi/Airfoil_2d/configFEA.cfg
+++ b/TestCases/disc_adj_fsi/Airfoil_2d/configFEA.cfg
@@ -31,11 +31,11 @@ DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT
 DISCADJ_LIN_PREC= ILU
 %
 % Convergence criteria ------------------------------------------------- %
-% interaction
 BGS_RELAXATION= FIXED_PARAMETER
 STAT_RELAX_PARAMETER= 1.0
-% solid
 INNER_ITER= 11
+% Incremental load on to test that it does not break the discadj, it should only affect the direct solution.
+INCREMENTAL_LOAD= YES
 %
 % In\Out --------------------------------------------------------------- %
 MESH_FILENAME= mesh.su2
diff --git a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_2surf_1obj.cfg b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_2surf_1obj.cfg
index 110bfbca15f8..0b75e05ca9e5 100644
--- a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_2surf_1obj.cfg
+++ b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_2surf_1obj.cfg
@@ -181,7 +181,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj.cfg b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj.cfg
index 0b41b93e7dcc..fc27543f6fea 100644
--- a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj.cfg
+++ b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj.cfg
@@ -182,7 +182,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_1surf.cfg b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_1surf.cfg
index d555757702fa..d4481ff19844 100644
--- a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_1surf.cfg
+++ b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_1surf.cfg
@@ -182,7 +182,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_combo.cfg b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_combo.cfg
index c7f6707002ae..ee65e8e106b0 100644
--- a/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_combo.cfg
+++ b/TestCases/optimization_euler/multiobjective_wedge/inv_wedge_ROE_multiobj_combo.cfg
@@ -184,7 +184,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 2da7dd7763a4..4f51709df64f 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -58,7 +58,7 @@ def main():
     discadj_cylinder3D.cfg_dir   = "disc_adj_euler/cylinder3D"
     discadj_cylinder3D.cfg_file  = "inv_cylinder3D.cfg"
     discadj_cylinder3D.test_iter = 5
-    discadj_cylinder3D.test_vals = [-3.724803, -3.838647, 0.000000, 0.000000] #last 4 columns
+    discadj_cylinder3D.test_vals = [-3.755976, -3.858734, 0.000000, 0.000000] #last 4 columns
     discadj_cylinder3D.su2_exec  = "parallel_computation.py -f"
     discadj_cylinder3D.timeout   = 1600
     discadj_cylinder3D.tol       = 0.00001
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 9dd4a5102350..89f8e893b716 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -58,7 +58,7 @@ def main():
     discadj_cylinder3D.cfg_dir   = "disc_adj_euler/cylinder3D"
     discadj_cylinder3D.cfg_file  = "inv_cylinder3D.cfg"
     discadj_cylinder3D.test_iter = 5
-    discadj_cylinder3D.test_vals = [-3.728728, -3.844274, -0.000000, 0.000000] #last 4 columns
+    discadj_cylinder3D.test_vals = [-3.759637, -3.864023, -0.000000, 0.000000] #last 4 columns
     discadj_cylinder3D.su2_exec  = "SU2_CFD_AD"
     discadj_cylinder3D.timeout   = 1600
     discadj_cylinder3D.tol       = 0.00001
diff --git a/config_template.cfg b/config_template.cfg
index 23d866db7767..0ce363c8b638 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1099,8 +1099,8 @@ CONV_NUM_METHOD_ADJFLOW= JST
 % Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT)
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
-% Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+% Relaxation coefficient (also for discrete adjoint problems)
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8