Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 23267)
+++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp	(revision 23268)
@@ -17,5 +17,5 @@
 
 #if defined (_HAVE_M1QN3_) && defined(_HAVE_AD_)
-/*m1qn3 prototypes*/
+/*m1qn3 prototypes {{{*/
 extern "C" void *ctonbe_; // DIS mode : Conversion
 extern "C" void *ctcabe_; // DIS mode : Conversion
@@ -36,4 +36,5 @@
 	int*         i;
 } m1qn3_struct;
+/*}}}*/
 
 /*m1qm3 functions*/
@@ -91,4 +92,6 @@
 		tape_codi.reset();
 		//alloc_profiler.Tag(FinishInit, true);
+		#else
+		tape_codi.reset();
 		#endif
 
@@ -181,9 +184,6 @@
 
 	FemModel* femmodel = input_struct->femmodel;
-	int num_responses,num_controls,numberofvertices,solution_type;
+	int num_responses,num_controls,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*/
@@ -197,5 +197,7 @@
 
 	/*Recover some parameters*/
-	IssmDouble* scaling_factors = NULL;
+	IssmDouble *scaling_factors = NULL;
+	int        *N               = NULL;
+	int        *control_enum    = NULL;
 	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
@@ -203,5 +205,5 @@
 	femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
 	femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
-	numberofvertices=femmodel->vertices->NumberOfVertices();
+	int numberofvertices=femmodel->vertices->NumberOfVertices();
 
 	/*Constrain input vector and update controls*/
@@ -211,5 +213,5 @@
 	GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
 
-	N_add = 0;
+	int N_add = 0;
 	for (int c=0;c<num_controls;c++){
 		for(int i=0;i<numberofvertices*N[c];i++){
@@ -304,17 +306,9 @@
 	//if(my_rank==0) _printf_("\nTIME: "<<now<<"\n");
 
-	/*diverse: */
-	int  dummy;
-	int  num_independents=0;
-
 	/*intermediary: */
+	int          num_independents=intn;
 	IssmPDouble *aWeightVector=NULL;
 	IssmPDouble *weightVectorTimesJac=NULL;
-
-	/*output: */
-	IssmPDouble *totalgradient=NULL;
-
-	/*retrieve parameters: */
-	num_independents = intn;
+	IssmPDouble *totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
 
 	/*if no dependents, no point in running a driver: */
@@ -341,6 +335,4 @@
 	/* 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: */
-	totalgradient=xNewZeroInit<IssmPDouble>(num_independents_old);
-
 	for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
 
@@ -379,31 +371,27 @@
 	/*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);
+
+	/* 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)  tape_codi.setGradient(codi_global.output_indices[aDepIndex],1.0);
+		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) {
+			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];
+		}
+		/*free resources :*/
+		xDelete(weightVectorTimesJac);
+		xDelete(aWeightVector);
+	}
 	/*}}}*/
 	#else
@@ -426,10 +414,4 @@
 	Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
 
-/*
-	IssmDouble* test = xNew<IssmDouble>(intn);
-	femmodel->parameters->FindParam(&test,NULL,InversionXBestEnum);
-	for(int i=0;i<10;i++)_printf_("X "<<test[i]<<"\n");
-	xDelete<IssmDouble>(intn);
-*/
 	if(*indic==0){
 		/*dry run, no gradient required*/
Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 23267)
+++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 23268)
@@ -12,5 +12,5 @@
 
 #if defined (_HAVE_M1QN3_) && !defined(_HAVE_AD_)
-/*m1qn3 prototypes*/
+/*m1qn3 prototypes {{{*/
 extern "C" void *ctonbe_; // DIS mode : Conversion
 extern "C" void *ctcabe_; // DIS mode : Conversion
@@ -34,6 +34,7 @@
 	int*        i;
 } m1qn3_struct;
-
-void controlm1qn3_core(FemModel* femmodel){
+/*}}}*/
+
+void controlm1qn3_core(FemModel* femmodel){/*{{{*/
 
 	/*Intermediaries*/
@@ -173,8 +174,6 @@
 	xDelete<double>(mystruct.Jlist);
 	xDelete<int>(mystruct.i);
-}
-
-/*Cost function definition*/
-void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
+}/*}}}*/
+void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
 
 	/*Recover Arguments*/
@@ -278,9 +277,7 @@
 	xDelete<IssmDouble>(XL);
 	xDelete<IssmDouble>(scaling_factors);
-}
+}/*}}}*/
 
 #else
-void controlm1qn3_core(FemModel* femmodel){
-	_error_("M1QN3 not installed");
-}
+void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
 #endif //_HAVE_M1QN3_