Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 26888)
+++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 26889)
@@ -15,5 +15,5 @@
 #include <sstream> // for output of the CoDiPack tape
 #include <fenv.h>
-void transient_ad(FemModel* femmodel);
+double transient_ad(FemModel* femmodel, double* G,double* Jlist);
 #endif
 
@@ -222,4 +222,5 @@
 	int    *N = NULL;
 	int    *control_enum    = NULL;
+	int     checkpoint_frequency;
 	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
@@ -228,4 +229,5 @@
 	femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
 	femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
+	femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum);
 
 	/*Constrain input vector and update controls*/
@@ -236,5 +238,5 @@
 
 	int offset = 0;
-	for (int c=0;c<num_controls;c++){
+	for(int c=0;c<num_controls;c++){
 		for(int i=0;i<M[c]*N[c];i++){
 			int index = offset+i;
@@ -246,78 +248,88 @@
 	}
 
-	/*Start Tracing*/
-	simul_starttrace(femmodel);
-	/*Set X as our new control input and as INDEPENDENT*/
+	/*Special case: do we need to run AD with checkpointing?*/
+	#ifdef _HAVE_CODIPACK_
+	if(checkpoint_frequency && solution_type == TransientSolutionEnum){
+		SetControlInputsFromVectorx(femmodel,X);
+		*pf = transient_ad(femmodel, G, &Jlist[(*Jlisti)*JlistN]);
+	}
+	else
+	#endif
+	  {
+
+		/*Start Tracing*/
+		simul_starttrace(femmodel);
+		/*Set X as our new control input and as INDEPENDENT*/
 #ifdef _HAVE_AD_
-	IssmDouble* aX=xNew<IssmDouble>(intn,"t");
+		IssmDouble* aX=xNew<IssmDouble>(intn,"t");
 #else
-	IssmDouble* aX=xNew<IssmDouble>(intn);
+		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_)
-
-	/*Get tape*/
-	#if _CODIPACK_MAJOR_==2
-	auto& tape_codi = IssmDouble::getTape();
-	#elif _CODIPACK_MAJOR_==1
-	auto& tape_codi = IssmDouble::getGlobalTape();
-	#else
-	#error "_CODIPACK_MAJOR_ not supported"
-	#endif
-
-	codi_global.input_indices.clear();
-	if(my_rank==0){
-		for (int i=0;i<intn;i++) {
-			aX[i]=X[i];
-			tape_codi.registerInput(aX[i]);
-			#if _CODIPACK_MAJOR_==2
-			codi_global.input_indices.push_back(aX[i].getIdentifier());
-			#elif _CODIPACK_MAJOR_==1
-			codi_global.input_indices.push_back(aX[i].getGradientData());
-			#else
-			#error "_CODIPACK_MAJOR_ not supported"
-			#endif
-
-		}
-	}
-	#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);
-
-	/*Get Dependents*/
-	IssmDouble   output_value;
-	int          num_dependents;
-	IssmPDouble *dependents;
-	IssmDouble   J = 0.;
-	DataSet     *dependent_objects = ((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
-	femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
-
-	/*Go through our dependent variables, and compute the response:*/
-	dependents=xNew<IssmPDouble>(num_dependents);
-	#if defined(_HAVE_CODIPACK_)
-	codi_global.output_indices.clear();
-	#endif
-	int i=-1;
-	for(Object* & object:dependent_objects->objects){
-		i++;
-		DependentObject* dep=xDynamicCast<DependentObject*>(object);
-		if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
-		if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
-		if(my_rank==0) {
+		#if defined(_HAVE_ADOLC_)
+		if(my_rank==0){
+			for(int i=0;i<intn;i++){
+				aX[i]<<=X[i];
+			}
+		}
+		#elif defined(_HAVE_CODIPACK_)
+
+		/*Get tape*/
+		#if _CODIPACK_MAJOR_==2
+		auto& tape_codi = IssmDouble::getTape();
+		#elif _CODIPACK_MAJOR_==1
+		auto& tape_codi = IssmDouble::getGlobalTape();
+		#else
+		#error "_CODIPACK_MAJOR_ not supported"
+		#endif
+
+		codi_global.input_indices.clear();
+		if(my_rank==0){
+			for (int i=0;i<intn;i++) {
+				aX[i]=X[i];
+				tape_codi.registerInput(aX[i]);
+				#if _CODIPACK_MAJOR_==2
+				codi_global.input_indices.push_back(aX[i].getIdentifier());
+				#elif _CODIPACK_MAJOR_==1
+				codi_global.input_indices.push_back(aX[i].getGradientData());
+				#else
+				#error "_CODIPACK_MAJOR_ not supported"
+				#endif
+
+			}
+		}
+		#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);
+
+		/*Get Dependents*/
+		IssmDouble   output_value;
+		int          num_dependents;
+		IssmPDouble *dependents;
+		IssmDouble   J = 0.;
+		DataSet     *dependent_objects = ((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
+		femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+
+		/*Go through our dependent variables, and compute the response:*/
+		dependents=xNew<IssmPDouble>(num_dependents);
+		#if defined(_HAVE_CODIPACK_)
+		codi_global.output_indices.clear();
+		#endif
+		int i=-1;
+		IssmDouble TEMP = 0.;
+		for(Object* & object:dependent_objects->objects){
+			i++;
+			DependentObject* dep=xDynamicCast<DependentObject*>(object);
+			if(solution_type==TransientSolutionEnum) output_value = dep->GetValue();
+			if(solution_type!=TransientSolutionEnum) dep->Responsex(&output_value,femmodel);
 
 			#if defined(_HAVE_CODIPACK_)
@@ -340,141 +352,145 @@
 			J+=output_value;
 		}
-	}
-
-	/*Turning off trace tape*/
-	simul_stoptrace();
-
-	/*intermediary: */
-	int          num_independents=intn;
-	IssmPDouble *aWeightVector=NULL;
-	IssmPDouble *weightVectorTimesJac=NULL;
-	IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
-
-	/*if no dependents, no point in running a driver: */
-	if(!(num_dependents*num_independents)) _error_("this is not allowed");
-
-	/*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
-	int num_dependents_old   = num_dependents;
-	int num_independents_old = num_independents;
-
-	#if defined(_HAVE_ADOLC_)
-	/*Get gradient for ADOLC {{{*/
-	if(my_rank!=0){
-		num_dependents   = 0;
-		num_independents = 0;
-	}
-
-	/*get the EDF pointer:*/
-	ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
-
-	/* these are always needed regardless of the interpreter */
-	anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
-	anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
-
-	/* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
-	 * as to generate num_dependents gradients: */
-	for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
-
-		/*initialize direction index in the weights vector: */
-		aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
-		if (my_rank==0) aWeightVector[aDepIndex]=1.;
-
-		/*initialize output gradient: */
-		weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
-
-		/*set the forward method function pointer: */
-		#ifdef _HAVE_GSL_
-		anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
-		#endif
-		#ifdef _HAVE_MUMPS_
-		anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
-		#endif
-
-		anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
-		anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
-
-		/*call driver: */
-		fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
-
-		/*Add to totalgradient: */
-		if(my_rank==0) for(int i=0;i<num_independents;i++) {
-			totalgradient[i]+=weightVectorTimesJac[i];
-		}
-
-		/*free resources :*/
-		xDelete(weightVectorTimesJac);
-		xDelete(aWeightVector);
-	}
-	/*}}}*/
-	#elif defined(_HAVE_CODIPACK_)
-	/*Get gradient for CoDiPack{{{*/
-	if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
-
-	/* call the fos_reverse in a loop on the index, from 0 to num_dependents, so
-	 * as to generate num_dependents gradients: */
-	for(int dep_index=0;dep_index<num_dependents_old;dep_index++){
-
-		/*initialize direction index in the weights vector: */
-		if(my_rank==0){
-			if(dep_index<0 || dep_index>=num_dependents || codi_global.output_indices.size() <= dep_index){
-				_error_("index value for dependent index should be in [0,num_dependents-1]");
+
+		/*Turning off trace tape*/
+		simul_stoptrace();
+
+		/*intermediary: */
+		int          num_independents=intn;
+		IssmPDouble *aWeightVector=NULL;
+		IssmPDouble *weightVectorTimesJac=NULL;
+		IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
+
+		/*if no dependents, no point in running a driver: */
+		if(!(num_dependents*num_independents)) _error_("this is not allowed");
+
+		/*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
+		int num_dependents_old   = num_dependents;
+		int num_independents_old = num_independents;
+
+		#if defined(_HAVE_ADOLC_)
+		/*Get gradient for ADOLC {{{*/
+		if(my_rank!=0){
+			num_dependents   = 0;
+			num_independents = 0;
+		}
+
+		/*get the EDF pointer:*/
+		ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
+
+		/* these are always needed regardless of the interpreter */
+		anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
+		anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
+
+		/* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
+		 * as to generate num_dependents gradients: */
+		for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
+
+			/*initialize direction index in the weights vector: */
+			aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
+			if (my_rank==0) aWeightVector[aDepIndex]=1.;
+
+			/*initialize output gradient: */
+			weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
+
+			/*set the forward method function pointer: */
+			#ifdef _HAVE_GSL_
+			anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
+			#endif
+			#ifdef _HAVE_MUMPS_
+			anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
+			#endif
+
+			anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
+			anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
+
+			/*call driver: */
+			fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
+
+			/*Add to totalgradient: */
+			if(my_rank==0) for(int i=0;i<num_independents;i++) {
+				totalgradient[i]+=weightVectorTimesJac[i];
 			}
-			tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
-		}
-		//feclearexcept(FE_ALL_EXCEPT);
-		//feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-		tape_codi.evaluate();
-
-		/*Get gradient for this dependent */
-		weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
-		auto in_size = codi_global.input_indices.size();
-		for(size_t i = 0; i < in_size; ++i){
-			_assert_(i<num_independents);
-			weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
-		}
-		if(my_rank==0) for(int i=0;i<num_independents;i++){
-			totalgradient[i]+=weightVectorTimesJac[i];
-		}
-		xDelete(weightVectorTimesJac);
-	}
-
-	/*Clear tape*/
-	tape_codi.reset();
-	/*}}}*/
-	#else
-	_error_("not suppoted");
-	#endif
-
-	/*Broadcast gradient to other ranks*/
-	ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-	/*Check size of Jlist to avoid crashes*/
-	_assert_((*Jlisti)<JlistM);
-	_assert_(JlistN==num_responses+1);
-
-	/*Compute objective function and broadcast it to other cpus*/
-	*pf = reCast<double>(J);
-	ISSM_MPI_Bcast(pf,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-	_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
-
-	/*Record cost function values and delete Jtemp*/
-	for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = dependents[i];
-	Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
-
-	if(*indic==0){
-		/*dry run, no gradient required*/
-
-		/*Retrieve objective functions independently*/
-		_printf0_("            N/A |\n");
-		for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
-		_printf0_("\n");
-
-		*Jlisti = (*Jlisti) +1;
-		xDelete<double>(XU);
-		xDelete<double>(XL);
-		return;
-	}
-
-	/*Compute gradient*/
-	for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
+
+			/*free resources :*/
+			xDelete(weightVectorTimesJac);
+			xDelete(aWeightVector);
+		}
+		/*}}}*/
+		#elif defined(_HAVE_CODIPACK_)
+		/*Get gradient for CoDiPack{{{*/
+		if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
+
+		/* call the fos_reverse in a loop on the index, from 0 to num_dependents, so
+		 * as to generate num_dependents gradients: */
+		for(int dep_index=0;dep_index<num_dependents_old;dep_index++){
+
+			/*initialize direction index in the weights vector: */
+			if(my_rank==0){
+				if(dep_index<0 || dep_index>=num_dependents || codi_global.output_indices.size() <= dep_index){
+					_error_("index value for dependent index should be in [0,num_dependents-1]");
+				}
+				tape_codi.setGradient(codi_global.output_indices[dep_index],1.0);
+			}
+			//feclearexcept(FE_ALL_EXCEPT);
+			//feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+			tape_codi.evaluate();
+
+			/*Get gradient for this dependent */
+			weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
+			auto in_size = codi_global.input_indices.size();
+			for(size_t i = 0; i < in_size; ++i){
+				_assert_(i<num_independents);
+				weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
+			}
+			if(my_rank==0) for(int i=0;i<num_independents;i++){
+				totalgradient[i]+=weightVectorTimesJac[i];
+			}
+			xDelete(weightVectorTimesJac);
+		}
+
+		/*Clear tape*/
+		tape_codi.reset();
+		/*}}}*/
+		#else
+		_error_("not suppoted");
+		#endif
+
+		/*Broadcast gradient to other ranks*/
+		ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+		/*Check size of Jlist to avoid crashes*/
+		_assert_((*Jlisti)<JlistM);
+		_assert_(JlistN==num_responses+1);
+
+		/*Compute objective function and broadcast it to other cpus*/
+		*pf = reCast<double>(J);
+		ISSM_MPI_Bcast(pf,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+		/*Record cost function values and delete Jtemp*/
+		for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = dependents[i];
+		Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
+
+		if(*indic==0){
+			/*dry run, no gradient required*/
+
+			/*Retrieve objective functions independently*/
+			_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+			_printf0_("            N/A |\n");
+			for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
+			_printf0_("\n");
+
+			*Jlisti = (*Jlisti) +1;
+			xDelete<double>(XU);
+			xDelete<double>(XL);
+			return;
+		}
+
+		/*Compute gradient*/
+		for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
+
+		xDelete<IssmPDouble>(dependents);
+		xDelete<IssmPDouble>(totalgradient);
+	  } /*====????*/
 
 	/*Constrain Gradient*/
@@ -496,4 +512,5 @@
 
 	/*Print info*/
+	_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
 	_printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
 	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
@@ -509,21 +526,7 @@
 	xDelete<int>(N);
 	xDelete<double>(scaling_factors);
-	xDelete<IssmPDouble>(dependents);
-	xDelete<IssmPDouble>(totalgradient);
 }/*}}}*/
 void controladm1qn3_core(FemModel* femmodel){/*{{{*/
 
-	int checkpoint_frequency;
-	femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum);
-	if(checkpoint_frequency){
-		#ifdef _HAVE_CODIPACK_
-		transient_ad(femmodel);
-		femmodel->OutputControlsx(&femmodel->results);
-		printf("No optimization implemented yet, skipping\n");
-		return;
-		#else
-		_error_("checkpointing not implemented for ADOLC");
-		#endif
-	}
 
 	/*Intermediaries*/
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 26888)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 26889)
@@ -16,4 +16,8 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+#ifdef _HAVE_CODIPACK_
+extern CoDi_global codi_global;
+#endif
 
 /*Prototypes*/
@@ -275,5 +279,5 @@
 
 #ifdef _HAVE_CODIPACK_
-void transient_ad(FemModel* femmodel){/*{{{*/
+double transient_ad(FemModel* femmodel, double* G, double* Jlist){/*{{{*/
 
 	/*parameters: */
@@ -282,5 +286,5 @@
 	bool       isoceancoupling;
 	int        step,timestepping;
-	int        checkpoint_frequency;
+	int        checkpoint_frequency,num_responses;
 
 	/*Get rank*/
@@ -293,4 +297,5 @@
 	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
 	femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 	femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum); _assert_(checkpoint_frequency>0);
 
@@ -343,11 +348,9 @@
 
 		/*Go through our dependent variables, and compute the response:*/
-		if(finaltime==time){
-			DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
-			for(Object* & object:dependent_objects->objects){
-				DependentObject* dep=(DependentObject*)object;
-				dep->Responsex(&output_value,femmodel);
-				dep->AddValue(output_value);
-			}
+		DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
+		for(Object* & object:dependent_objects->objects){
+			DependentObject* dep=(DependentObject*)object;
+			dep->Responsex(&output_value,femmodel);
+			dep->AddValue(output_value);
 		}
 
@@ -395,15 +398,37 @@
 	SetControlInputsFromVectorx(femmodel,X);
 
-	IssmDouble J = 0.;
+	IssmDouble J     = 0.;
+	int        count = 0;
 	DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
 	for(Object* & object:dependent_objects->objects){
 		DependentObject* dep=(DependentObject*)object;
-		J += dep->GetValue();
-	}
-	tape_codi.registerOutput(J);
+		IssmDouble       output_value = dep->GetValue();
+
+
+		J += output_value;
+
+		tape_codi.registerOutput(J);
+		#if _CODIPACK_MAJOR_==2
+		codi_global.output_indices.push_back(J.getIdentifier());
+		#elif _CODIPACK_MAJOR_==1
+		codi_global.output_indices.push_back(J.getGradientData());
+		#else
+		#error "_CODIPACK_MAJOR_ not supported"
+		#endif
+
+		/*Keep track of output for printing*/
+		Jlist[count] = output_value.getValue();
+      count++;
+	}
+	Jlist[count] = J.getValue();
+	_assert_(count == num_responses);
 
 	tape_codi.setPassive();
-	if(IssmComm::GetRank()==0) J.gradient() = 1.0;
-	tape_codi.evaluate();
+
+	if(VerboseAutodiff())_printf0_("   CoDiPack fos_reverse\n");
+	for(int i=0;i<num_responses;i++){
+		if(my_rank==0) tape_codi.setGradient(codi_global.output_indices[i],1.0);
+		tape_codi.evaluate();
+	}
 
 	/*Initialize Xb and Yb*/
@@ -449,11 +474,9 @@
 
 			/*Go through our dependent variables, and compute the response:*/
-			if(thistime==finaltime){
-				DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
-				for(Object* & object:dependent_objects->objects){
-					DependentObject* dep=(DependentObject*)object;
-					dep->Responsex(&output_value,femmodel);
-					dep->AddValue(output_value);
-				}
+			DataSet* dependent_objects=((DataSetParam*)femmodel->parameters->FindParamObject(AutodiffDependentObjectsEnum))->value;
+			for(Object* & object:dependent_objects->objects){
+				DependentObject* dep=(DependentObject*)object;
+				dep->Responsex(&output_value,femmodel);
+				dep->AddValue(output_value);
 			}
 
@@ -484,22 +507,17 @@
 		for(int i=0; i<Xsize;   i++) Xb[i] += X[i].gradient();
 	}
+
 	/*Clear tape*/
 	tape_codi.reset();
 
 	/*Broadcast gradient to other ranks (make sure to sum all gradients)*/
-	double* gradient=xNew<double>(Xsize);
-	ISSM_MPI_Allreduce(Xb,gradient,Xsize,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-
-	/*Now set final gradient*/
-	IssmDouble* aG=xNew<IssmDouble>(Xsize);
-	for(int i=0;i<Xsize;i++) aG[i] = reCast<IssmDouble>(gradient[i]);
-	xDelete<double>(gradient);
-	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
-	xDelete<IssmDouble>(aG);
-
+	ISSM_MPI_Allreduce(Xb,G,Xsize,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+
+	/*Cleanup and return misfit*/
 	xDelete<IssmDouble>(X);
 	xDelete<double>(Xb);
 	xDelete<double>(Yb);
 	xDelete<int>(Yin);
+   return J.getValue();
 }/*}}}*/
 #endif
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp	(revision 26888)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp	(revision 26889)
@@ -311,27 +311,28 @@
 				/*cfsurfacelogvel variables: */
 				int          num_cfsurfacelogvels;
-				char**       cfsurfacelogvel_name						= NULL;    
-				char**		 cfsurfacelogvel_definitionstring		= NULL;    
-				IssmDouble** cfsurfacelogvel_vxobs			= NULL;
-				IssmDouble** cfsurfacelogvel_vyobs			= NULL;
-				char**		 cfsurfacelogvel_vxobs_string	= NULL;
-				char**		 cfsurfacelogvel_vyobs_string	= NULL;
-				int*         cfsurfacelogvel_observation_M			= NULL;
-				int*         cfsurfacelogvel_observation_N			= NULL;
-				IssmDouble** cfsurfacelogvel_weights					= NULL;
-				int*         cfsurfacelogvel_weights_M				= NULL;
-				int*         cfsurfacelogvel_weights_N				= NULL;
-				char**       cfsurfacelogvel_weightstring		= NULL;
-				IssmDouble*  cfsurfacelogvel_datatime				= NULL;
-
-				/*Fetch name, modeltring, observation, observationtring, etc ... (see src/m/classes/cfsurfacelogvel.m): */
-				iomodel->FetchMultipleData(&cfsurfacelogvel_name,&num_cfsurfacelogvels,                                                        "md.cfsurfacelogvel.name");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_definitionstring,&num_cfsurfacelogvels,                                            "md.cfsurfacelogvel.definitionstring");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vxobs");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs_string,&num_cfsurfacelogvels,                                          "md.cfsurfacelogvel.vxobs_string");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vyobs");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels,                                          "md.cfsurfacelogvel.vyobs_string");			iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels,             "md.cfsurfacelogvel.weights");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels,                                              "md.cfsurfacelogvel.weights_string");
-				iomodel->FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels,																	 "md.cfsurfacelogvel.datatime");
+				char       **cfsurfacelogvel_name             = NULL;
+				char       **cfsurfacelogvel_definitionstring = NULL;
+				IssmDouble **cfsurfacelogvel_vxobs            = NULL;
+				IssmDouble **cfsurfacelogvel_vyobs            = NULL;
+				char       **cfsurfacelogvel_vxobs_string     = NULL;
+				char       **cfsurfacelogvel_vyobs_string     = NULL;
+				int         *cfsurfacelogvel_observation_M    = NULL;
+				int         *cfsurfacelogvel_observation_N    = NULL;
+				IssmDouble **cfsurfacelogvel_weights          = NULL;
+				int         *cfsurfacelogvel_weights_M        = NULL;
+				int         *cfsurfacelogvel_weights_N        = NULL;
+				char       **cfsurfacelogvel_weightstring     = NULL;
+				IssmDouble  *cfsurfacelogvel_datatime         = NULL;
+
+            /*Fetch name, modeltring, observation, observationtring, etc ... (see src/m/classes/cfsurfacelogvel.m): */
+            iomodel->FetchMultipleData(&cfsurfacelogvel_name,&num_cfsurfacelogvels,"md.cfsurfacelogvel.name");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_definitionstring,&num_cfsurfacelogvels,"md.cfsurfacelogvel.definitionstring");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vxobs");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs_string,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vxobs_string");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs,NULL,NULL,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vyobs");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels,"md.cfsurfacelogvel.vyobs_string");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels,"md.cfsurfacelogvel.weights");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels,"md.cfsurfacelogvel.weights_string");
+            iomodel->FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels,"md.cfsurfacelogvel.datatime");
 
 				for(j=0;j<num_cfsurfacelogvels;j++){