diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml
index cef9d68216ed..73533f6ce341 100644
--- a/.github/workflows/regression.yml
+++ b/.github/workflows/regression.yml
@@ -23,7 +23,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP, ReverseOMP, ForwardOMP]
+        config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, ReverseTagNoMPI, BaseOMP, ReverseOMP, ForwardOMP]
         include:
           - config_set: BaseMPI
             flags: '-Denable-pywrapper=true -Denable-coolprop=true -Denable-mpp=true -Dinstall-mpp=true -Denable-mlpcpp=true -Denable-tests=true --warnlevel=2'
@@ -37,6 +37,8 @@ jobs:
             flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror'
           - config_set: ForwardNoMPI
             flags: '-Denable-directdiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-tests=true --warnlevel=3 --werror'
+          - config_set: ReverseTagNoMPI
+            flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror -Dcodi-tape=Tag'
           - config_set: BaseOMP
             flags: '-Dwith-omp=true -Denable-mixedprec=true -Denable-pywrapper=true -Denable-tecio=false --warnlevel=3 --werror'
           - config_set: ReverseOMP
@@ -215,6 +217,55 @@ jobs:
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
+  discadj_tape_tests:
+    if: inputs.runner != 'ARM64'
+    runs-on: ${{ inputs.runner || 'ubuntu-latest' }}
+    name: Tape Tests
+    needs: build
+    strategy:
+      fail-fast: false
+      matrix:
+        testscript: ['serial_regression_AD.py']
+        include:
+          - testscript: 'serial_regression_AD.py'
+            tag: TagNoMPI
+    steps:
+      - name: Pre Cleanup
+        uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
+        with:
+          entrypoint: /bin/rm
+          args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
+      - name: Download All artifacts
+        uses: actions/download-artifact@v4
+      - name: Uncompress and Move Binaries
+        run: |
+          BIN_FOLDER="$PWD/install/bin"
+          mkdir -p $BIN_FOLDER
+          ls -lah $BIN_FOLDER
+          for type in Base Reverse Forward; do
+            TYPE_FOLDER="${type}${{matrix.tag}}"
+            echo "Processing '$TYPE_FOLDER' ..."
+            if [ -d $TYPE_FOLDER ]; then
+              pushd $TYPE_FOLDER
+              ls -lah
+              tar -zxvf install_bin.tgz
+              ls -lah install/bin/
+              cp -r install/* $BIN_FOLDER/../
+              popd;
+            fi
+          done
+          chmod a+x $BIN_FOLDER/*
+          ls -lahR $BIN_FOLDER
+      - name: Run Tests in Container
+        uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
+        with:
+          # -t <Tutorials-branch> -c <Testcases-branch>
+          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+      - name: Cleanup
+        uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
+        with:
+          entrypoint: /bin/rm
+          args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
   thread_sanitizer_tests:
     if: inputs.runner != 'ARM64'
     runs-on: ${{ inputs.runner || 'ubuntu-latest' }}
diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 1d98b92534e7..084fc3e487c1 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -597,8 +597,10 @@ class CConfig {
   //bool ContactResistance = false; /*!< \brief Apply contact resistance for conjugate heat transfer. */
   unsigned short nID_DV;   /*!< \brief ID for the region of FEM when computed using direct differentiation. */
 
-  bool AD_Mode;             /*!< \brief Algorithmic Differentiation support. */
-  bool AD_Preaccumulation;  /*!< \brief Enable or disable preaccumulation in the AD mode. */
+  bool AD_Mode;                         /*!< \brief Algorithmic Differentiation support. */
+  bool AD_Preaccumulation;              /*!< \brief Enable or disable preaccumulation in the AD mode. */
+  CHECK_TAPE_TYPE AD_CheckTapeType;           /*!< \brief Type of tape that is checked in a tape debug run. */
+  CHECK_TAPE_VARIABLES AD_CheckTapeVariables; /*!< \brief Type of variables that are checked in a tape debug run. */
   STRUCT_COMPRESS Kind_Material_Compress;  /*!< \brief Determines if the material is compressible or incompressible (structural analysis). */
   STRUCT_MODEL Kind_Material;              /*!< \brief Determines the material model to be used (structural analysis). */
   STRUCT_DEFORMATION Kind_Struct_Solver;   /*!< \brief Determines the geometric condition (small or large deformations) for structural analysis. */
@@ -822,6 +824,7 @@ class CConfig {
   SurfSens_FileName,             /*!< \brief Output file for the sensitivity on the surface (discrete adjoint). */
   VolSens_FileName,              /*!< \brief Output file for the sensitivity in the volume (discrete adjoint). */
   ObjFunc_Hess_FileName;         /*!< \brief Hessian approximation obtained by the Sobolev smoothing solver. */
+  bool Multizone_Adapt_FileName; /*!< \brief Append zone number to solution and restart file names. */
 
   bool
   Wrt_Performance,           /*!< \brief Write the performance summary at the end of a calculation.  */
@@ -1055,7 +1058,8 @@ class CConfig {
   long ParMETIS_pointWgt;           /*!< \brief Load balancing weight given to points. */
   long ParMETIS_edgeWgt;            /*!< \brief Load balancing weight given to edges. */
   unsigned short DirectDiff;        /*!< \brief Direct Differentation mode. */
-  bool DiscreteAdjoint;                /*!< \brief AD-based discrete adjoint mode. */
+  bool DiscreteAdjoint,                /*!< \brief AD-based discrete adjoint mode. */
+  DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
   su2double Const_DES;                 /*!< \brief Detached Eddy Simulation Constant. */
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
@@ -1897,6 +1901,12 @@ class CConfig {
    */
   su2double GetPressure_FreeStreamND(void) const { return Pressure_FreeStreamND; }
 
+  /*!
+   * \brief Get a reference to the non-dimensionalized freestream pressure (used for AD tracking).
+   * \return Reference to non-dimensionalized freestream pressure.
+   */
+  su2double& GetPressure_FreeStreamND(void) { return Pressure_FreeStreamND; }
+
   /*!
    * \brief Get the value of the thermodynamic pressure.
    * \return Thermodynamic pressure.
@@ -1922,6 +1932,12 @@ class CConfig {
    */
   su2double GetTemperature_FreeStreamND(void) const { return Temperature_FreeStreamND; }
 
+  /*!
+   * \brief Get a reference to the non-dimensionalized freestream temperature (used for AD tracking).
+   * \return Reference to non-dimensionalized freestream temperature.
+   */
+  su2double& GetTemperature_FreeStreamND(void) { return Temperature_FreeStreamND; }
+
   /*!
    * \brief Get the value of the non-dimensionalized vibrational-electronic freestream temperature.
    * \return Non-dimensionalized vibrational-electronic freestream temperature.
@@ -5384,19 +5400,25 @@ class CConfig {
   bool GetWrt_Volume_Overwrite(void) const { return Wrt_Volume_Overwrite; }
 
   /*!
-   * \brief Provides the number of varaibles.
+   * \brief Get whether filenames are appended the zone number automatically (multiphysics solver).
+   * \return Flag for appending zone numbers to restart and solution filenames. If Flag=true, zone numer is appended.
+   */
+  bool GetMultizone_AdaptFilename(void) const { return Multizone_Adapt_FileName; }
+
+  /*!
+   * \brief Provides the number of variables.
    * \return Number of variables.
    */
   unsigned short GetnVar(void);
 
   /*!
-   * \brief Provides the number of varaibles.
+   * \brief Provides the number of variables.
    * \return Number of variables.
    */
   unsigned short GetnZone(void) const { return nZone; }
 
   /*!
-   * \brief Provides the number of varaibles.
+   * \brief Provides the number of variables.
    * \return Number of variables.
    */
   unsigned short GetiZone(void) const { return iZone; }
@@ -8779,6 +8801,12 @@ class CConfig {
    */
   bool GetDiscrete_Adjoint(void) const { return DiscreteAdjoint; }
 
+  /*!
+   * \brief Get the indicator whether a debug run for the discrete adjoint solver will be started.
+   * \return the discrete adjoint debug indicator.
+   */
+  bool GetDiscrete_Adjoint_Debug(void) const { return DiscreteAdjointDebug; }
+
   /*!
    * \brief Get the number of subiterations while a ramp is applied.
    * \return Number of FSI subiters.
@@ -9216,6 +9244,16 @@ class CConfig {
    */
   su2double GetConst_DES(void) const { return Const_DES; }
 
+  /*!
+   * \brief Get the type of tape that will be checked in a tape debug run.
+   */
+  CHECK_TAPE_TYPE GetAD_CheckTapeType(void) const { return AD_CheckTapeType; }
+
+  /*!
+   * \brief Get the type of variables that will be checked for in a tape debug run.
+   */
+  CHECK_TAPE_VARIABLES GetAD_CheckTapeVariables(void) const { return AD_CheckTapeVariables; }
+
   /*!
    * \brief Get if AD preaccumulation should be performed.
    */
@@ -9605,6 +9643,11 @@ class CConfig {
    */
   unsigned short GetnVolumeOutputFiles() const { return nVolumeOutputFiles; }
 
+  /*!
+   * \brief GetnVolumeOutputFrequencies
+   */
+  unsigned short GetnVolumeOutputFrequencies() const { return nVolumeOutputFrequencies; }
+
   /*!
    * \brief GetVolumeOutputFrequency
    * \param[in] iFile: index of file number for which the writing frequency needs to be returned.
diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp
index 0bd979b9663a..5a02b23c2c67 100644
--- a/Common/include/basic_types/ad_structure.hpp
+++ b/Common/include/basic_types/ad_structure.hpp
@@ -38,6 +38,9 @@
  */
 namespace AD {
 #ifndef CODI_REVERSE_TYPE
+
+using Identifier = int;
+
 /*!
  * \brief Start the recording of the operations and involved variables.
  * If called, the computational graph of all operations occuring after the call will be stored,
@@ -45,6 +48,17 @@ namespace AD {
  */
 inline void StartRecording() {}
 
+/*!
+ * \brief Pause the recording of the operations and involved variables.
+ * If called, all operations occuring after the call will not be stored on the computational graph.
+ */
+inline void PauseRecording() {}
+
+/*!
+ * \brief Resume the recording of the operations and variables after the recording had been paused.
+ */
+inline void ResumeRecording() {}
+
 /*!
  * \brief Stops the recording of the operations and variables.
  */
@@ -101,14 +115,20 @@ inline void EndUseAdjoints() {}
  * \param[in] index - Position in the adjoint vector.
  * \param[in] val - adjoint value to be set.
  */
-inline void SetDerivative(int index, const double val) {}
+inline void SetDerivative(Identifier index, const double val) {}
 
 /*!
  * \brief Extracts the adjoint value at index
  * \param[in] index - position in the adjoint vector where the derivative will be extracted.
  * \return Derivative value.
  */
-inline double GetDerivative(int index) { return 0.0; }
+inline double GetDerivative(Identifier index) { return 0.0; }
+
+/*!
+ * \brief Returns the identifier that represents an inactive variable.
+ * \return Passive index.
+ */
+inline Identifier GetPassiveIndex() { return 0; }
 
 /*!
  * \brief Clears the currently stored adjoints but keeps the computational graph.
@@ -259,7 +279,50 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {}
  * \param[in] data - variable whose gradient information will be extracted.
  * \param[in] index - where obtained gradient information will be stored.
  */
-inline void SetIndex(int& index, const su2double& data) {}
+inline void SetIndex(Identifier& index, const su2double& data) {}
+
+/*!
+ * \brief Sets the tag tape to a specific tag.
+ * \param[in] tag - the number to which the tag is set.
+ */
+inline void SetTag(int tag) {}
+
+/*!
+ * \brief Sets the tag of a variable to 0.
+ * \param[in] v - the variable whose tag is cleared.
+ */
+inline void ClearTagOnVariable(su2double& v) {}
+
+/*!
+ * \brief Struct to store information about errors during a tag debug run.
+ */
+struct ErrorReport {};
+
+/*!
+ * \brief Set a reference to the output file of an ErrorReport.
+ * \param[in] report - the ErrorReport whose output file is set.
+ * \param[in] output_file - pointer to the output file.
+ */
+inline void SetDebugReportFile(ErrorReport& report, std::ostream* output_file) {}
+
+/*!
+ * \brief Set the ErrorReport to which error information from a tag debug recording is written.
+ * \param[in] report - the ErrorReport to which error information is written.
+ */
+inline void SetTagErrorCallback(ErrorReport& report) {}
+
+/*!
+ * \brief Reset the error counter in an ErrorReport.
+ * \param[in] report - the ErrorReport whose error counter is resetted.
+ */
+inline void ResetErrorCounter(ErrorReport& report) {}
+
+/*!
+ * \brief Get the error count of an ErrorReport.
+ * \param[in] report - the ErrorReport whose pointer to its error counter is returned.
+ * \return Value of the error counter.
+ */
+inline unsigned long GetErrorCount(const ErrorReport& report) { return 0; }
 
 /*!
  * \brief Pushes back the current tape position to the tape position's vector.
@@ -304,6 +367,7 @@ inline void EndNoSharedReading() {}
 using CheckpointHandler = codi::ExternalFunctionUserData;
 
 using Tape = su2double::Tape;
+using Identifier = su2double::Identifier;
 
 #ifdef HAVE_OPDI
 using ExtFuncHelper = codi::OpenMPExternalFunctionHelper<su2double>;
@@ -349,6 +413,10 @@ FORCEINLINE void ResetInput(su2double& data) { data = data.getValue(); }
 
 FORCEINLINE void StartRecording() { AD::getTape().setActive(); }
 
+FORCEINLINE void PauseRecording() { AD::getTape().setPassive(); }
+
+FORCEINLINE void ResumeRecording() { AD::getTape().setActive(); }
+
 FORCEINLINE void StopRecording() { AD::getTape().setPassive(); }
 
 FORCEINLINE bool TapeActive() { return AD::getTape().isActive(); }
@@ -470,15 +538,15 @@ FORCEINLINE void BeginUseAdjoints() { AD::getTape().beginUseAdjointVector(); }
 
 FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); }
 
-FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); }
+FORCEINLINE void SetIndex(Identifier& index, const su2double& data) { index = data.getIdentifier(); }
 
 // 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().
 // This method does not perform locking either.
 // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints().
-FORCEINLINE void SetDerivative(int index, const double val) {
-  if (index == 0)  // Allow multiple threads to "set the derivative" of passive variables without causing data races.
-    return;
+FORCEINLINE void SetDerivative(Identifier index, const double val) {
+  // Allow multiple threads to "set the derivative" of passive variables without causing data races.
+  if (!AD::getTape().isIdentifierActive(index)) return;
 
   AD::getTape().setGradient(index, val, codi::AdjointsManagement::Manual);
 }
@@ -488,13 +556,11 @@ FORCEINLINE void SetDerivative(int index, const double val) {
 // Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
 // This method does not perform locking either.
 // It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints().
-FORCEINLINE double GetDerivative(int index) {
+FORCEINLINE double GetDerivative(Identifier index) {
   return AD::getTape().getGradient(index, codi::AdjointsManagement::Manual);
 }
 
-FORCEINLINE bool IsIdentifierActive(su2double const& value) {
-  return getTape().isIdentifierActive(value.getIdentifier());
-}
+FORCEINLINE Identifier GetPassiveIndex() { return AD::getTape().getPassiveIndex(); }
 
 /*--- Base case for parameter pack expansion. ---*/
 FORCEINLINE void SetPreaccIn() {}
@@ -502,7 +568,7 @@ FORCEINLINE void SetPreaccIn() {}
 template <class T, class... Ts, su2enable_if<std::is_same<T, su2double>::value> = 0>
 FORCEINLINE void SetPreaccIn(const T& data, Ts&&... moreData) {
   if (!PreaccActive) return;
-  if (IsIdentifierActive(data)) PreaccHelper.addInput(data);
+  PreaccHelper.addInput(data);
   SetPreaccIn(moreData...);
 }
 
@@ -515,9 +581,7 @@ template <class T>
 FORCEINLINE void SetPreaccIn(const T& data, const int size) {
   if (PreaccActive) {
     for (int i = 0; i < size; i++) {
-      if (IsIdentifierActive(data[i])) {
-        PreaccHelper.addInput(data[i]);
-      }
+      PreaccHelper.addInput(data[i]);
     }
   }
 }
@@ -527,9 +591,7 @@ FORCEINLINE void SetPreaccIn(const T& data, const int size_x, const int size_y)
   if (!PreaccActive) return;
   for (int i = 0; i < size_x; i++) {
     for (int j = 0; j < size_y; j++) {
-      if (IsIdentifierActive(data[i][j])) {
-        PreaccHelper.addInput(data[i][j]);
-      }
+      PreaccHelper.addInput(data[i][j]);
     }
   }
 }
@@ -547,7 +609,7 @@ FORCEINLINE void SetPreaccOut() {}
 template <class T, class... Ts, su2enable_if<std::is_same<T, su2double>::value> = 0>
 FORCEINLINE void SetPreaccOut(T& data, Ts&&... moreData) {
   if (!PreaccActive) return;
-  if (IsIdentifierActive(data)) PreaccHelper.addOutput(data);
+  PreaccHelper.addOutput(data);
   SetPreaccOut(moreData...);
 }
 
@@ -555,9 +617,7 @@ template <class T>
 FORCEINLINE void SetPreaccOut(T&& data, const int size) {
   if (PreaccActive) {
     for (int i = 0; i < size; i++) {
-      if (IsIdentifierActive(data[i])) {
-        PreaccHelper.addOutput(data[i]);
-      }
+      PreaccHelper.addOutput(data[i]);
     }
   }
 }
@@ -567,9 +627,7 @@ FORCEINLINE void SetPreaccOut(T&& data, const int size_x, const int size_y) {
   if (!PreaccActive) return;
   for (int i = 0; i < size_x; i++) {
     for (int j = 0; j < size_y; j++) {
-      if (IsIdentifierActive(data[i][j])) {
-        PreaccHelper.addOutput(data[i][j]);
-      }
+      PreaccHelper.addOutput(data[i][j]);
     }
   }
 }
@@ -676,6 +734,40 @@ FORCEINLINE void ResumePreaccumulation(bool wasActive) {
   SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
 }
 
+struct ErrorReport {
+  unsigned long ErrorCounter = 0;
+  std::ostream* out = &std::cout;
+};
+
+FORCEINLINE void ResetErrorCounter(ErrorReport& report) { report.ErrorCounter = 0; }
+
+FORCEINLINE void SetDebugReportFile(ErrorReport& report, std::ostream* output_file) { report.out = output_file; }
+
+FORCEINLINE unsigned long GetErrorCount(const ErrorReport& report) { return report.ErrorCounter; }
+
+#ifdef CODI_TAG_TAPE
+
+FORCEINLINE void SetTag(int tag) { AD::getTape().setCurTag(tag); }
+FORCEINLINE void ClearTagOnVariable(su2double& v) { AD::getTape().clearTagOnVariable(v); }
+
+static void tagErrorCallback(const int& correctTag, const int& wrongTag, void* userData) {
+  auto* report = static_cast<ErrorReport*>(userData);
+
+  report->ErrorCounter += 1;
+  *(report->out) << "Use of variable with bad tag '" << wrongTag << "', should be '" << correctTag << "'." << std::endl;
+}
+
+FORCEINLINE void SetTagErrorCallback(ErrorReport& report) {
+  AD::getTape().setTagErrorCallback(tagErrorCallback, &report);
+}
+
+#else
+FORCEINLINE void SetTag(int tag) {}
+FORCEINLINE void ClearTagOnVariable(su2double& v) {}
+FORCEINLINE void SetTagErrorCallback(ErrorReport report) {}
+
+#endif  // CODI_TAG_TAPE
+
 #endif  // CODI_REVERSE_TYPE
 
 void Initialize();
diff --git a/Common/include/basic_types/datatype_structure.hpp b/Common/include/basic_types/datatype_structure.hpp
index e16c223e8a69..9d5eb4559b02 100644
--- a/Common/include/basic_types/datatype_structure.hpp
+++ b/Common/include/basic_types/datatype_structure.hpp
@@ -89,8 +89,15 @@ void SetDerivative(su2double& data, const passivedouble& val);
 
 FORCEINLINE void SetValue(su2double& data, const passivedouble& val) { data.setValue(val); }
 
+FORCEINLINE passivedouble GetValue(const double& data) { return data; }
+
 FORCEINLINE passivedouble GetValue(const su2double& data) { return data.getValue(); }
 
+template <typename Expr>
+FORCEINLINE passivedouble GetValue(const codi::ExpressionInterface<double, Expr>& data) {
+  return data.cast().getValue();
+}
+
 FORCEINLINE void SetSecondary(su2double& data, const passivedouble& val) { data.setGradient(val); }
 
 FORCEINLINE void SetDerivative(su2double& data, const passivedouble& val) { data.setGradient(val); }
diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp
index 48625f20bf44..5509dd8d525f 100644
--- a/Common/include/code_config.hpp
+++ b/Common/include/code_config.hpp
@@ -110,6 +110,8 @@ using su2double = codi::RealReversePrimal;
 using su2double = codi::RealReversePrimalIndexGen<double, double, codi::ReuseIndexManager<int> >;
 #elif defined(CODI_PRIMAL_MULTIUSE_TAPE)
 using su2double = codi::RealReversePrimalIndex;
+#elif defined(CODI_TAG_TAPE)
+using su2double = codi::RealReverseTag;
 #else
 #error "Please define a CoDiPack tape."
 #endif
diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp
index f3912aa2a2de..05d64e25482a 100644
--- a/Common/include/geometry/dual_grid/CPoint.hpp
+++ b/Common/include/geometry/dual_grid/CPoint.hpp
@@ -113,9 +113,10 @@ class CPoint {
   su2activevector MaxLength;          /*!< \brief The maximum cell-center to cell-center length. */
   su2activevector RoughnessHeight;    /*!< \brief Roughness of the nearest wall. */
 
-  su2matrix<int> AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */
-  su2matrix<int>
-      AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */
+  su2matrix<AD::Identifier>
+      AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector before solver iteration. */
+  su2matrix<AD::Identifier>
+      AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after solver iteration. */
 
   /*!
    * \brief Allocate fields required by the minimal constructor.
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 223a53da742b..6f118ce29197 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -2543,6 +2543,29 @@ static const MapType<std::string, ENUM_DIRECTDIFF_VAR> DirectDiff_Var_Map = {
   MakePair("ELECTRIC_FIELD", D_EFIELD)
 };
 
+/*!
+ * \brief Types of tapes that can be checked in a tape debug run.
+ */
+enum class CHECK_TAPE_TYPE {
+  OBJECTIVE_FUNCTION,    /*!< \brief Tape that only includes dependencies and objective function calculation. */
+  FULL_SOLVER            /*!< \brief Tape that includes dependencies and all solvers. */
+};
+static const MapType<std::string, CHECK_TAPE_TYPE> CheckTapeType_Map = {
+    MakePair("OBJECTIVE_FUNCTION_TAPE", CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION)
+    MakePair("FULL_SOLVER_TAPE", CHECK_TAPE_TYPE::FULL_SOLVER)
+};
+
+/*!
+ * \brief Types of variables that can be checked for in a tape debug run.
+ */
+enum class CHECK_TAPE_VARIABLES {
+  SOLVER_VARIABLES,     /*!< \brief A (debug) tag will be assigned to solver/conservative variables. */
+  MESH_COORDINATES      /*!< \brief A (debug) tag will be assigned to solver/conservative variables and mesh coordinates. */
+};
+static const MapType<std::string, CHECK_TAPE_VARIABLES> CheckTapeVariables_Map = {
+    MakePair("SOLVER_VARIABLES", CHECK_TAPE_VARIABLES::SOLVER_VARIABLES)
+    MakePair("SOLVER_VARIABLES_AND_MESH_COORDINATES", CHECK_TAPE_VARIABLES::MESH_COORDINATES)
+};
 
 enum class RECORDING {
   CLEAR_INDICES,
@@ -2550,6 +2573,10 @@ enum class RECORDING {
   MESH_COORDS,
   MESH_DEFORM,
   SOLUTION_AND_MESH,
+  TAG_INIT_SOLVER_VARIABLES,
+  TAG_CHECK_SOLVER_VARIABLES,
+  TAG_INIT_SOLVER_AND_MESH,
+  TAG_CHECK_SOLVER_AND_MESH
 };
 
 /*!
diff --git a/Common/include/parallelization/vectorization.hpp b/Common/include/parallelization/vectorization.hpp
index f08017a9bee9..1e4b8b77f327 100644
--- a/Common/include/parallelization/vectorization.hpp
+++ b/Common/include/parallelization/vectorization.hpp
@@ -90,11 +90,11 @@ class Array : public CVecExpr<Array<Scalar_t, N>, Scalar_t> {
  public:
   using Scalar = Scalar_t;
   enum : size_t { Size = N };
-  enum : size_t { Align = Size * sizeof(Scalar) };
+  enum : size_t { Align = Size * alignof(Scalar) };
   static constexpr bool StoreAsRef = true;
 
  private:
-  alignas(Size * sizeof(Scalar)) Scalar x_[N];
+  alignas(Size * alignof(Scalar)) Scalar x_[N];
 
  public:
 #define ARRAY_BOILERPLATE                                                  \
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 9222ed7f8928..0d9e4a1fdd60 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1108,10 +1108,16 @@ void CConfig::SetConfig_Options() {
   addBoolOption("MULTIZONE", Multizone_Problem, NO);
   /*!\brief PHYSICAL_PROBLEM \n DESCRIPTION: Physical governing equations \n Options: see \link Solver_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("MULTIZONE_SOLVER", Kind_MZSolver, Multizone_Map, ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL);
-#ifdef CODI_REVERSE_TYPE
+#if defined (CODI_REVERSE_TYPE)
   const bool discAdjDefault = true;
+# if defined (CODI_TAG_TAPE)
+    DiscreteAdjointDebug = true;
+# else
+    DiscreteAdjointDebug = false;
+# endif
 #else
   const bool discAdjDefault = false;
+  DiscreteAdjointDebug = false;
 #endif
   /*!\brief MATH_PROBLEM  \n DESCRIPTION: Mathematical problem \n  Options: DIRECT, ADJOINT \ingroup Config*/
   addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault);
@@ -1186,6 +1192,8 @@ void CConfig::SetConfig_Options() {
   addBoolOption("WRT_VOLUME_OVERWRITE", Wrt_Volume_Overwrite, true);
   /*!\brief SYSTEM_MEASUREMENTS \n DESCRIPTION: System of measurements \n OPTIONS: see \link Measurements_Map \endlink \n DEFAULT: SI \ingroup Config*/
   addEnumOption("SYSTEM_MEASUREMENTS", SystemMeasurements, Measurements_Map, SI);
+  /*!\brief MULTIZONE_ADAPT_FILENAME \n DESCRIPTION: Append zone number to restart and solution filenames. \ingroup Config*/
+  addBoolOption("MULTIZONE_ADAPT_FILENAME", Multizone_Adapt_FileName, YES);
 
   /*!\par CONFIG_CATEGORY: FluidModel \ingroup Config*/
   /*!\brief FLUID_MODEL \n DESCRIPTION: Fluid model \n OPTIONS: See \link FluidModel_Map \endlink \n DEFAULT: STANDARD_AIR \ingroup Config*/
@@ -2793,6 +2801,12 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Preaccumulation in the AD mode. */
   addBoolOption("PREACC", AD_Preaccumulation, YES);
 
+  /* DESCRIPTION: Specify the tape which is checked in a tape debug run. */
+  addEnumOption("CHECK_TAPE_TYPE", AD_CheckTapeType, CheckTapeType_Map, CHECK_TAPE_TYPE::FULL_SOLVER);
+
+  /* DESCRIPTION: Specify the tape which is checked in a tape debug run. */
+  addEnumOption("CHECK_TAPE_VARIABLES", AD_CheckTapeVariables, CheckTapeVariables_Map, CHECK_TAPE_VARIABLES::SOLVER_VARIABLES);
+
   /*--- options that are used in the python optimization scripts. These have no effect on the c++ toolsuite ---*/
   /*!\par CONFIG_CATEGORY:Python Options\ingroup Config*/
 
@@ -3413,7 +3427,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
     /*---  Using default frequency of 250 for all files when steady, and 1 for unsteady. ---*/
     for (auto iVolumeFreq = 0; iVolumeFreq < nVolumeOutputFrequencies; iVolumeFreq++){
-      VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250;
+      if (Multizone_Problem && DiscreteAdjoint) {
+        VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : nOuterIter;
+      }
+      else {
+        VolumeOutputFrequencies[iVolumeFreq] = Time_Domain ? 1 : 250;
+      }
     }
   } else if (nVolumeOutputFrequencies < nVolumeOutputFiles) {
     /*--- If there are fewer frequencies than files, repeat the last frequency.
@@ -8378,7 +8397,7 @@ string CConfig::GetFilename(string filename, const string& ext, int timeIter) co
   filename = filename + string(ext);
 
   /*--- Append the zone number if multizone problems ---*/
-  if (Multizone_Problem)
+  if (Multizone_Problem && Multizone_Adapt_FileName)
     filename = GetMultizone_FileName(filename, GetiZone(), ext);
 
   /*--- Append the zone number if multiple instance problems ---*/
diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp
index 675fa232ab84..16da2e70b605 100644
--- a/Common/src/geometry/dual_grid/CPoint.cpp
+++ b/Common/src/geometry/dual_grid/CPoint.cpp
@@ -67,8 +67,8 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) {
   }
 
   if (config->GetDiscrete_Adjoint()) {
-    AD_InputIndex.resize(npoint, nDim) = 0;
-    AD_OutputIndex.resize(npoint, nDim) = 0;
+    AD_InputIndex.resize(npoint, nDim) = AD::GetPassiveIndex();
+    AD_OutputIndex.resize(npoint, nDim) = AD::GetPassiveIndex();
   }
 
   /*--- Multigrid structures. ---*/
diff --git a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
index dde873162087..8eec26d2e82b 100644
--- a/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
+++ b/SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
@@ -69,10 +69,10 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
   };
 
   /*!
-   * \brief Kinds of recordings (three different ones).
+   * \brief Kinds of recordings.
    */
   enum class Kind_Tape {
-    FULL_TAPE,                /*!< \brief Entire derivative information for a coupled adjoint
+    FULL_SOLVER_TAPE,         /*!< \brief Entire derivative information for a coupled adjoint
                                           solution update. */
     OBJECTIVE_FUNCTION_TAPE,  /*!< \brief Record only the dependence of the objective function
                                           w.r.t. solver variables (from all zones). */
@@ -96,7 +96,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
 
   bool eval_transfer = false;     /*!< \brief Evaluate the transfer section of the tape. */
   su2double ObjFunc;              /*!< \brief Value of the objective function. */
-  int ObjFunc_Index;              /*!< \brief Index of the value of the objective function. */
+  AD::Identifier ObjFunc_Index;   /*!< \brief Index of the value of the objective function. */
 
   CIteration*** direct_iteration;       /*!< \brief Array of pointers to the direct iterations. */
   COutput** direct_output;              /*!< \brief Array of pointers to the direct outputs. */
@@ -141,6 +141,16 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
    */
   void StartSolver() override;
 
+  /*!
+   * \brief [Overload] Launch the tape test mode for the discrete adjoint multizone solver.
+   */
+  void TapeTest ();
+
+  /*!
+   * \brief [Overload] Get error numbers after a tape test run of the discrete adjoint multizone solver.
+   */
+  int TapeTestGatherErrors(AD::ErrorReport& error_report) const;
+
   /*!
    * \brief Preprocess the multizone iteration
    */
diff --git a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp
index 7071c31750f9..d196eb1bff4f 100644
--- a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp
+++ b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp
@@ -57,7 +57,7 @@ class CDiscAdjSolver final : public CSolver {
   su2double Total_Sens_BPress;   /*!< \brief Total sensitivity to outlet pressure. */
   su2double Total_Sens_Density;  /*!< \brief Total sensitivity to initial density (incompressible). */
   su2double Total_Sens_ModVel;   /*!< \brief Total sensitivity to inlet velocity (incompressible). */
-  su2double Mach, Alpha, Beta, Pressure, Temperature, BPressure, ModVel;
+  su2double Mach, Alpha, Beta, Temperature, BPressure, ModVel;
   su2double TemperatureRad, Total_Sens_Temp_Rad;
 
   CDiscAdjVariable* nodes = nullptr;  /*!< \brief The highest level in the variable hierarchy this solver can safely use. */
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index d3673e1b763a..9ffb59baf2e8 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -671,12 +671,11 @@ void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& confi
 
     StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint));
     AD::SetPreaccOut(StrainMag(iPoint));
+    AD::EndPreacc();
 
     /*--- Max is not differentiable, so we not register them for preacc. ---*/
-    strainMax = max(strainMax, StrainMag(iPoint));
-    omegaMax = max(omegaMax, GeometryToolbox::Norm(3, Vorticity));
-
-    AD::EndPreacc();
+    strainMax = SU2_TYPE::GetValue(max(strainMax, StrainMag(iPoint)));
+    omegaMax = SU2_TYPE::GetValue(max(omegaMax, GeometryToolbox::Norm(3, Vorticity)));
   }
   END_SU2_OMP_FOR
 
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index 7088d2ea5a6e..3ff773038e11 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -92,8 +92,8 @@ class CVariable {
 
   MatrixType Solution_BGS_k;     /*!< \brief Old solution container for BGS iterations. */
 
-  su2matrix<int> AD_InputIndex;    /*!< \brief Indices of Solution variables in the adjoint vector. */
-  su2matrix<int> AD_OutputIndex;   /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */
+  su2matrix<AD::Identifier> AD_InputIndex;    /*!< \brief Indices of Solution variables in the adjoint vector before solver iteration. */
+  su2matrix<AD::Identifier> AD_OutputIndex;   /*!< \brief Indices of Solution variables in the adjoint vector after solver iteration. */
 
   VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables.
                                         Currently only streamwise periodic pressure-drop for massflow prescribed flows. */
@@ -118,7 +118,7 @@ class CVariable {
     assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class.");
   }
 
-  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<int>* ad_index = nullptr) {
+  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<AD::Identifier>* ad_index = nullptr) {
     const auto nPoint = variable.rows();
     SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads()))
     for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
@@ -133,11 +133,11 @@ class CVariable {
     END_SU2_OMP_FOR
   }
 
-  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<int>& ad_index) {
+  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<AD::Identifier>& ad_index) {
     RegisterContainer(input, variable, &ad_index);
   }
 
-  void RegisterContainer(bool input, su2activevector& variable, su2vector<int>* ad_index = nullptr) {
+  void RegisterContainer(bool input, su2activevector& variable, su2vector<AD::Identifier>* ad_index = nullptr) {
     const auto nPoint = variable.rows();
     SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads()))
     for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
@@ -2187,7 +2187,7 @@ class CVariable {
   }
 
   inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const {
-    int index = 0;
+    AD::Identifier index = AD::GetPassiveIndex();
     for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) {
       AD::SetIndex(index, Solution_time_n(iPoint, iVar));
       adj_sol[iVar] = AD::GetDerivative(index);
@@ -2195,7 +2195,7 @@ class CVariable {
   }
 
   inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const {
-    int index = 0;
+    AD::Identifier index = AD::GetPassiveIndex();
     for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++) {
       AD::SetIndex(index, Solution_time_n1(iPoint, iVar));
       adj_sol[iVar] = AD::GetDerivative(index);
diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
index c8dd0b1e98c3..42d2de2c1c59 100644
--- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
@@ -159,6 +159,18 @@ void CDiscAdjMultizoneDriver::Preprocess(unsigned long TimeIter) {
 
 void CDiscAdjMultizoneDriver::StartSolver() {
 
+  /*--- Start the debug recording mode for the discrete adjoint solver. ---*/
+
+  if (driver_config->GetDiscrete_Adjoint_Debug()) {
+
+    if (rank == MASTER_NODE) {
+      cout << "\nSU2_CFD_AD is compiled for debug mode recording. To resume the discrete adjoint solver, adjust -Dcodi-tape (-Dcodi-tape=JacobianLinear by default) and recompile." << endl;
+    }
+    Preprocess(0);
+    TapeTest();
+    return;
+  }
+
   const bool time_domain = driver_config->GetTime_Domain();
 
   /*--- Main external loop of the solver. Runs for the number of time steps required. ---*/
@@ -217,6 +229,96 @@ void CDiscAdjMultizoneDriver::StartSolver() {
 
 }
 
+void CDiscAdjMultizoneDriver::TapeTest() {
+
+  if (rank == MASTER_NODE) {
+    cout <<"\n---------------------------- Start Debug Run ----------------------------" << endl;
+  }
+
+  int total_errors = 0;
+  AD::ErrorReport error_report;
+  AD::SetTagErrorCallback(error_report);
+  std::ofstream out1("run1_process" + to_string(rank) + ".out");
+  std::ofstream out2("run2_process" + to_string(rank) + ".out");
+
+  AD::ResetErrorCounter(error_report);
+  AD::SetDebugReportFile(error_report, &out1);
+
+  /*--- This recording will assign the initial (same) tag to each registered variable.
+   *    During the recording, each dependent variable will be assigned the same tag. ---*/
+
+  if(driver_config->GetAD_CheckTapeType() == CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION) {
+    if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) {
+      if (rank == MASTER_NODE) cout << "\nChecking OBJECTIVE_FUNCTION_TAPE for SOLVER_VARIABLES_AND_MESH_COORDINATES." << endl;
+      SetRecording(RECORDING::TAG_INIT_SOLVER_AND_MESH, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0);
+    }
+    else {
+      if (rank == MASTER_NODE) cout << "\nChecking OBJECTIVE_FUNCTION_TAPE for SOLVER_VARIABLES." << endl;
+      SetRecording(RECORDING::TAG_INIT_SOLVER_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0);
+    }
+  }
+  else {
+    if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES) {
+      if (rank == MASTER_NODE) cout << "\nChecking FULL_SOLVER_TAPE for SOLVER_VARIABLES_AND_MESH_COORDINATES." << endl;
+      SetRecording(RECORDING::TAG_INIT_SOLVER_AND_MESH, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
+    }
+    else {
+      if (rank == MASTER_NODE) cout << "\nChecking FULL_SOLVER_TAPE for SOLVER_VARIABLES." << endl;
+      SetRecording(RECORDING::TAG_INIT_SOLVER_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
+    }
+  }
+  total_errors = TapeTestGatherErrors(error_report);
+
+  AD::ResetErrorCounter(error_report);
+  AD::SetDebugReportFile(error_report, &out2);
+
+  /*--- This recording repeats the initial recording with a different tag.
+   *    If a variable was used before it became dependent on the inputs, this variable will still carry the tag
+   *    from the initial recording and a mismatch with the "check" recording tag will throw an error.
+   *    In such a case, a possible reason could be that such a variable is set by a post-processing routine while
+   *    for a mathematically correct recording this dependency must be included earlier. ---*/
+
+  if(driver_config->GetAD_CheckTapeType() == CHECK_TAPE_TYPE::OBJECTIVE_FUNCTION) {
+    if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES)
+      SetRecording(RECORDING::TAG_CHECK_SOLVER_AND_MESH, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0);
+    else
+      SetRecording(RECORDING::TAG_CHECK_SOLVER_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0);
+  }
+  else {
+    if(driver_config->GetAD_CheckTapeVariables() == CHECK_TAPE_VARIABLES::MESH_COORDINATES)
+      SetRecording(RECORDING::TAG_CHECK_SOLVER_AND_MESH, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
+    else
+      SetRecording(RECORDING::TAG_CHECK_SOLVER_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
+  }
+  total_errors += TapeTestGatherErrors(error_report);
+
+  if (rank == MASTER_NODE) {
+    cout << "\n------------------------- Tape Test Run Summary -------------------------" << endl;
+    cout << "\nTotal number of tape inconsistencies: " << total_errors << endl;
+    cout << "\n--------------------------- End Tape Test Run ---------------------------" << endl;
+  }
+}
+
+int CDiscAdjMultizoneDriver::TapeTestGatherErrors(AD::ErrorReport& error_report) const {
+
+  int num_errors = AD::GetErrorCount(error_report);
+  int total_errors = 0;
+  std::vector<int> process_error(size);
+  SU2_MPI::Allreduce(&num_errors, &total_errors, 1, MPI_INT, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Gather(&num_errors, 1, MPI_INT, process_error.data(), 1, MPI_INT, MASTER_NODE, SU2_MPI::GetComm());
+
+  if (rank == MASTER_NODE) {
+    std::cout << "\nTotal number of detected tape inconsistencies: " << total_errors << std::endl;
+    if(total_errors > 0 && size > 1) {
+      std::cout << "\n";
+      for (int irank = 0; irank < size; irank++) {
+        std::cout << "Number of inconsistencies from process " << irank << ": " << process_error[irank] << std::endl;
+      }
+    }
+  }
+  return total_errors;
+}
+
 bool CDiscAdjMultizoneDriver::Iterate(unsigned short iZone, unsigned long iInnerIter, bool KrylovMode) {
 
   config_container[iZone]->SetInnerIter(iInnerIter);
@@ -310,6 +412,21 @@ void CDiscAdjMultizoneDriver::Run() {
   const bool time_domain = driver_config->GetTime_Domain();
   driver_config->Set_StartTime(SU2_MPI::Wtime());
 
+  /*--- Temporary warning because we need to test writing intermediate output to file (requires re-recording). ---*/
+  for(iZone = 0; iZone < nZone; iZone++) {
+    for (auto iVolumeFreq = 0; iVolumeFreq < config_container[iZone]->GetnVolumeOutputFrequencies(); iVolumeFreq++){
+      if (config_container[iZone]->GetVolumeOutputFrequency(iVolumeFreq) < nOuterIter) {
+        if (rank == MASTER_NODE) {
+          cout << "\nWARNING (iZone = " << iZone
+               << "): "
+                  "Writing out restart files during solver iterations is not tested for the discrete adjoint multizone solver.\n"
+                  "It is recommended to remove OUTPUT_WRT_FREQ from the config file, output files will be written when solver finalizes." << std::endl;
+        }
+        break;
+      }
+    }
+  }
+
   /*--- If the gradient of the objective function is 0 so are the adjoint variables.
    * Unless in unsteady problems where there are other contributions to the RHS. ---*/
 
@@ -352,8 +469,8 @@ void CDiscAdjMultizoneDriver::Run() {
      *    here. Otherwise, the whole tape of a coupled run will be created. ---*/
 
     if (RecordingState != RECORDING::SOLUTION_VARIABLES) {
-      SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_TAPE, ZONE_0);
-      SetRecording(RECORDING::SOLUTION_VARIABLES, Kind_Tape::FULL_TAPE, ZONE_0);
+      SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
+      SetRecording(RECORDING::SOLUTION_VARIABLES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
     }
 
     /*-- Start loop over zones. ---*/
@@ -499,11 +616,11 @@ void CDiscAdjMultizoneDriver::EvaluateSensitivities(unsigned long Iter, bool for
   /*--- SetRecording stores the computational graph on one iteration of the direct problem. Calling it with NONE
    *    as argument ensures that all information from a previous recording is removed. ---*/
 
-  SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_TAPE, ZONE_0);
+  SetRecording(RECORDING::CLEAR_INDICES, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
 
   /*--- Store the computational graph of one direct iteration with the mesh coordinates as input. ---*/
 
-  SetRecording(RECORDING::MESH_COORDS, Kind_Tape::FULL_TAPE, ZONE_0);
+  SetRecording(RECORDING::MESH_COORDS, Kind_Tape::FULL_SOLVER_TAPE, ZONE_0);
 
   /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution
    *    of the current iteration. The values are passed to the AD tool. ---*/
@@ -589,6 +706,10 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t
     case RECORDING::CLEAR_INDICES:      cout << "Clearing the computational graph." << endl; break;
     case RECORDING::MESH_COORDS:        cout << "Storing computational graph wrt MESH COORDINATES." << endl; break;
     case RECORDING::SOLUTION_VARIABLES: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << endl; break;
+    case RECORDING::TAG_INIT_SOLVER_VARIABLES:    cout << "Simulating recording with tag 1 on conservative variables." << endl; AD::SetTag(1); break;
+    case RECORDING::TAG_CHECK_SOLVER_VARIABLES:   cout << "Checking first recording with tag 2 on conservative variables." << endl; AD::SetTag(2); break;
+    case RECORDING::TAG_INIT_SOLVER_AND_MESH:     cout << "Simulating recording with tag 1 on conservative variables and mesh coordinates." << endl; AD::SetTag(1); break;
+    case RECORDING::TAG_CHECK_SOLVER_AND_MESH:    cout << "Checking first recording with tag 2 on conservative variables and mesh coordinates." << endl; AD::SetTag(2); break;
     default: break;
     }
   }
@@ -741,12 +862,12 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) {
   if (rank == MASTER_NODE) {
     AD::RegisterOutput(ObjFunc);
     AD::SetIndex(ObjFunc_Index, ObjFunc);
-    if (kind_recording == RECORDING::SOLUTION_VARIABLES) {
-      cout << " Objective function                   : " << ObjFunc;
-      if (driver_config->GetWrt_AD_Statistics()){
-        cout << " (" << ObjFunc_Index << ")\n";
-      }
-      cout << endl;
+    if (kind_recording == RECORDING::SOLUTION_VARIABLES ||
+        kind_recording == RECORDING::TAG_INIT_SOLVER_VARIABLES ||
+        kind_recording == RECORDING::TAG_CHECK_SOLVER_VARIABLES ||
+        kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH ||
+        kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH) {
+      cout << " Objective function                   : " << ObjFunc << endl;
     }
   }
 }
diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
index 7e51ff89c4b7..e0134adc0124 100644
--- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
@@ -404,7 +404,13 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge
 
   SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
-  if (kind_recording == RECORDING::SOLUTION_VARIABLES || kind_recording == RECORDING::SOLUTION_AND_MESH) {
+  if (kind_recording == RECORDING::SOLUTION_VARIABLES ||
+      kind_recording == RECORDING::TAG_INIT_SOLVER_VARIABLES ||
+      kind_recording == RECORDING::TAG_CHECK_SOLVER_VARIABLES ||
+      kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH ||
+      kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH ||
+      kind_recording == RECORDING::SOLUTION_AND_MESH) {
+
     /*--- Register flow and turbulent variables as input ---*/
 
     if (config[iZone]->GetFluidProblem()) {
@@ -429,9 +435,12 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge
     }
   }
 
-  if (kind_recording == RECORDING::MESH_COORDS || kind_recording == RECORDING::SOLUTION_AND_MESH) {
-    /*--- Register node coordinates as input ---*/
+  if (kind_recording == RECORDING::MESH_COORDS ||
+      kind_recording == RECORDING::SOLUTION_AND_MESH ||
+      kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH ||
+      kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH) {
 
+    /*--- Register node coordinates as input ---*/
     geometry0->RegisterCoordinates();
   }
 
diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
index 7a7e172a84d1..07caefe91c3c 100644
--- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
+++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
@@ -178,17 +178,21 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo
 
   /*--- Register farfield values as input ---*/
 
-  if((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS && !config->GetBoolTurbomachinery())) {
+  if((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS && !config->GetBoolTurbomachinery())){
 
     su2double Velocity_Ref = config->GetVelocity_Ref();
     Alpha                  = config->GetAoA()*PI_NUMBER/180.0;
     Beta                   = config->GetAoS()*PI_NUMBER/180.0;
     Mach                   = config->GetMach();
-    Pressure               = config->GetPressure_FreeStreamND();
-    Temperature            = config->GetTemperature_FreeStreamND();
+    /*--- Pressure and Temperature can be registered directly via their config file value
+     * (no further treatment required here). ---*/
+    su2double& Pressure    = config->GetPressure_FreeStreamND();
+    su2double& Temperature = config->GetTemperature_FreeStreamND();
 
     su2double SoundSpeed = 0.0;
 
+    // Treat Velocity_FreeStreamND config value as non-dependent (in debug mode)
+    AD::ClearTagOnVariable(config->GetVelocity_FreeStreamND()[0]);
     if (nDim == 2) { SoundSpeed = config->GetVelocity_FreeStreamND()[0]*Velocity_Ref/(cos(Alpha)*Mach); }
     if (nDim == 3) { SoundSpeed = config->GetVelocity_FreeStreamND()[0]*Velocity_Ref/(cos(Alpha)*cos(Beta)*Mach); }
 
@@ -211,11 +215,8 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo
       config->GetVelocity_FreeStreamND()[2] = sin(Alpha)*cos(Beta)*Mach*SoundSpeed/Velocity_Ref;
     }
 
-    config->SetTemperature_FreeStreamND(Temperature);
     direct_solver->SetTemperature_Inf(Temperature);
-    config->SetPressure_FreeStreamND(Pressure);
     direct_solver->SetPressure_Inf(Pressure);
-
   }
 
   if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && config->GetBoolTurbomachinery()){
@@ -399,6 +400,9 @@ void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *conf
   if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && !config->GetBoolTurbomachinery()) {
     su2double Local_Sens_Press, Local_Sens_Temp, Local_Sens_AoA, Local_Sens_Mach;
 
+    su2double& Pressure    = config->GetPressure_FreeStreamND();
+    su2double& Temperature = config->GetTemperature_FreeStreamND();
+
     Local_Sens_Mach  = SU2_TYPE::GetDerivative(Mach);
     Local_Sens_AoA   = SU2_TYPE::GetDerivative(Alpha);
     Local_Sens_Temp  = SU2_TYPE::GetDerivative(Temperature);
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index f2e6e4b6d43f..796e009f181b 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -4108,7 +4108,7 @@ void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *conf
 
   unsigned short iMarker, iDim;
   unsigned long iVertex, iPoint;
-  int index;
+  AD::Identifier index;
 
   /*--- Loop over all the markers ---*/
   for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp
index 75bf88822826..3370fbb7a104 100644
--- a/SU2_CFD/src/variables/CVariable.cpp
+++ b/SU2_CFD/src/variables/CVariable.cpp
@@ -69,8 +69,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
       External.resize(nPoint,nVar) = su2double(0.0);
 
     if (!adjoint) {
-      AD_InputIndex.resize(nPoint,nVar) = -1;
-      AD_OutputIndex.resize(nPoint,nVar) = -1;
+      AD_InputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex();
+      AD_OutputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex();
     }
   }
 
diff --git a/TestCases/TestCase.py b/TestCases/TestCase.py
index 03b431ba3a7d..2da91f390462 100644
--- a/TestCases/TestCase.py
+++ b/TestCases/TestCase.py
@@ -46,6 +46,7 @@ def parse_args(description: str):
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument('--tsan', action='store_true', help='Run thread sanitizer tests. Requires a tsan-enabled SU2 build.')
     parser.add_argument('--asan', action='store_true', help='Run address sanitizer tests. Requires an asan-enabled SU2 build.')
+    parser.add_argument('--tapetests', action='store_true', help='Run discrete adjoint tests in tape debug mode. Requires a SU2_CFD_AD build with Tag reverse type.')
     return parser.parse_args()
 
 class TestCase:
@@ -107,10 +108,14 @@ def __init__(self,tag_in):
         # multizone problem
         self.multizone = False
 
+        # Indicate tapetest mode
+        self.enabled_with_tapetests = False
+
         # The test condition. These must be set after initialization
         self.test_iter = 1
         self.ntest_vals = 4
         self.test_vals = []
+        self.tapetest_vals = []
         self.test_vals_aarch64 = []
         self.cpu_arch = platform.machine().casefold()
         self.enabled_on_cpu_arch = ["x86_64","amd64","aarch64","arm64"]
@@ -119,6 +124,7 @@ def __init__(self,tag_in):
         self.command = self.Command()
         self.timeout = 0
         self.tol = 0.0
+        self.tapetest_tol = 0
         self.tol_file_percent = 0.0
         self.comp_threshold = 0.0
 
@@ -127,9 +133,10 @@ def __init__(self,tag_in):
         self.reference_file_aarch64 = ""
         self.test_file      = "of_grad.dat"
 
-    def run_test(self, with_tsan=False, with_asan=False):
+    def run_test(self, with_tsan=False, with_asan=False, with_tapetests=False):
 
-        if not self.is_enabled(with_tsan, with_asan):
+        # Check whether this test is valid and can be continued
+        if not self.is_enabled(with_tsan, with_asan, with_tapetests):
             return True
 
         print('==================== Start Test: %s ===================='%self.tag)
@@ -138,6 +145,7 @@ def run_test(self, with_tsan=False, with_asan=False):
         timed_out    = False
         iter_missing = True
         start_solver = True
+        tapetest_out = True
 
         # if root, add flag to mpirun
         self.command.allow_mpi_as_root()
@@ -188,10 +196,10 @@ def run_test(self, with_tsan=False, with_asan=False):
         delta_vals = []
         sim_vals = []
 
-        if not with_tsan and not with_asan: # sanitizer findings result in non-zero return code, no need to examine the output
+        if not with_tsan and not with_asan and not with_tapetests: # Sanitizer findings result in non-zero return code, no need to examine the output. Tapetest output is examined separately.
             # Examine the output
-            f = open(logfilename,'r')
-            output = f.readlines()
+            with open(logfilename,'r') as f:
+                output = f.readlines()
             if not timed_out and len(self.test_vals) != 0:
                 start_solver = False
                 for line in output:
@@ -240,6 +248,32 @@ def run_test(self, with_tsan=False, with_asan=False):
         #for j in output:
         #  print(j)
 
+        if with_tapetests and self.enabled_with_tapetests: # examine the tapetest output
+            with open(logfilename,'r') as f:
+                output = f.readlines()
+            if not timed_out and len(self.tapetest_vals) != 0:
+                tapetest_out = False
+                for line in output:
+                    if not tapetest_out: # Don't bother parsing anything before "Total number of tape inconsistencies"
+                        if line.find('Total number of tape inconsistencies:') > -1:
+                            tapetest_out = True
+                            raw_data = line.split(':') # Split line into description and the string representing the number of errors
+                            data = raw_data[1].strip() # Clear the string representing the number of errors (for now, expecting a single integer, but zone-wise error numbers are planned)
+
+                            if not len(self.tapetest_vals)==len(data):   # something went wrong... probably bad input
+                                print("Error in tapetest_vals!")
+                                passed = False
+                                break
+                            for j in range(len(data)):
+                                sim_vals.append( int(data[j]) )
+                                delta_vals.append( abs(int(data[j])-self.tapetest_vals[j]) )
+                                if delta_vals[j] > self.tapetest_tol:
+                                    exceed_tol = True
+                                    passed     = False
+
+                if not tapetest_out:
+                    passed = False
+
 
         process.communicate()
         if process.returncode != 0:
@@ -259,15 +293,18 @@ def run_test(self, with_tsan=False, with_asan=False):
         if exceed_tol:
             print('ERROR: Difference between computed input and test_vals exceeded tolerance. TOL=%f'%self.tol)
 
-        if not start_solver:
+        if not start_solver and not with_tapetests:
             print('ERROR: The code was not able to get to the "Begin solver" section.')
 
-        if not with_tsan and not with_asan and iter_missing:
+        if not tapetest_out and with_tapetests:
+            print('ERROR: The code was not able to get to the "Total number of tape inconsistencies" line.')
+
+        if not with_tsan and not with_asan and not with_tapetests and iter_missing:
             print('ERROR: The iteration number %d could not be found.'%self.test_iter)
 
         print('CPU architecture=%s' % self.cpu_arch)
 
-        if len(self.test_vals) != 0:
+        if len(self.test_vals) != 0 and not with_tapetests:
             print('test_iter=%d' % self.test_iter)
 
             print_vals(self.test_vals, name="test_vals (stored)")
@@ -283,9 +320,9 @@ def run_test(self, with_tsan=False, with_asan=False):
         os.chdir(workdir)
         return passed
 
-    def run_filediff(self, with_tsan=False, with_asan=False):
+    def run_filediff(self, with_tsan=False, with_asan=False, with_tapetests=False):
 
-        if not self.is_enabled(with_tsan, with_asan):
+        if not self.is_enabled(with_tsan, with_asan, with_tapetests):
             return True
 
         print('==================== Start Test: %s ===================='%self.tag)
@@ -480,9 +517,9 @@ def run_filediff(self, with_tsan=False, with_asan=False):
         os.chdir(workdir)
         return passed
 
-    def run_opt(self, with_tsan=False, with_asan=False):
+    def run_opt(self, with_tsan=False, with_asan=False, with_tapetests=False):
 
-        if not self.is_enabled(with_tsan, with_asan):
+        if not self.is_enabled(with_tsan, with_asan, with_tapetests):
             return True
 
         print('==================== Start Test: %s ===================='%self.tag)
@@ -605,9 +642,9 @@ def run_opt(self, with_tsan=False, with_asan=False):
         os.chdir(workdir)
         return passed
 
-    def run_geo(self, with_tsan=False, with_asan=False):
+    def run_geo(self, with_tsan=False, with_asan=False, with_tapetests=False):
 
-        if not self.is_enabled(with_tsan, with_asan):
+        if not self.is_enabled(with_tsan, with_asan, with_tapetests):
             return True
 
         print('==================== Start Test: %s ===================='%self.tag)
@@ -696,7 +733,7 @@ def run_geo(self, with_tsan=False, with_asan=False):
                         delta_vals.append( abs(float(data[j])-self.test_vals[j]) )
                         if delta_vals[j] > self.tol:
                             exceed_tol = True
-                            passed     = False
+                            passed = False
                 else:
                     iter_missing = True
 
@@ -745,9 +782,9 @@ def run_geo(self, with_tsan=False, with_asan=False):
         os.chdir(workdir)
         return passed
 
-    def run_def(self, with_tsan=False, with_asan=False):
+    def run_def(self, with_tsan=False, with_asan=False, with_tapetests=False):
 
-        if not self.is_enabled(with_tsan, with_asan):
+        if not self.is_enabled(with_tsan, with_asan, with_tapetests):
             return True
 
         print('==================== Start Test: %s ===================='%self.tag)
@@ -970,19 +1007,24 @@ def disable_restart(self):
 
         return
 
-    def is_enabled(self, with_tsan=False, with_asan=False):
+    def is_enabled(self, with_tsan=False, with_asan=False, with_tapetests=False):
         is_enabled_on_arch = self.cpu_arch in self.enabled_on_cpu_arch
 
         if not is_enabled_on_arch:
             print('Ignoring test "%s" because it is not enabled for the current CPU architecture: %s' % (self.tag, self.cpu_arch))
 
+        # A test case is valid to be continued if neither of the special modes (tsan, asan, tapetests) is active, or if so, the corresponding test case is enabled, too.
         tsan_compatible = not with_tsan or self.enabled_with_tsan
         asan_compatible = not with_asan or self.enabled_with_asan
+        tapetests_compatible = not with_tapetests or self.enabled_with_tapetests
 
         if not tsan_compatible:
             print('Ignoring test "%s" because it is not enabled to run with the thread sanitizer.' % self.tag)
 
-        return is_enabled_on_arch and tsan_compatible and asan_compatible
+        if not tapetests_compatible:
+            print('Ignoring test "%s" because it is not enabled to run a test of the tape.' % self.tag)
+
+        return is_enabled_on_arch and tsan_compatible and asan_compatible and tapetests_compatible and tapetests_compatible
 
     def adjust_test_data(self):
 
diff --git a/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg
new file mode 100644
index 000000000000..85d1d03c5a24
--- /dev/null
+++ b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_discadj_multizone.cfg
@@ -0,0 +1,112 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Adjoint transonic inviscid flow around a NACA0012 airfoil  %
+% Author: Thomas D. Economon                                                   %
+% Institution: Stanford University                                             %
+% Date: 2011.11.02                                                             %
+% File Version 8.2.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+MULTIZONE= YES
+MULTIZONE_ADAPT_FILENAME= NO
+SOLVER= EULER
+MATH_PROBLEM= DISCRETE_ADJOINT
+RESTART_SOL= NO
+READ_BINARY_RESTART= NO
+
+% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
+%
+MACH_NUMBER= 0.8
+AOA= 1.25
+FREESTREAM_PRESSURE= 101325.0
+FREESTREAM_TEMPERATURE= 288.15
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+
+% ----------------------- BOUNDARY CONDITION DEFINITION -----------------------%
+%
+MARKER_EULER= ( airfoil )
+MARKER_FAR= ( farfield )
+MARKER_PLOTTING= ( airfoil )
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------%
+%
+NUM_METHOD_GRAD= GREEN_GAUSS
+OBJECTIVE_FUNCTION= DRAG
+CFL_NUMBER= 5.0
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+OUTER_ITER= 150
+INNER_ITER= 1
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= JST
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% ---------------- ADJOINT-FLOW NUMERICAL METHOD DEFINITION -------------------%
+%
+CONV_NUM_METHOD_ADJFLOW= ROE
+MUSCL_ADJFLOW= YES
+SLOPE_LIMITER_ADJFLOW= NONE
+ADJ_SHARP_LIMITER_COEFF= 3.0
+ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
+CFL_REDUCTION_ADJFLOW= 0.5
+TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
+
+% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------%
+%
+DV_KIND= HICKS_HENNE
+DV_MARKER= ( airfoil )
+DV_PARAM= ( 1, 0.5 )
+DV_VALUE= 0.01
+
+% ------------------------ GRID DEFORMATION PARAMETERS ------------------------%
+%
+DEFORM_STIFFNESS_TYPE= WALL_DISTANCE
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_RESIDUAL_MINVAL= -12
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= mesh_NACA0012_inv.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out.su2
+SOLUTION_FILENAME= solution_flow.dat
+SOLUTION_ADJ_FILENAME= solution_adj.dat
+TABULAR_FORMAT= TECPLOT
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow.dat
+RESTART_ADJ_FILENAME= restart_adj.dat
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad.dat
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 250
+SCREEN_OUTPUT= (OUTER_ITER, BGS_RES[0])
+OUTPUT_FILES=(RESTART_ASCII)
+
+% --------------------- OPTIMAL SHAPE DESIGN DEFINITION -----------------------%
+%
+OPT_OBJECTIVE= DRAG * 0.001
+OPT_CONSTRAINT= ( LIFT > 0.327 ) * 0.001; ( MOMENT_Z > 0.0 ) * 0.001; ( AIRFOIL_THICKNESS > 0.12 ) * 0.001
+DEFINITION_DV= ( 30, 1.0 | airfoil | 0, 0.05 ); ( 30, 1.0 | airfoil | 0, 0.10 ); ( 30, 1.0 | airfoil | 0, 0.15 ); ( 30, 1.0 | airfoil | 0, 0.20 ); ( 30, 1.0 | airfoil | 0, 0.25 ); ( 30, 1.0 | airfoil | 0, 0.30 ); ( 30, 1.0 | airfoil | 0, 0.35 ); ( 30, 1.0 | airfoil | 0, 0.40 ); ( 30, 1.0 | airfoil | 0, 0.45 ); ( 30, 1.0 | airfoil | 0, 0.50 ); ( 30, 1.0 | airfoil | 0, 0.55 ); ( 30, 1.0 | airfoil | 0, 0.60 ); ( 30, 1.0 | airfoil | 0, 0.65 ); ( 30, 1.0 | airfoil | 0, 0.70 ); ( 30, 1.0 | airfoil | 0, 0.75 ); ( 30, 1.0 | airfoil | 0, 0.80 ); ( 30, 1.0 | airfoil | 0, 0.85 ); ( 30, 1.0 | airfoil | 0, 0.90 ); ( 30, 1.0 | airfoil | 0, 0.95 ); ( 30, 1.0 | airfoil | 1, 0.05 ); ( 30, 1.0 | airfoil | 1, 0.10 ); ( 30, 1.0 | airfoil | 1, 0.15 ); ( 30, 1.0 | airfoil | 1, 0.20 ); ( 30, 1.0 | airfoil | 1, 0.25 ); ( 30, 1.0 | airfoil | 1, 0.30 ); ( 30, 1.0 | airfoil | 1, 0.35 ); ( 30, 1.0 | airfoil | 1, 0.40 ); ( 30, 1.0 | airfoil | 1, 0.45 ); ( 30, 1.0 | airfoil | 1, 0.50 ); ( 30, 1.0 | airfoil | 1, 0.55 ); ( 30, 1.0 | airfoil | 1, 0.60 ); ( 30, 1.0 | airfoil | 1, 0.65 ); ( 30, 1.0 | airfoil | 1, 0.70 ); ( 30, 1.0 | airfoil | 1, 0.75 ); ( 30, 1.0 | airfoil | 1, 0.80 ); ( 30, 1.0 | airfoil | 1, 0.85 ); ( 30, 1.0 | airfoil | 1, 0.90 ); ( 30, 1.0 | airfoil | 1, 0.95 )
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 89f64d54c22b..687a459cb870 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -53,6 +53,16 @@ def main():
     discadj_naca0012.test_vals = [-3.562611, -8.932639, -0.000000, 0.005608]
     test_list.append(discadj_naca0012)
 
+    # Inviscid NACA0012 (via discadj multizone driver)
+    discadj_naca0012_via_mz = TestCase('discadj_naca0012_via_mz')
+    discadj_naca0012_via_mz.cfg_dir = "cont_adj_euler/naca0012"
+    discadj_naca0012_via_mz.cfg_file = "inv_NACA0012_discadj_multizone.cfg"
+    discadj_naca0012_via_mz.test_iter = 100
+    discadj_naca0012_via_mz.test_vals = [-3.563784, -5.975640, -6.326231, -8.929567]
+    discadj_naca0012_via_mz.enabled_with_tapetests = True
+    discadj_naca0012_via_mz.tapetest_vals = [0]
+    test_list.append(discadj_naca0012_via_mz)
+
     # Inviscid Cylinder 3D (multiple markers)
     discadj_cylinder3D           = TestCase('discadj_cylinder3D')
     discadj_cylinder3D.cfg_dir   = "disc_adj_euler/cylinder3D"
@@ -226,7 +236,7 @@ def main():
         if test.tol == 0.0:
             test.tol = 0.00001
 
-    pass_list = [ test.run_test(args.tsan, args.asan) for test in test_list ]
+    pass_list = [ test.run_test(args.tsan, args.asan, args.tapetests) for test in test_list ]
 
     ###################################
     ### Coupled RHT-CFD Adjoint     ###
@@ -243,7 +253,7 @@ def main():
     discadj_rht.reference_file_aarch64 = "of_grad_cd_aarch64.csv.ref"
     discadj_rht.test_file      = "of_grad_cd.csv"
     discadj_rht.enabled_with_asan = False
-    pass_list.append(discadj_rht.run_filediff(args.tsan, args.asan))
+    pass_list.append(discadj_rht.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(discadj_rht)
 
     ######################################
@@ -261,7 +271,7 @@ def main():
     discadj_euler_py.reference_file_aarch64 = "of_grad_cd_disc_aarch64.dat.ref"
     discadj_euler_py.test_file = "of_grad_cd.dat"
     discadj_euler_py.enabled_with_asan = False
-    pass_list.append(discadj_euler_py.run_filediff(args.tsan, args.asan))
+    pass_list.append(discadj_euler_py.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(discadj_euler_py)
 
     # test discrete_adjoint with multiple ffd boxes
@@ -275,7 +285,7 @@ def main():
     discadj_multiple_ffd_py.reference_file_aarch64 = "of_grad_cd_aarch64.dat.ref"
     discadj_multiple_ffd_py.test_file = "of_grad_cd.dat"
     discadj_multiple_ffd_py.enabled_with_asan = False
-    pass_list.append(discadj_multiple_ffd_py.run_filediff(args.tsan, args.asan))
+    pass_list.append(discadj_multiple_ffd_py.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(discadj_multiple_ffd_py)
 
     # test direct_differentiation.py
@@ -289,7 +299,7 @@ def main():
     directdiff_euler_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref"
     directdiff_euler_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat"
     directdiff_euler_py.enabled_with_asan = False
-    pass_list.append(directdiff_euler_py.run_filediff(args.tsan, args.asan))
+    pass_list.append(directdiff_euler_py.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(directdiff_euler_py)
 
     # test direct_differentiation.py with multiple ffd boxes
@@ -303,7 +313,7 @@ def main():
     directdiff_multiple_ffd_py.reference_file_aarch64 = "of_grad_directdiff_aarch64.dat.ref"
     directdiff_multiple_ffd_py.test_file = "DIRECTDIFF/of_grad_directdiff.dat"
     directdiff_multiple_ffd_py.enabled_with_asan = False
-    pass_list.append(directdiff_multiple_ffd_py.run_filediff(args.tsan, args.asan))
+    pass_list.append(directdiff_multiple_ffd_py.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(directdiff_multiple_ffd_py)
 
     # test continuous_adjoint.py, with multiple objectives
@@ -330,7 +340,7 @@ def main():
     pywrapper_FEA_AD_FlowLoad.new_output    = False
     pywrapper_FEA_AD_FlowLoad.enabled_with_asan = False
     test_list.append(pywrapper_FEA_AD_FlowLoad)
-    pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test(args.tsan, args.asan))
+    pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test(args.tsan, args.asan, args.tapetests))
 
     # Flow AD Mesh Displacement Sensitivity
     pywrapper_CFD_AD_MeshDisp               = TestCase('pywrapper_CFD_AD_MeshDisp')
@@ -344,7 +354,7 @@ def main():
     pywrapper_CFD_AD_MeshDisp.new_output    = False
     pywrapper_CFD_AD_MeshDisp.enabled_with_asan = False
     test_list.append(pywrapper_CFD_AD_MeshDisp)
-    pass_list.append(pywrapper_CFD_AD_MeshDisp.run_test(args.tsan, args.asan))
+    pass_list.append(pywrapper_CFD_AD_MeshDisp.run_test(args.tsan, args.asan, args.tapetests))
 
 
     ###################################
@@ -361,7 +371,7 @@ def main():
     grad_smooth_naca0012.reference_file_aarch64 = "of_hess_aarch64.dat.ref"
     grad_smooth_naca0012.test_file = "of_hess.dat"
     grad_smooth_naca0012.enabled_with_asan = False
-    pass_list.append(grad_smooth_naca0012.run_filediff(args.tsan, args.asan))
+    pass_list.append(grad_smooth_naca0012.run_filediff(args.tsan, args.asan, args.tapetests))
     test_list.append(grad_smooth_naca0012)
 
     # Tests summary
diff --git a/config_template.cfg b/config_template.cfg
index 32e6369e6e33..4e2d456170fe 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -2041,13 +2041,20 @@ TOTAL_DV_PENALTY= 0.0
 %
 % Parameters for the corresponding OF (allowed stress and KS multiplier).
 STRESS_PENALTY_PARAM= (1.0, 10.0)
-%
 
 % ---------------- AUTOMATIC DIFFERENTIATION -------------------%
 %
 % Preaccumulation in the AD mode.
 PREACC= YES
 
+% ---------------- AUTOMATIC DIFFERENTIATION (TAPE TEST DEBUG MODE) -------------------%
+%
+% Specify the kind of tape that is checked (OBJECTIVE_FUNCTION_TAPE, FULL_SOLVER_TAPE)
+CHECK_TAPE_TYPE = FULL_SOLVER_TAPE
+%
+% Specify the variables for which the tape is checked (SOLVER_VARIABLES, SOLVER_VARIABLES_AND_MESH_COORDINATES)
+CHECK_TAPE_VARIABLES = SOLVER_VARIABLES
+
 % ---------------- MESH DEFORMATION PARAMETERS (NEW SOLVER) -------------------%
 %
 % Use the reformatted pseudo-elastic solver for grid deformation
@@ -2376,6 +2383,9 @@ TABULAR_FORMAT= CSV
 % Set <ofstream>.precision(value) to specified value for SU2_DOT and HISTORY output. Useful for exact gradient validation.
 OUTPUT_PRECISION= 10
 %
+% For multizone problems, extend solution and restart filenames automatically by zone number
+MULTIZONE_ADAPT_FILENAME= YES
+%
 % Files to output
 % Possible formats : (TECPLOT_ASCII, TECPLOT, SURFACE_TECPLOT_ASCII,
 %  SURFACE_TECPLOT, CSV, SURFACE_CSV, PARAVIEW_ASCII, PARAVIEW_LEGACY, SURFACE_PARAVIEW_ASCII,
diff --git a/meson.build b/meson.build
index 443cb1f4e543..1ed26b6809e1 100644
--- a/meson.build
+++ b/meson.build
@@ -101,6 +101,8 @@ if get_option('enable-autodiff') and not omp
     codi_rev_args += '-DCODI_PRIMAL_REUSE_TAPE'
   elif get_option('codi-tape') == 'PrimalMultiUse'
     codi_rev_args += '-DCODI_PRIMAL_MULTIUSE_TAPE'
+  elif get_option('codi-tape') == 'Tag'
+    codi_rev_args += '-DCODI_TAG_TAPE'
   else
     error('Invalid CoDiPack tape choice @0@'.format(get_option('codi-tape')))
   endif
diff --git a/meson_options.txt b/meson_options.txt
index 830b5f6574ee..58305769f163 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -24,7 +24,7 @@ option('enable-coolprop',  type : 'boolean', value : false, description: 'enable
 option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable MLPCpp support')
 option('enable-gprof', type : 'boolean', value : false, description: 'enable profiling through gprof')
 option('opdi-backend', type : 'combo', choices : ['auto', 'macro', 'ompt'], value : 'auto', description: 'OpDiLib backend choice')
-option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse'], value : 'JacobianLinear', description: 'CoDiPack tape choice')
+option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse', 'Tag'], value : 'JacobianLinear', description: 'CoDiPack tape choice')
 option('opdi-shared-read-opt', type : 'boolean', value : true, description : 'OpDiLib shared reading optimization')
 option('librom_root', type : 'string', value : '', description: 'libROM base directory')
 option('enable-librom', type : 'boolean', value : false, description: 'enable LLNL libROM support')