Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23250)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23251)
@@ -17,4 +17,9 @@
 #include "../shared/Enum/Enum.h"
 #include "../analyses/analyses.h"
+
+#if _HAVE_CODIPACK_
+#include <sstream> // for output of the CoDiPack tape
+extern CoDi_global codi_global;
+#endif
 
 /*module includes: {{{*/
@@ -1851,4 +1856,8 @@
 			dependents=xNew<IssmPDouble>(num_dependents);
 
+			#if defined(_HAVE_CODIPACK_)
+			auto& tape_codi = IssmDouble::getGlobalTape();
+			#endif
+
 			/*Go through our dependent variables, and compute the response:*/
 			for(int i=0;i<dependent_objects->Size();i++){
@@ -1856,5 +1865,11 @@
 				dep->Responsex(&output_value,this);
 				if (my_rank==0) {
-					output_value>>=dependents[i];
+					#if defined(_HAVE_CODIPACK_)
+						tape_codi.registerOutput(output_value);
+						dependents[i] = output_value.getValue();
+						codi_global.output_indices.push_back(output_value.getGradientData());
+					#else
+						output_value>>=dependents[i];
+					#endif
 				}
 			}
Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23250)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23251)
@@ -27,4 +27,9 @@
 #include "../shared/io/io.h"
 #include "../shared/shared.h"
+
+#ifdef _HAVE_CODIPACK_
+extern CoDi_global codi_global;
+#include <sstream> // for output of the CoDiPack tape
+#endif
 
 /*IoConstant class and methods*/
@@ -431,5 +436,6 @@
 	if(trace || (autodiff && !iscontrol)){
 
-		#ifdef _HAVE_ADOLC_
+		#ifdef _HAVE_AD_
+		// FIXME codi here we should be able to execute codi version as normal
 		this->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
 		if(num_independent_objects){
@@ -460,5 +466,5 @@
 		#else
 		/*if we asked for AD computations, we have a problem!: */
-		_error_("Cannot carry out AD mode computations without support of ADOLC compiled in!");
+		_error_("Cannot carry out AD mode computations without support of ADOLC or CoDiPack compiled in!");
 		#endif
 	}
@@ -1702,5 +1708,5 @@
 	if(X) Xcount=*pXcount;
 
-	#ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
+	#ifdef _HAVE_AD_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
 
 	/*output: */
@@ -1719,10 +1725,22 @@
 		/*Now, before we even broadcast this to other nodes, declare the scalar  as an independent variable!. If we 
 		 *have been supplied an X vector, use it instead of what we just read: */
-		if(X){
-			scalar<<=X[Xcount];
-		}
-		else{
-			scalar<<=pscalar;
-		}
+		#if defined(_HAVE_CODIPACK_)
+			// FIXME codi here we just assign instead of using "operator <<="
+			if(X){
+				scalar=X[Xcount];
+			} else {
+				scalar=pscalar;
+			}
+			auto& tape_codi = IssmDouble::getGlobalTape();
+			tape_codi.registerInput(scalar);
+			codi_global.input_indices.push_back(scalar.getGradientData());
+		#else
+			if(X){
+				scalar<<=X[Xcount];
+			}
+			else{
+				scalar<<=pscalar;
+			}
+		#endif
 	}
 
@@ -1785,10 +1803,29 @@
 			/*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!
 			  If we have been supplied an X vector, use it instead of what we just read: */
-			if(X){
-				for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
-			}
-			else{
-				for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
-			}
+			#if defined(_HAVE_CODIPACK_)
+				// FIXME codi here we just assign instead of using "operator <<="
+				auto& tape_codi = IssmDouble::getGlobalTape();
+				if(X){
+					for (int i=0;i<M*N;i++) {
+						matrix[i]=X[Xcount+i];
+						tape_codi.registerInput(matrix[i]);
+						codi_global.input_indices.push_back(matrix[i].getGradientData());
+					}
+				}
+				else{
+					for (int i=0;i<M*N;i++) {
+						matrix[i]=buffer[i];
+						tape_codi.registerInput(matrix[i]);
+						codi_global.input_indices.push_back(matrix[i].getGradientData());
+					}
+				}
+			#else /*ADOLC*/
+				if(X){
+					for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
+				}
+				else{
+					for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
+				}
+			#endif
 		}
 		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
@@ -1798,4 +1835,6 @@
 	else _error_("cannot declare the independent variable \"" << data_name <<  "\" if it's empty!");
 
+	// FIXME codi is that at all relevant to CoDiPack or can we simply assume the same?
+	
 	/*Add to data as independent*/
 	this->AddDataIndependent(new IoData(matrix,code,layout,M,N,data_name));
@@ -2562,5 +2601,6 @@
 
 	if(trace || (autodiff && !iscontrol)){
-		#ifdef _HAVE_ADOLC_
+
+		#if defined(_HAVE_ADOLC_)
 		/*Retrieve parameters: */
 		this->FetchData(&keep,"md.autodiff.keep");
@@ -2579,7 +2619,35 @@
 		int skipFileDeletion=1;
 		trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
+
+		#elif defined(_HAVE_CODIPACK_)
+		//fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
+		/*
+		 * FIXME codi
+		 * - ADOL-C variant uses fine grained tracing with various arguments
+		 * - ADOL-C variant sets a garbage collection parameter for its tape
+		 * -> These parameters are not read for the CoDiPack ISSM version!
+		 */
+		auto& tape_codi = IssmDouble::getGlobalTape();
+		tape_codi.setActive();
+		#if _AD_TAPE_ALLOC_
+		//alloc_profiler.Tag(StartInit, true);
+		IssmDouble x_t(1.0), y_t(1.0);
+		tape_codi.registerInput(y_t);
+		int codi_allocn = 0;
+		this->FetchData(&codi_allocn, AutodiffTapeAllocEnum);
+		for(int i = 0;i < codi_allocn;++i) {
+			x_t = y_t * y_t;
+		}
+		/*
+		std::stringstream out_s;
+		IssmDouble::getGlobalTape().printStatistics(out_s);
+		_printf0_("StartTrace::Tape Statistics	   : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
+		*/
+		tape_codi.reset();
+		//alloc_profiler.Tag(FinishInit, true);
 		#endif
-	}
-
-}
-/*}}}*/
+		#endif
+	}
+
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 23250)
+++ /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 23251)
@@ -14,4 +14,5 @@
 #include "../../classes/Params/Parameters.h"
 #include "../adolc/adolcincludes.h"
+#include "../codipack/codipackincludes.h"
 #include "./gslincludes.h"
 
@@ -251,2 +252,141 @@
 /*}}}*/
 #endif
+
+#ifdef _HAVE_CODIPACK_
+void SolverxSeq_codi_b(codi::DataStore* data) {/*{{{*/
+  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
+
+  IssmDouble::Real* valueATrans;
+  IssmDouble::GradientData* indexATrans;
+  IssmDouble::GradientData* indexB;
+  IssmDouble::Real* valueX;
+  IssmDouble::GradientData* indexX;
+  int n;
+
+  data->getData(valueATrans);
+  data->getData(indexATrans);
+  data->getData(indexB);
+  data->getData(valueX);
+  data->getData(indexX);
+  data->getData(n);
+
+
+  // create the adjoint vector for x and reset the adjoint values on the tape
+  IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
+  getVectorAdjoint(tape, indexX, adjX, n);
+
+  IssmDouble::GradientValue* sol  = xNew<IssmDouble::GradientValue>(n);
+  SolverxSeq(sol, valueATrans, adjX, 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]);
+    }
+  }
+
+  xDelete(sol);
+  xDelete(adjX);
+}
+/*}}}*/
+void SolverxSeq_codi_delete(codi::DataStore* data) {/*{{{*/
+  IssmDouble::Real* valueATrans;
+  IssmDouble::GradientData* indexATrans;
+  IssmDouble::GradientData* indexB;
+  IssmDouble::Real* valueX;
+  IssmDouble::GradientData* indexX;
+  int n;
+
+  data->getData(valueATrans);
+  data->getData(indexATrans);
+  data->getData(indexB);
+  data->getData(valueX);
+  data->getData(indexX);
+  data->getData(n);
+
+  xDelete(valueATrans);
+  xDelete(indexATrans);
+  xDelete(indexB);
+  xDelete(valueX);
+  xDelete(indexX);
+}
+/*}}}*/
+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
+  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
+  codi::DataStore* dataHandler = NULL;
+
+  if(tape.isActive()) {
+    dataHandler = new codi::DataStore();
+
+    // 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);
+
+    // read the data for matrix in a transposed fashion
+	  for (int i=0; i<n; ++i) {
+      for (int j=0; j<n; ++j) {
+        getPrimalAndGradData(A[i*n+j], valueATrans[j*n+i], indexATrans[j*n+i]);
+      }
+    }
+
+    // read the data from B (primal values are not required vor B
+    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    getVectorGradData(B, indexB, n);
+
+    dataHandler->addData(valueATrans);
+    dataHandler->addData(indexATrans);
+    dataHandler->addData(indexB);
+  }
+
+  // unpack the primal values from the matrix and the vector
+  IssmDouble::Real* valueA = xNew<IssmDouble::Real>(n*n);
+  IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
+  // read the data from A and B
+  getVectorPrimal(A, valueA, n*n);
+  getVectorPrimal(B, valueB, n);
+
+  // create the placeholder for X and solve the system
+  IssmDouble::Real* valueX = xNew<IssmDouble::Real>(n);
+  SolverxSeq(valueX, valueA, valueB, n);
+
+  // pack the values into x
+  setVectorPrimal(X, valueX, n);
+
+  if(tape.isActive()) {
+    // create the index vector X and register x as active variables
+    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    registerVector(X, indexX, n);
+
+    dataHandler->addData(valueX);
+    dataHandler->addData(indexX);
+
+    // store other arguments
+    dataHandler->addData(n);
+
+    tape.pushExternalFunction(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete);
+  } else {
+    // if the tape is active valueX is stored in the dataHandler and deleted in the reverse sweep
+    xDelete(valueX);
+  }
+
+  xDelete(valueB);
+  xDelete(valueA);
+}
+/*}}}*/
+void DenseGslSolve(/*output*/ IssmDouble** px,/*stiffness matrix:*/ IssmDouble* Kff, int Kff_M, int Kff_N, /*right hand side load vector: */ IssmDouble* pf, int pf_M, Parameters* parameters){ /*{{{*/
+
+	/*Intermediary: */
+
+	if(Kff_N!=pf_M)_error_("Right hand side vector of size " << pf_M << ", when matrix is of size " << Kff_M << "-" << Kff_N << " !");
+	if(Kff_M!=Kff_N)_error_("Stiffness matrix should be square!");
+
+	IssmDouble *x  = xNew<IssmDouble>(Kff_N,"t");
+
+	SolverxSeq(x,Kff,pf,Kff_N,parameters);
+
+	/*allocate output pointers: */
+	*px=x;
+}
+/*}}}*/
+#endif
Index: /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.cpp	(revision 23250)
+++ /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.cpp	(revision 23251)
@@ -60,5 +60,5 @@
 # endif
 #else
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -97,5 +97,5 @@
 #else
   assert(sendcount==recvcounts[0]); // we handle only identical representations 
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -129,5 +129,5 @@
 # endif
 #else
-# ifdef _HAVE_ADOLC_
+#ifdef _HAVE_AD_
   if (datatype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -215,12 +215,16 @@
 int ISSM_MPI_Finalize(void){  /*{{{*/
 
-  int rc=0;
-#ifdef _HAVE_MPI_
-# ifdef _HAVE_AMPI_
-  rc=AMPI_Finalize_NT();
-# else
-  rc=MPI_Finalize();
-# endif
-#endif
+	int rc=0;
+	#ifdef _HAVE_MPI_
+		#ifdef _HAVE_AMPI_
+			#ifdef _HAVE_ADJOINTMPI_
+				rc=AMPI_Finalize();
+			#else
+				rc=AMPI_Finalize_NT();
+			#endif
+		#else
+		  rc=MPI_Finalize();
+		#endif
+	#endif
   return rc;
 }/*}}}*/
@@ -250,5 +254,5 @@
 # endif
 #else
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -290,5 +294,5 @@
 #else
   assert(sendcnt==recvcnts[0]); // we handle only identical representations 
-# ifdef _HAVE_ADOLC_
+#ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -304,12 +308,16 @@
 int ISSM_MPI_Init(int *argc, char ***argv){  /*{{{*/
 
-  int rc=0;
-#ifdef _HAVE_MPI_
-# ifdef _HAVE_AMPI_
-  rc=AMPI_Init_NT(argc,argv);
-# else
-  rc=MPI_Init(argc,argv);
-# endif
-#endif
+	int rc=0;
+	#ifdef _HAVE_MPI_
+		#ifdef _HAVE_AMPI_
+			#ifdef _HAVE_ADJOINTMPI_
+				rc=AMPI_Init(argc,argv);
+			#else
+				rc=AMPI_Init_NT(argc,argv);
+			#endif
+		#else
+			rc=MPI_Init(argc,argv);
+		#endif
+	#endif
   return rc;
 }/*}}}*/
@@ -324,5 +332,7 @@
 	       source,
 	       tag,
+			 #ifndef _HAVE_ADJOINTMPI_
 	       AMPI_FROM_SEND, // as long as there are no other variants
+			 #endif
 	       comm,
 	       status);
@@ -365,5 +375,5 @@
 # endif
 #else
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (datatype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -402,5 +412,5 @@
 # endif
 #else
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
@@ -442,5 +452,5 @@
 #else
   assert(sendcnts[0]==recvcnt); // we handle only identical representations 
-# ifdef _HAVE_ADOLC_
+# ifdef _HAVE_AD_
   if (sendtype==ISSM_MPI_DOUBLE) { 
     IssmDouble* activeSendBuf=(IssmDouble*)(sendbuf)+displs[0];
@@ -464,5 +474,7 @@
 	       dest,
 	       tag,
+			 #ifndef _HAVE_ADJOINTMPI_
 	       AMPI_TO_RECV, // as long as there are no other variants
+			 #endif
 	       comm);
 # else
@@ -492,4 +504,6 @@
 #ifdef _HAVE_ADOLC_
   ensureContiguousLocations(aSize);
+#else
+  fprintf(stderr, "*** Codipack ISSM_MPI_ContiguousInAdolc()\n");
 #endif
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h	(revision 23250)
+++ /issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h	(revision 23251)
@@ -19,5 +19,9 @@
 	/*Include header files: {{{*/
 	#ifdef _HAVE_AMPI_
-		#include <ampi/ampi.h>
+		#ifdef _HAVE_ADJOINTMPI_
+			#include <ampi_tape.hpp>
+		#else
+			#include <ampi/ampi.h>
+		#endif
 	#elif  _HAVE_PETSC_MPI_ // Petsc now hides there MPI header. It can be reached through Petsc.
 		#include <petsc.h>
@@ -35,6 +39,8 @@
 	// data types
 	#define ISSM_MPI_CHAR       MPI_CHAR
-	#ifdef _HAVE_AMPI_
+	#if defined(_HAVE_AMPI_)
 		#define ISSM_MPI_DOUBLE    AMPI_ADOUBLE // corresponds to IssmDouble
+	#elif defined(_HAVE_ADJOINTMPI_)
+		#define ISSM_MPI_DOUBLE    AMPI_DOUBLE // corresponds to IssmDouble
 	#else 
 		#define ISSM_MPI_DOUBLE    MPI_DOUBLE  // corresponds to IssmDouble 
Index: /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 23250)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 23251)
@@ -18,4 +18,5 @@
 #include "../mpi/issmmpi.h"
 #include "../adolc/adolcincludes.h"
+#include "../codipack/codipackincludes.h"
 #include "../issm/SparseRow.h"
 #include "./mumpsincludes.h"
@@ -98,5 +99,5 @@
 }
 
-#ifdef _HAVE_ADOLC_
+#ifdef _HAVE_AD_
 // prototype for active variant
 void MumpsSolve(int n,
@@ -354,6 +355,5 @@
 
 #ifdef _HAVE_ADOLC_
-
-int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) {
+int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y){/*{{{*/
   // unpack parameters
   int n=iArr[0];
@@ -381,7 +381,6 @@
   xDelete(local_irn);
   return 0;
-}
-
-void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){
+}/*}}}*/
+void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
   int packedDimsSparseArrLength=1+1+1+local_nnz+local_nnz;
   int *packedDimsSparseArr=xNew<int>(packedDimsSparseArrLength);
@@ -425,10 +424,6 @@
   xDelete(pack_A_rhs);
   xDelete(packedDimsSparseArr);
-}
-
-int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, 
-			      int m, IssmPDouble *dp_U,
-			      int nPlusNz, IssmPDouble *dp_Z,
-			      IssmPDouble *dp_x, IssmPDouble *dp_y) {
+}/*}}}*/
+int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, int m, IssmPDouble *dp_U, int nPlusNz, IssmPDouble *dp_Z, IssmPDouble *dp_x, IssmPDouble *dp_y) {/*{{{*/
   // unpack parameters
   int n=iArr[0];
@@ -481,5 +476,152 @@
   xDelete(local_irn);
   return 3;
-}
-
-#endif
+}/*}}}*/
+#endif
+
+#ifdef _HAVE_CODIPACK_
+void MumpsSolve_codi_b(codi::DataStore* data) {/*{{{*/
+  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
+
+  IssmDouble::Real* valueATrans;
+  IssmDouble::GradientData* indexATrans;
+  int* irnATrans;
+  int* jcnATrans;
+  IssmDouble::GradientData* indexB;
+  IssmDouble::Real* valueX;
+  IssmDouble::GradientData* indexX;
+  int n;
+  int nnz;
+  int local_nnz;
+  Parameters* parameters;
+
+  data->getData(valueATrans);
+  data->getData(indexATrans);
+  data->getData(irnATrans);
+  data->getData(jcnATrans);
+  data->getData(indexB);
+  data->getData(valueX);
+  data->getData(indexX);
+  data->getData(n);
+  data->getData(nnz);
+  data->getData(local_nnz);
+  data->getData(parameters);
+
+  // create the adjoint vector for x and reset the adjoint values on the tape
+  IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(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);
+  for(int i=0; i<local_nnz; ++i) {
+    // 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]);
+  }
+
+  xDelete(adjX);
+}
+/*}}}*/
+void MumpsSolve_codi_delete(codi::DataStore* data) {/*{{{*/
+  IssmDouble::Real* valueATrans;
+  IssmDouble::GradientData* indexATrans;
+  int* irnATrans;
+  int* jcnATrans;
+  IssmDouble::GradientData* indexB;
+  IssmDouble::Real* valueX;
+  IssmDouble::GradientData* indexX;
+  int n;
+  int nnz;
+  int local_nnz;
+  Parameters* parameters;
+
+  data->getData(valueATrans);
+  data->getData(indexATrans);
+  data->getData(irnATrans);
+  data->getData(jcnATrans);
+  data->getData(indexB);
+  data->getData(valueX);
+  data->getData(indexX);
+  data->getData(n);
+  data->getData(nnz);
+  data->getData(local_nnz);
+  data->getData(parameters);
+
+  xDelete(valueATrans);
+  xDelete(indexATrans);
+  xDelete(irnATrans);
+  xDelete(jcnATrans);
+  xDelete(indexB);
+  xDelete(valueX);
+  xDelete(indexX);
+}
+/*}}}*/
+void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
+  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
+  codi::DataStore* dataHandler = NULL;
+
+  if(tape.isActive()) {
+    dataHandler = new codi::DataStore();
+
+    // 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);
+    int* irnATrans = xNew<int>(local_nnz);
+    int* jcnATrans = xNew<int>(local_nnz);
+
+    // read the data for the matrix A in a transposed fashion
+	  for (int i=0; i<local_nnz; ++i) {
+        getPrimalAndGradData(a_loc[i], valueATrans[i], indexATrans[i]);
+        irnATrans[i]=jcn_loc[i];  // transposed store
+        jcnATrans[i]=irn_loc[i];  // transposed store
+    }
+
+    // create the index vector for a (primal values are not needed for a)
+    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
+    getVectorGradData(rhs, indexB, n);
+
+    dataHandler->addData(valueATrans);
+    dataHandler->addData(indexATrans);
+    dataHandler->addData(irnATrans);
+    dataHandler->addData(jcnATrans);
+    dataHandler->addData(indexB);
+  }
+
+  // unpack the primal values from the matrix and the vector
+  IssmDouble::Real* valueA = xNew<IssmDouble::Real>(local_nnz);
+  IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
+  // read the data from A and B
+  getVectorPrimal(a_loc, valueA, local_nnz);
+  getVectorPrimal(rhs, valueB, n);
+
+  MumpsSolve(n, nnz, local_nnz, irn_loc, jcn_loc, valueA, valueB, parameters);
+  // valueB contains now the solution
+
+  // pack the values into rhs
+  setVectorPrimal(rhs, valueB, n);
+
+  if(tape.isActive()) {
+    // create the index vector X and register x as active variables
+    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
+    registerVector(rhs, indexX, n);
+
+    dataHandler->addData(valueB); // contains the values from x
+    dataHandler->addData(indexX);
+
+    // store other arguments
+    dataHandler->addData(n);
+    dataHandler->addData(nnz);
+    dataHandler->addData(local_nnz);
+    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);
+  } else {
+    // if the tape is active valueB is stored in the dataHandler and deleted in the reverse sweep
+    xDelete(valueB);
+  }
+
+  xDelete(valueA);
+}
+/*}}}*/
+#endif