Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 23266)
+++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 23267)
@@ -11,5 +11,10 @@
 #include "../solutionsequences/solutionsequences.h"
 
-#if defined (_HAVE_M1QN3_) && defined(_HAVE_ADOLC_)
+#ifdef _HAVE_CODIPACK_
+extern CoDi_global codi_global;
+#include <sstream> // for output of the CoDiPack tape
+#endif
+
+#if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
 /*m1qn3 prototypes*/
 extern "C" void *ctonbe_; // DIS mode : Conversion
@@ -35,4 +40,5 @@
 void simul_starttrace(FemModel* femmodel){/*{{{*/
 
+	#if defined(_HAVE_ADOLC_)
 	/*Retrive ADOLC parameters*/
 	IssmDouble gcTriggerRatio;
@@ -57,117 +63,46 @@
 	int my_rank=IssmComm::GetRank();
 	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;
+		femmodel->parameters->FindParam(&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
+
+	#else
+	_error_("not implemented");
+	#endif
 }/*}}}*/
-void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
-
-	/*Get rank*/
-	int my_rank=IssmComm::GetRank();
-
-	/*Recover Arguments*/
-	m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
-
-	FemModel* femmodel = input_struct->femmodel;
-	int num_responses,num_controls,numberofvertices,solution_type;
-	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	int* N = NULL;
-	int N_add = 0;
-	int* control_enum = NULL;
-
-	if (solution_type == TransientSolutionEnum){
-		femmodel = input_struct->femmodel->copy();
-	}
-
-	IssmPDouble*  Jlist        = input_struct->Jlist;
-	int           JlistM       = input_struct->M;
-	int           JlistN       = input_struct->N;
-	int*          Jlisti       = input_struct->i;
-	int           intn         = (int)*n;
-
-	/*Recover some parameters*/
-	IssmDouble* scaling_factors = NULL;
-	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
-	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
-	femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
-	femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
-	numberofvertices=femmodel->vertices->NumberOfVertices();
-
-	/*Constrain input vector and update controls*/
-	double  *XL = NULL;
-	double  *XU = NULL;
-	GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
-	GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
-
-	N_add = 0;
-	for (int c=0;c<num_controls;c++){
-		for(int i=0;i<numberofvertices*N[c];i++){
-			int index = N_add*numberofvertices+i;
-			X[index] = X[index]*reCast<double>(scaling_factors[c]);
-			if(X[index]>XU[index]) X[index]=XU[index];
-			if(X[index]<XL[index]) X[index]=XL[index];
-		}
-		N_add+=N[c];
-	}
-
-	/*Start Tracing*/
-	simul_starttrace(femmodel);
-	/*Set X as our new control input and as INDEPENDENT*/
-#ifdef _HAVE_AD_
-	IssmDouble* aX=xNew<IssmDouble>(intn,"t");
-#else
-	IssmDouble* aX=xNew<IssmDouble>(intn);
-#endif
-	if(my_rank==0){
-		for(int i=0;i<intn;i++){
-			aX[i]<<=X[i];
-		}
-	}
-
-	ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	SetControlInputsFromVectorx(femmodel,aX);
-	xDelete<IssmDouble>(aX);
-
-	/*Compute solution (forward)*/
-	void (*solutioncore)(FemModel*)=NULL;
-	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
-	solutioncore(femmodel);
-
-	/*Reset the time to zero for next optimization*/
-	if(solution_type==TransientSolutionEnum){
-		IssmDouble restart_time;
-		femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
-		femmodel->parameters->SetParam(restart_time,TimeEnum);
-
-	}
-
-	/*Get Dependents*/
-	IssmDouble  output_value;
-	int         num_dependents;
-	IssmPDouble *dependents;
-	DataSet*    dependent_objects=NULL;
-	IssmDouble	J=0.;
-	femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-	femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
-
-	/*Go through our dependent variables, and compute the response:*/
-	dependents=xNew<IssmPDouble>(num_dependents);
-	for(int i=0;i<dependent_objects->Size();i++){
-		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-		if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
-		if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
-		if (my_rank==0) {
-			output_value>>=dependents[i];
-			J+=output_value;
-		}
-	}
-
-	/*Turning off trace tape*/
+void simul_stoptrace(){/*{{{*/
+
+	#if defined(_HAVE_ADOLC_)
 	trace_off();
-	//time_t now = time(NULL);
-	//if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
-
-	/*Print tape statistics so that user can kill this run if something is off already:*/
 	if(VerboseAutodiff()){ /*{{{*/
 
 		#ifdef _HAVE_ADOLC_
+		int my_rank=IssmComm::GetRank();
 		size_t  tape_stats[15];
 		tapestats(my_rank,tape_stats); //reading of tape statistics
@@ -222,4 +157,151 @@
 	} /*}}}*/
 
+	#elif defined(_HAVE_CODIPACK_)
+	auto& tape_codi = IssmDouble::getGlobalTape();
+	tape_codi.setPassive();
+	if(VerboseAutodiff()){
+		int my_rank=IssmComm::GetRank();
+		if(my_rank == 0) {
+			// FIXME codi "just because" for now
+			tape_codi.printStatistics(std::cout);
+			codi_global.print(std::cout);
+		}
+	}
+	#else
+	_error_("not implemented");
+	#endif
+}/*}}}*/
+void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
+
+	/*Get rank*/
+	int my_rank=IssmComm::GetRank();
+
+	/*Recover Arguments*/
+	m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
+
+	FemModel* femmodel = input_struct->femmodel;
+	int num_responses,num_controls,numberofvertices,solution_type;
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	int* N = NULL;
+	int N_add = 0;
+	int* control_enum = NULL;
+
+	/*In transient, we need to make sure we do not modify femmodel at each iteration, make a copy*/
+	if(solution_type == TransientSolutionEnum) femmodel = input_struct->femmodel->copy();
+
+	IssmPDouble*  Jlist        = input_struct->Jlist;
+	int           JlistM       = input_struct->M;
+	int           JlistN       = input_struct->N;
+	int*          Jlisti       = input_struct->i;
+	int           intn         = (int)*n;
+
+	/*Recover some parameters*/
+	IssmDouble* scaling_factors = NULL;
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
+	femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
+	femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
+	numberofvertices=femmodel->vertices->NumberOfVertices();
+
+	/*Constrain input vector and update controls*/
+	double  *XL = NULL;
+	double  *XU = NULL;
+	GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+
+	N_add = 0;
+	for (int c=0;c<num_controls;c++){
+		for(int i=0;i<numberofvertices*N[c];i++){
+			int index = N_add*numberofvertices+i;
+			X[index] = X[index]*reCast<double>(scaling_factors[c]);
+			if(X[index]>XU[index]) X[index]=XU[index];
+			if(X[index]<XL[index]) X[index]=XL[index];
+		}
+		N_add+=N[c];
+	}
+
+	/*Start Tracing*/
+	simul_starttrace(femmodel);
+	/*Set X as our new control input and as INDEPENDENT*/
+#ifdef _HAVE_AD_
+	IssmDouble* aX=xNew<IssmDouble>(intn,"t");
+#else
+	IssmDouble* aX=xNew<IssmDouble>(intn);
+#endif
+
+	#if defined(_HAVE_ADOLC_)
+	if(my_rank==0){
+		for(int i=0;i<intn;i++){
+			aX[i]<<=X[i];
+		}
+	}
+	#elif defined(_HAVE_CODIPACK_)
+	auto& tape_codi = IssmDouble::getGlobalTape();
+	if(my_rank==0){
+		for (int i=0;i<intn;i++) {
+			aX[i]=X[i];
+			tape_codi.registerInput(aX[i]);
+			codi_global.input_indices.push_back(aX[i].getGradientData());
+		}
+	}
+	#else
+	_error_("not suppoted");
+	#endif
+
+	ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	SetControlInputsFromVectorx(femmodel,aX);
+	xDelete<IssmDouble>(aX);
+
+	/*Compute solution (forward)*/
+	void (*solutioncore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
+
+	/*Reset the time to zero for next optimization*/
+	if(solution_type==TransientSolutionEnum){
+		IssmDouble restart_time;
+		femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
+		femmodel->parameters->SetParam(restart_time,TimeEnum);
+
+	}
+
+	/*Get Dependents*/
+	IssmDouble  output_value;
+	int         num_dependents;
+	IssmPDouble *dependents;
+	DataSet*    dependent_objects=NULL;
+	IssmDouble	J=0.;
+	femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+	femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
+
+	/*Go through our dependent variables, and compute the response:*/
+	dependents=xNew<IssmPDouble>(num_dependents);
+	for(int i=0;i<dependent_objects->Size();i++){
+		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+		if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
+		if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
+		if(my_rank==0) {
+
+			#if defined(_HAVE_CODIPACK_)
+			tape_codi.registerOutput(output_value);
+			dependents[i] = output_value.getValue();
+			codi_global.output_indices.push_back(output_value.getGradientData());
+
+			#elif defined(_HAVE_ADOLC_)
+			output_value>>=dependents[i];
+
+			#else
+			_error_("not suppoted");
+			#endif
+			J+=output_value;
+		}
+	}
+
+	/*Turning off trace tape*/
+	simul_stoptrace();
+	//time_t now = time(NULL);
+	//if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
+
 	/*diverse: */
 	int  dummy;
@@ -246,4 +328,7 @@
 		num_independents = 0;
 	}
+
+	#if defined(_HAVE_ADOLC_)
+	/*Get gradient for ADOLC {{{*/
 
 	/*get the EDF pointer:*/
@@ -290,4 +375,39 @@
 		xDelete(aWeightVector);
 	}
+	/*}}}*/
+	#elif defined(_HAVE_CODIPACK_)
+	/*Get gradient for CoDiPack{{{*/
+	if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
+	int     aDepIndex=0; /*FIXME: do we really need this?*/
+
+	/*retrieve direction index: */
+	femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
+	if (my_rank==0) {
+		if (aDepIndex<0 || aDepIndex>=num_dependents || codi_global.output_indices.size() <= aDepIndex){
+			_error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
+		}
+		tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0);
+	}
+	tape_codi.evaluate();
+
+	weightVectorTimesJac=xNew<double>(num_independents);
+	/*call driver: */
+	auto in_size = codi_global.input_indices.size();
+	for(size_t i = 0; i < in_size; ++i) {
+		weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
+	}
+
+	/*Add to totalgradient: */
+	totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old);
+	if(my_rank==0) for(int i=0;i<num_independents;i++) {
+		totalgradient[i]+=weightVectorTimesJac[i];
+	}
+
+	/*free resources :*/
+	xDelete(weightVectorTimesJac);
+	/*}}}*/
+	#else
+	_error_("not suppoted");
+	#endif
 
 	/*Broadcast gradient to other ranks*/
@@ -551,4 +671,4 @@
 
 #else
-void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC not installed");}
+void controladm1qn3_core(FemModel* femmodel){_error_("M1QN3 or ADOLC/CoDiPack not installed");}
 #endif //_HAVE_M1QN3_
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp	(revision 23266)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp	(revision 23267)
@@ -103,6 +103,6 @@
 
 	if(isautodiff){
-#if _HAVE_ADOLC_
-		/*Copy some parameters from IoModel to parameters dataset: {{{*/
+		#if defined(_HAVE_ADOLC_)
+		/*Copy some parameters from IoModel to parameters dataset*/
 		parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
@@ -111,6 +111,12 @@
 		parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
 		parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
-		/*}}}*/
-#endif
+
+		#elif defined(_HAVE_CODIPACK_)
+		parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum));
+
+		#else
+		_error_("not supported yet");
+		#endif
+
 		/*retrieve driver: {{{*/
 		iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");