Index: /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp	(revision 26432)
+++ /issm/trunk-jpl/src/c/toolkits/codipack/CoDiPackCommon.hpp	(revision 26433)
@@ -55,5 +55,5 @@
 template<typename Real, typename Data>
 inline void registerVector(Real* vec, Data* dataVec, int n) {
-  typename Real::TapeType& tape = Real::getTape();
+  typename Real::Tape& tape = Real::getTape();
 
   for(int i = 0; i < n; ++i) {
Index: /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 26432)
+++ /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 26433)
@@ -254,16 +254,16 @@
 
 #ifdef _HAVE_CODIPACK_
-void SolverxSeq_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
+using Tape = typename IssmDouble::Tape;
+using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
+void SolverxSeq_codi_b(Tape* tape,void* data_in, AccessInterface* ra) {/*{{{*/
 
 	/*recast data_in and tape*/
-	codi::DataStore* data = (codi::DataStore*)data_in;
-	//IssmDouble::Tape& tape = (IssmDouble::Tape&)tape_in;
-	IssmDouble::Tape& tape = IssmDouble::getTape();
+	codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexATrans;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
 
@@ -275,17 +275,16 @@
   data->getData(n);
 
-
   // create the adjoint vector for x and reset the adjoint values on the tape
   IssmDouble::Gradient* adjX = xNew<IssmDouble::Gradient>(n);
-  getVectorAdjoint(tape, indexX, adjX, n);
+  getVectorAdjoint(*tape, indexX, adjX, n);
 
   IssmDouble::Gradient* sol  = xNew<IssmDouble::Gradient>(n);
   SolverxSeq(sol, valueATrans, adjX, n);
 
-  updateVectorAdjoint(tape, indexB, sol, n);
+  updateVectorAdjoint(*tape, indexB, sol, n);
   for(int i=0; i<n; ++i) {
     for (int j=0; j<n; ++j) {
       // we access the transposed matrix here because we stored the indices in a transposed way
-      updateAdjoint(tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
+      updateAdjoint(*tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
     }
   }
@@ -295,14 +294,14 @@
 }
 /*}}}*/
-void SolverxSeq_codi_delete(void* tape_in,void* data_in) {/*{{{*/
+void SolverxSeq_codi_delete(Tape* tape,void* data_in) {/*{{{*/
 
 	/*recast data_in*/
-	codi::DataStore* data = (codi::DataStore*)data_in;
+	codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexATrans;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
 
@@ -325,12 +324,12 @@
 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
   IssmDouble::Tape& tape = IssmDouble::getTape();
-  codi::DataStore* dataHandler = NULL;
+  codi::ExternalFunctionUserData* dataHandler = NULL;
 
   if(tape.isActive()) {
-    dataHandler = new codi::DataStore();
+    dataHandler = new codi::ExternalFunctionUserData();
 
     // create the index vector and the double data for A and B
     IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(n*n);
-    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(n*n);
+    IssmDouble::Identifier* indexATrans = xNew<IssmDouble::Identifier>(n*n);
 
     // read the data for matrix in a transposed fashion
@@ -342,5 +341,5 @@
 
     // read the data from B (primal values are not required vor B
-    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexB = xNew<IssmDouble::Identifier>(n);
     getVectorGradData(B, indexB, n);
 
@@ -366,5 +365,5 @@
   if(tape.isActive()) {
     // create the index vector X and register x as active variables
-    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexX = xNew<IssmDouble::Identifier>(n);
     registerVector(X, indexX, n);
 
@@ -375,5 +374,5 @@
     dataHandler->addData(n);
 
-    tape.pushExternalFunctionHandle(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete);
+	 tape.pushExternalFunction(codi::ExternalFunction<Tape>::create(&SolverxSeq_codi_b,(void*)dataHandler, &SolverxSeq_codi_delete));
   }
   else{
Index: /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 26432)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 26433)
@@ -485,19 +485,19 @@
 
 #ifdef _HAVE_CODIPACK_
-void MumpsSolve_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
+using Tape = typename IssmDouble::Tape;
+using AccessInterface = codi::VectorAccessInterface<typename Tape::Real, typename Tape::Identifier>;
+
+void MumpsSolve_codi_b(Tape* tape, void* data_in, AccessInterface* ra) {/*{{{*/
 
 	/*recast data_in and tape*/
-  codi::DataStore* data = (codi::DataStore*)data_in;
-  //IssmDouble::Tape& tape = (IssmDouble::Tape&)tape_in;
-  IssmDouble::Tape& tape = IssmDouble::getTape();
-
+  codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
+  IssmDouble::Identifier* indexATrans;
   int* irnATrans;
   int* jcnATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
   int nnz;
@@ -519,10 +519,10 @@
   // create the adjoint vector for x and reset the adjoint values on the tape
   IssmDouble::Gradient* adjX = xNew<IssmDouble::Gradient>(n);
-  getVectorAdjoint(tape, indexX, adjX, n);
+  getVectorAdjoint(*tape, indexX, adjX, n);
 
   MumpsSolve(n, nnz, local_nnz, irnATrans, jcnATrans, valueATrans, adjX, parameters);
   // adjX contains now the solution
 
-  updateVectorAdjoint(tape, indexB, adjX, n);
+  updateVectorAdjoint(*tape, indexB, adjX, n);
 
   // bcast dp_y (the solution of the forward system)
@@ -534,5 +534,5 @@
     // we access the transposed matrix here because we stored the indices in a transposed way
     // -1 is substracted because jcn and irn are stored with fortran indexing
-    updateAdjoint(tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
+    updateAdjoint(*tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
   }
 
@@ -540,16 +540,16 @@
 }
 /*}}}*/
-void MumpsSolve_codi_delete(void* tape_in,void* data_in) {/*{{{*/
+void MumpsSolve_codi_delete(Tape* tape,void* data_in) {/*{{{*/
 
 	/*recast data_in*/
-	codi::DataStore* data = (codi::DataStore*)data_in;
+	codi::ExternalFunctionUserData* data = (codi::ExternalFunctionUserData*)data_in;
 
   IssmDouble::Real* valueATrans;
-  IssmDouble::GradientData* indexATrans;
+  IssmDouble::Identifier* indexATrans;
   int* irnATrans;
   int* jcnATrans;
-  IssmDouble::GradientData* indexB;
+  IssmDouble::Identifier* indexB;
   IssmDouble::Real* valueX;
-  IssmDouble::GradientData* indexX;
+  IssmDouble::Identifier* indexX;
   int n;
   int nnz;
@@ -582,12 +582,12 @@
 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
   IssmDouble::Tape& tape = IssmDouble::getTape();
-  codi::DataStore* dataHandler = NULL;
+  codi::ExternalFunctionUserData* dataHandler = NULL;
 
   if(tape.isActive()) {
-    dataHandler = new codi::DataStore();
+    dataHandler = new codi::ExternalFunctionUserData();
 
     // create the index and double vector for the matrix
     IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(local_nnz);
-    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(local_nnz);
+    IssmDouble::Identifier* indexATrans = xNew<IssmDouble::Identifier>(local_nnz);
     int* irnATrans = xNew<int>(local_nnz);
     int* jcnATrans = xNew<int>(local_nnz);
@@ -601,5 +601,5 @@
 
     // create the index vector for a (primal values are not needed for a)
-    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexB = xNew<IssmDouble::Identifier>(n);
     getVectorGradData(rhs, indexB, n);
 
@@ -626,5 +626,5 @@
   if(tape.isActive()) {
     // create the index vector X and register x as active variables
-    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    IssmDouble::Identifier* indexX = xNew<IssmDouble::Identifier>(n);
     registerVector(rhs, indexX, n);
 
@@ -638,6 +638,5 @@
     dataHandler->addData(parameters); // we assume here that parameters is still intact when the reverse run is called
 
-	 //tape.pushExternalFunction(&MumpsSolve_codi_b, dataHandler, &MumpsSolve_codi_delete);
-    tape.pushExternalFunctionHandle(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete);
+	 tape.pushExternalFunction(codi::ExternalFunction<Tape>::create(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete));
   }
   else{