Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 25475)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 25476)
@@ -277,4 +277,25 @@
 }
 /*}}}*/
+void FemModel::CheckPointAD(int step){/*{{{*/
+
+	/*Get rank*/
+	int my_rank = IssmComm::GetRank();
+
+	/*Get string sizes*/
+	int rank_length = (my_rank == 0 ? 1 : int(log10(static_cast<double>(my_rank))+1));
+	int step_length = (step    == 0 ? 1 : int(log10(static_cast<double>(step))   +1));
+
+	/*Create restart file*/
+	char* restartfilename  = xNew<char>(strlen("AD_step_")+step_length+strlen("_rank_")+rank_length+strlen(".ckpt")+1);
+	sprintf(restartfilename,"%s%i%s%i%s","AD_step_",step,"_rank_",my_rank,".ckpt");
+	this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename));
+
+	/*Write files*/
+	this->CheckPoint();
+
+	/*Clean up and return*/
+	xDelete<char>(restartfilename);
+
+}/*}}}*/
 void FemModel::CleanUp(void){/*{{{*/
 
@@ -600,4 +621,24 @@
 	xDelete<char>(restartfilename);
 	xDelete<char>(femmodel_buffer);
+}/*}}}*/
+void FemModel::RestartAD(int step){ /*{{{*/
+
+	/*Get rank*/
+	int my_rank = IssmComm::GetRank();
+
+	/*Get string sizes*/
+	int rank_length = (my_rank == 0 ? 1 : int(log10(static_cast<double>(my_rank))+1));
+	int step_length = (step    == 0 ? 1 : int(log10(static_cast<double>(step))   +1));
+
+	/*Create restart file*/
+	char* restartfilename  = xNew<char>(strlen("AD_step_")+step_length+strlen("_rank_")+rank_length+strlen(".ckpt")+1);
+	sprintf(restartfilename,"%s%i%s%i%s","AD_step_",step,"_rank_",my_rank,".ckpt");
+	this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename));
+
+	/*Read files*/
+	this->Restart();
+
+	/*Clean up and return*/
+	xDelete<char>(restartfilename);
 }/*}}}*/
 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 25475)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 25476)
@@ -74,4 +74,5 @@
 		int  AnalysisIndex(int);
 		void CheckPoint(void);
+		void CheckPointAD(int step);
 		void CleanUp(void);
 		FemModel* copy();
@@ -82,4 +83,5 @@
 		void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction);
 		void Restart(void);
+		void RestartAD(int step);
 		void SetCurrentConfiguration(int configuration_type);
 		void SetCurrentConfiguration(int configuration_type,int analysis_type);
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25475)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25476)
@@ -272,18 +272,9 @@
 
 	/*parameters: */
-	IssmDouble finaltime,dt,yts;
-	bool       isoceancoupling,iscontrol,isautodiff,isslr;
-	int        timestepping;
-	int        output_frequency;
-	int        recording_frequency;
-	int        amr_frequency,amr_restart;
-
-	/*intermediary: */
-	int        step;
-	IssmDouble time;
-
-	/*first, figure out if there was a check point, if so, do a reset of the FemModel* femmodel structure. */
-	femmodel->parameters->FindParam(&recording_frequency,SettingsRecordingFrequencyEnum);
-	if(recording_frequency) femmodel->Restart();
+	IssmDouble output_value;
+	IssmDouble finaltime,dt,yts,time;
+	bool       isoceancoupling,isslr;
+	int        step,timestepping;
+	DataSet*   dependent_objects=NULL;
 
 	/*then recover parameters common to all solutions*/
@@ -292,16 +283,7 @@
 	femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
 	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
-	femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
 	femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
 	femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
-	femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
-	femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
-	femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
-
-	DataSet* dependent_objects=NULL;
-	if(iscontrol){
-		femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
-	}
-
+	femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
 	if(isslr) sealevelrise_core_geometry(femmodel);
 
@@ -309,16 +291,6 @@
 
 		/*Time Increment*/
-		switch(timestepping){
-			case AdaptiveTimesteppingEnum:
-				femmodel->TimeAdaptx(&dt);
-				if(time+dt>finaltime) dt=finaltime-time;
-				femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
-				break;
-			case FixedTimesteppingEnum:
-				femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-				break;
-			default:
-				_error_("Time stepping \""<<EnumToStringx(timestepping)<<"\" not supported yet");
-		}
+		if(timestepping!=FixedTimesteppingEnum) _error_("not supported yet, but easy to handle...");
+		femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 		step+=1;
 		time+=dt;
@@ -327,16 +299,14 @@
 		femmodel->parameters->SetParam(false,SaveResultsEnum);
 
+		/*Store Model State at the beginning of the step*/
+		if(VerboseSolution()) _printf0_("   checkpointing model \n");
+		femmodel->CheckPointAD(step);
+
 		/*Run transient step!*/
 		transient_step(femmodel);
-
-		if(recording_frequency && (step%recording_frequency==0)){
-			if(VerboseSolution()) _printf0_("   checkpointing model \n");
-			femmodel->CheckPoint();
-		}
 
 		/*Go through our dependent variables, and compute the response:*/
 		for(int i=0;i<dependent_objects->Size();i++){
 			DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
-			IssmDouble  output_value;
 			dep->Responsex(&output_value,femmodel);
 			dep->AddValue(output_value);
@@ -344,6 +314,91 @@
 	}
 
-	if(!iscontrol) femmodel->RequestedDependentsx();
-	if(iscontrol) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum);
+	/*__________________________________________________________________________________*/
+
+	/*Get X (control) and Y (model state)*/
+	IssmDouble *X = NULL; int Xsize;
+	GetVectorFromControlInputsx(&X,&Xsize,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+	IssmDouble *Y = NULL; int Ysize;
+	_error_("don't know how to get Y and Ysize here...");
+
+	/*Initialize Xb, Yb and Yin*/
+	double *Xb = xNewZeroInit<double>(Xsize);
+	double *Yb  = xNewZeroInit<double>(Ysize);
+	double *Yin = xNewZeroInit<double>(Ysize);
+
+	/*Start tracing*/
+	auto& tape_codi = IssmDouble::getGlobalTape();
+	tape_codi.setActive();
+
+	/*Reverse dependent (f)*/
+	for(int i=0; i < Ysize; i++) tape_codi.registerInput(Y[i]);
+	for(int i=0; i < Xsize; i++) tape_codi.registerInput(X[i]);
+	SetControlInputsFromVectorx(femmodel,X);
+
+	IssmDouble J = 0.;
+	for(int i=0;i<dependent_objects->Size();i++){
+		DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
+		J += dep->GetValue();
+	}
+
+	if(IssmComm::GetRank()==0) {
+		tape_codi.registerOutput(J);
+	}
+	tape_codi.setPassive();
+
+	J.gradient() = 1.0;
+	tape_codi.evaluate();
+
+	for(int i=0;i<Xsize;i++) Xb[i] += X[i].gradient();
+	for(int i=0;i<Ysize;i++) Yb[i]  = Y[i].gradient();
+
+	/*reverse loop for transient step (G)*/
+	for(int reverse_step = step;reverse_step>=0; reverse_step--){
+
+		/*Restore model from this step*/
+		tape_codi.reset();
+		femmodel->RestartAD(reverse_step);
+		tape_codi.setActive();
+
+		/*We need to store the CoDiPack identifier here, since y is overwritten.*/
+		for(int i=0; i<Ysize; i++) {
+			tape_codi.registerInput(Y[i]);
+			Yin[i] = Y[i].getGradientData();
+		}
+
+		/*Tell codipack that X is the independent*/
+		for(int i=0; i<Xsize; i++) tape_codi.registerInput(X[i]);
+		SetControlInputsFromVectorx(femmodel,X);
+
+		/*Get New state*/
+		transient_step(femmodel);
+		_error_("don't know how to get new Y again here...");
+		for(int i=0; i<Ysize; i++) tape_codi.registerOutput(Y[i]);
+
+		/*stop tracing*/
+		tape_codi.setPassive();
+
+		/*Reverse transient step (G)*/
+		/* Using y_b here to seed the next reverse iteration there y_b is always overwritten*/
+		for(int i=0; i<Ysize; i++) Y[i].gradient() = Yb[i]; 
+		tape_codi.evaluate();
+		/* here we access the gradient data via the stored identifiers.*/
+		for(int i=0; i<Ysize; i++) Xb[i]  = tape_codi.gradient(Yin[i]);
+		for(int i=0; i<Xsize; i++) Yb[i] += X[i].gradient();
+	}
+
+	/*Now set final gradient*/
+	IssmDouble* aG=xNew<IssmDouble>(Xsize);
+	for(int i=0;i<Xsize;i++){
+		aG[i] = reCast<IssmDouble>(Xb[i]);
+	}
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
+	xDelete<IssmDouble>(aG);
+
+	xDelete<IssmDouble>(X);
+	xDelete<double>(Xb);
+	xDelete<IssmDouble>(Y);
+	xDelete<double>(Yb);
+	xDelete<double>(Yin);
 }/*}}}*/
 #endif