diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp
index c3ba75d92042..b12d0068c50c 100644
--- a/Common/include/basic_types/ad_structure.hpp
+++ b/Common/include/basic_types/ad_structure.hpp
@@ -74,6 +74,11 @@ inline void RegisterInput(su2double& data, bool push_index = true) {}
  */
 inline void RegisterOutput(su2double& data) {}
 
+/*!
+ * \brief Resize the adjoint vector, for subsequent access without bounds checking.
+ */
+inline void ResizeAdjoints() {}
+
 /*!
  * \brief Sets the adjoint value at index to val
  * \param[in] index - Position in the adjoint vector.
@@ -369,11 +374,24 @@ FORCEINLINE void Reset() {
   }
 }
 
+FORCEINLINE void ResizeAdjoints() { AD::getTape().resizeAdjointVector(); }
+
 FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); }
 
-FORCEINLINE void SetDerivative(int index, const double val) { AD::getTape().setGradient(index, val); }
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// When using it, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE void SetDerivative(int index, const double val) {
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  AD::getTape().setGradient(index, val, BoundsChecking::False);
+}
 
-FORCEINLINE double GetDerivative(int index) { return AD::getTape().getGradient(index); }
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// If called after tape evaluations, the adjoints should exist.
+// Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE double GetDerivative(int index) {
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  return AD::getTape().getGradient(index, BoundsChecking::False);
+}
 
 FORCEINLINE bool IsIdentifierActive(su2double const& value) {
   return getTape().isIdentifierActive(value.getIdentifier());
@@ -523,26 +541,14 @@ FORCEINLINE void delete_handler(void* handler) {
 
 FORCEINLINE bool BeginPassive() {
   if (AD::getTape().isActive()) {
-    StopRecording();
+    AD::getTape().setPassive();
     return true;
   }
   return false;
 }
 
 FORCEINLINE void EndPassive(bool wasActive) {
-  if (wasActive) StartRecording();
-}
-
-FORCEINLINE bool PausePreaccumulation() {
-  const auto current = PreaccEnabled;
-  if (!current) return false;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
-  return true;
-}
-
-FORCEINLINE void ResumePreaccumulation(bool wasActive) {
-  if (!wasActive) return;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+  if (wasActive) AD::getTape().setActive();
 }
 
 FORCEINLINE void StartNoSharedReading() {
@@ -558,6 +564,19 @@ FORCEINLINE void EndNoSharedReading() {
   opdi::logic->addReverseBarrier();
 #endif
 }
+
+FORCEINLINE bool PausePreaccumulation() {
+  const auto current = PreaccEnabled;
+  if (!current) return false;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
+  return true;
+}
+
+FORCEINLINE void ResumePreaccumulation(bool wasActive) {
+  if (!wasActive) return;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+}
+
 #endif  // CODI_REVERSE_TYPE
 
 void Initialize();
diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp
index bbc5eebb76e9..958a82f0e389 100644
--- a/Common/include/code_config.hpp
+++ b/Common/include/code_config.hpp
@@ -96,7 +96,7 @@ FORCEINLINE Out su2staticcast_p(In ptr) {
 #include "codi/tools/data/externalFunctionUserData.hpp"
 
 #if defined(HAVE_OMP)
-using su2double = codi::RealReverseIndexOpenMP;
+using su2double = codi::RealReverseIndexOpenMPGen<double, double>;
 #else
 #if defined(CODI_INDEX_TAPE)
 using su2double = codi::RealReverseIndex;
diff --git a/Common/include/parallelization/omp_structure.hpp b/Common/include/parallelization/omp_structure.hpp
index de19ec1e544d..60baa3e44dda 100644
--- a/Common/include/parallelization/omp_structure.hpp
+++ b/Common/include/parallelization/omp_structure.hpp
@@ -188,12 +188,14 @@ void omp_finalize();
  * thread, with all threads and memory views synchronized both beforehand and afterwards.
  */
 
-#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  SU2_OMP_BARRIER                        \
+#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  SU2_OMP_BARRIER                                    \
+  if (omp_in_parallel()) AD::StartNoSharedReading(); \
   SU2_OMP_MASTER
 
-#define END_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  END_SU2_OMP_MASTER                   \
+#define END_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  END_SU2_OMP_MASTER                               \
+  if (omp_in_parallel()) AD::EndNoSharedReading(); \
   SU2_OMP_BARRIER
 
 #define SU2_OMP_SAFE_GLOBAL_ACCESS(...) BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS{__VA_ARGS__} END_SU2_OMP_SAFE_GLOBAL_ACCESS
diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
index fbaf2fd1def0..f561618b6ad8 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -192,7 +192,7 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
                                                   const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                   unsigned long m, ScalarType& residual, bool monitoring,
                                                   const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
 
@@ -241,16 +241,22 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master && (lin_sol_mode != LINEAR_SOLVER_MODE::MESH_DEFORM))
+      if (masterRank && (lin_sol_mode != LINEAR_SOLVER_MODE::MESH_DEFORM)) {
+        SU2_OMP_MASTER
         cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual ---*/
 
-    if (monitoring && master) {
-      WriteHeader("CG", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if (monitoring && masterRank) {
+      SU2_OMP_MASTER {
+        WriteHeader("CG", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -281,7 +287,11 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
 
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
 
     precond(r, z);
@@ -300,7 +310,11 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("CG", i, norm_r / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("CG", i, norm_r / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, A_x);
@@ -308,8 +322,10 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
       ScalarType true_res = r.norm();
 
       if (fabs(true_res - norm_r) > tol * 10.0) {
-        if (master) {
+        if (masterRank) {
+          SU2_OMP_MASTER
           WriteWarning(norm_r, true_res, tol);
+          END_SU2_OMP_MASTER
         }
       }
     }
@@ -325,7 +341,7 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
                                                       const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                       unsigned long m, ScalarType& residual, bool monitoring,
                                                       const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   const bool flexible = !precond.IsIdentity();
 
   /*---  Check the subspace size ---*/
@@ -386,7 +402,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   if ((beta < tol * norm0) || (beta < eps)) {
     /*--- System is already solved ---*/
 
-    if (master) cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl;
+    if (masterRank) {
+      SU2_OMP_MASTER
+      cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl;
+      END_SU2_OMP_MASTER
+    }
     residual = beta;
     return 0;
   }
@@ -403,9 +423,12 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   /*--- Output header information including initial residual ---*/
 
   unsigned long i = 0;
-  if ((monitoring) && (master)) {
-    WriteHeader("FGMRES", tol, beta);
-    WriteHistory(i, beta / norm0);
+  if ((monitoring) && (masterRank)) {
+    SU2_OMP_MASTER {
+      WriteHeader("FGMRES", tol, beta);
+      WriteHistory(i, beta / norm0);
+    }
+    END_SU2_OMP_MASTER
   }
 
   /*---  Loop over all search directions ---*/
@@ -444,7 +467,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
 
     /*---  Output the relative residual if necessary ---*/
 
-    if ((((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) && (master)) WriteHistory(i + 1, beta / norm0);
+    if ((((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0))) {
+      SU2_OMP_MASTER
+      WriteHistory(i + 1, beta / norm0);
+      END_SU2_OMP_MASTER
+    }
   }
 
   /*---  Solve the least-squares system and update solution ---*/
@@ -460,7 +487,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   /*---  Recalculate final (neg.) residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("FGMRES", i, beta / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("FGMRES", i, beta / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, W[0]);
@@ -468,8 +499,10 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
       ScalarType res = W[0].norm();
 
       if (fabs(res - beta) > tol * 10) {
-        if (master) {
+        if (masterRank) {
+          SU2_OMP_MASTER
           WriteWarning(beta, res, tol);
+          END_SU2_OMP_MASTER
         }
       }
     }
@@ -511,7 +544,7 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
                                                        const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                        unsigned long m, ScalarType& residual, bool monitoring,
                                                        const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
 
@@ -561,15 +594,22 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master) cout << "CSysSolve::BCGSTAB(): system solved by initial guess." << endl;
+      if (masterRank) {
+        SU2_OMP_MASTER
+        cout << "CSysSolve::BCGSTAB(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual ---*/
 
-    if ((monitoring) && (master)) {
-      WriteHeader("BCGSTAB", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if ((monitoring) && (masterRank)) {
+      SU2_OMP_MASTER {
+        WriteHeader("BCGSTAB", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -637,22 +677,32 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
 
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0) && (master)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
   }
 
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("BCGSTAB", i, norm_r / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("BCGSTAB", i, norm_r / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, A_x);
       r = b - A_x;
       ScalarType true_res = r.norm();
 
-      if ((fabs(true_res - norm_r) > tol * 10.0) && (master)) {
+      if ((fabs(true_res - norm_r) > tol * 10.0) && (masterRank)) {
+        SU2_OMP_MASTER
         WriteWarning(norm_r, true_res, tol);
+        END_SU2_OMP_MASTER
       }
     }
   }
@@ -667,7 +717,7 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
                                                         const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                         unsigned long m, ScalarType& residual, bool monitoring,
                                                         const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   const bool fix_iter_mode = tol < eps;
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
@@ -717,15 +767,22 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master) cout << "CSysSolve::Smoother_LinSolver(): system solved by initial guess." << endl;
+      if (masterRank) {
+        SU2_OMP_MASTER
+        cout << "CSysSolve::Smoother_LinSolver(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual. ---*/
 
-    if ((monitoring) && (master)) {
-      WriteHeader("Smoother", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if ((monitoring) && (masterRank)) {
+      SU2_OMP_MASTER {
+        WriteHeader("Smoother", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -759,14 +816,20 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
     if (!fix_iter_mode && config->GetComm_Level() == COMM_FULL) {
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
   }
 
   if (fix_iter_mode) norm_r = r.norm();
 
-  if ((monitoring) && (master) && (config->GetComm_Level() == COMM_FULL)) {
+  if ((monitoring) && (masterRank) && (config->GetComm_Level() == COMM_FULL)) {
+    SU2_OMP_MASTER
     WriteFinalResidual("Smoother", i, norm_r / norm0);
+    END_SU2_OMP_MASTER
   }
 
   residual = norm_r / norm0;
diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp
index a37971f0aeb8..0ecf23e17fc2 100644
--- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp
@@ -274,6 +274,7 @@ void CDiscAdjFEAIteration::InitializeAdjoint(CSolver***** solver, CGeometry****
                                              unsigned short iZone, unsigned short iInst) {
   /*--- Initialize the adjoints the conservative variables ---*/
 
+  AD::ResizeAdjoints();
   solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
 }
 
diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
index 5b4519078651..004d9cabfa47 100644
--- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
@@ -360,6 +360,8 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry***
   auto solvers0 = solver[iZone][iInst][MESH_0];
   auto geometry0 = geometry[iZone][iInst][MESH_0];
 
+  AD::ResizeAdjoints();
+
   SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Initialize the adjoints the conservative variables ---*/
diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
index a7935c92b19d..0d220db739e0 100644
--- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
@@ -157,6 +157,8 @@ void CDiscAdjHeatIteration::InitializeAdjoint(CSolver***** solver, CGeometry****
                                               unsigned short iZone, unsigned short iInst) {
 
   /*--- Initialize the adjoints the solution variables ---*/
+
+  AD::ResizeAdjoints();
   solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
 }
 
diff --git a/externals/codi b/externals/codi
index 21e4f6e084d6..4237fe38c90b 160000
--- a/externals/codi
+++ b/externals/codi
@@ -1 +1 @@
-Subproject commit 21e4f6e084d635e7834ca55dc9d83b667a5aca2b
+Subproject commit 4237fe38c90b83541f699bd84e00c656099b7d18
diff --git a/externals/opdi b/externals/opdi
index f33b507c24f7..3ff7bb16bd70 160000
--- a/externals/opdi
+++ b/externals/opdi
@@ -1 +1 @@
-Subproject commit f33b507c24f7448d4cf4df16ab5c53ea254b8774
+Subproject commit 3ff7bb16bd70aa8f74277704b5feccaaaf6838bd
diff --git a/meson_scripts/init.py b/meson_scripts/init.py
index 47f5841a74a6..f1e298a9cbbf 100755
--- a/meson_scripts/init.py
+++ b/meson_scripts/init.py
@@ -54,11 +54,11 @@ def init_submodules(
 
     # This information of the modules is used if projects was not cloned using git
     # The sha tag must be maintained manually to point to the correct commit
-    sha_version_codi = "21e4f6e084d635e7834ca55dc9d83b667a5aca2b"
+    sha_version_codi = "4237fe38c90b83541f699bd84e00c656099b7d18"
     github_repo_codi = "https://github.com/scicompkl/CoDiPack"
     sha_version_medi = "aafc2d1966ba1233640af737e71c77c1a86183fd"
     github_repo_medi = "https://github.com/SciCompKL/MeDiPack"
-    sha_version_opdi = "f33b507c24f7448d4cf4df16ab5c53ea254b8774"
+    sha_version_opdi = "3ff7bb16bd70aa8f74277704b5feccaaaf6838bd"
     github_repo_opdi = "https://github.com/SciCompKL/OpDiLib"
     sha_version_meson = "41c650a040d50e0912d268af7a903a9ce1456dfa"
     github_repo_meson = "https://github.com/mesonbuild/meson"