Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25463)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 25464)
@@ -17,9 +17,12 @@
 #include "../solutionsequences/solutionsequences.h"
 
-void transient_core(FemModel* femmodel){
+/*Prototypes*/
+void transient_step(FemModel* femmodel);
+
+void transient_core(FemModel* femmodel){/*{{{*/
 
 	/*parameters: */
 	IssmDouble finaltime,dt,yts,starttime;
-	bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isesa,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff;
+	bool       isoceancoupling,iscontrol,isautodiff,isgroundingline,isslr;
 	bool       save_results,dakota_analysis;
 	int        timestepping;
@@ -27,7 +30,6 @@
 	int        sb_coupling_frequency;
 	int        recording_frequency;
-	int        domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart;
+	int        groundingline_migration,smb_model,amr_frequency,amr_restart;
 	int        numoutputs;
-	Analysis  *analysis          = NULL;
 	char     **requested_outputs = NULL;
 
@@ -41,5 +43,4 @@
 
 	/*then recover parameters common to all solutions*/
-	femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
 	femmodel->parameters->FindParam(&step,StepEnum);
 	femmodel->parameters->FindParam(&time,TimeEnum);
@@ -51,19 +52,8 @@
 	femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
 	femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
-	femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
-	femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
-	femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
-	femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
-	femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
-	femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
 	femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
-	femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
 	femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
-	femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
 	femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
-	femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
-	femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
 	femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
-	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
 	femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
 	femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
@@ -128,105 +118,6 @@
 		femmodel->parameters->SetParam(save_results,SaveResultsEnum);
 
-	#if defined(_HAVE_OCEAN_ )
-	if(isoceancoupling) OceanExchangeDatax(femmodel,false);
-	#endif
-
-		if(isthermal && domaintype==Domain3DEnum){
-			if(issmb){
-				bool isenthalpy;
-				femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
-				femmodel->parameters->FindParam(&smb_model,SmbEnum);
-				if(isenthalpy){
-					if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
-						femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
-						ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
-					}
-				}
-				else{
-					if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
-						femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
-						ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
-					}
-				}
-			}
-			thermal_core(femmodel);
-		}
-		/* Using Hydrology dc  coupled we need to compute smb in the hydrology inner time loop*/
-		if(issmb) {
-			if(VerboseSolution()) _printf0_("   computing smb\n");
-			smb_core(femmodel);
-		}
-
-		if(ishydrology){
-			if(VerboseSolution()) _printf0_("   computing hydrology\n");
-			int hydrology_model;
-			hydrology_core(femmodel);
-			femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
-			if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
-		}
-
-		if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
-			if(VerboseSolution()) _printf0_("   computing stress balance\n");
-			stressbalance_core(femmodel);
-		}
-
-		if(isdamageevolution) {
-			if(VerboseSolution()) _printf0_("   computing damage\n");
-			damage_core(femmodel);
-		}
-
-		if(ismovingfront)	{
-			if(VerboseSolution()) _printf0_("   computing moving front\n");
-			movingfront_core(femmodel);
-		}
-
-		/* from here on, prepare geometry for next time step*/
-		//if(issmb) smb_core(femmodel);
-
-		if(ismasstransport){
-			if(VerboseSolution()) _printf0_("   computing mass transport\n");
-			bmb_core(femmodel);
-			masstransport_core(femmodel);
-			femmodel->UpdateVertexPositionsx();
-		}
-
-		if(isgroundingline){
-
-			/*Start profiler*/
-			femmodel->profiler->Start(GROUNDINGLINECORE);
-
-			if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
-			GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-
-			femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
-			extrudefrombase_core(femmodel);
-			femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
-			extrudefrombase_core(femmodel);
-			femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
-			extrudefrombase_core(femmodel);
-
-			/*Stop profiler*/
-			femmodel->profiler->Stop(GROUNDINGLINECORE);
-
-			if(save_results){
-				int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
-				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
-			}
-		}
-
-		if(isgia){
-			if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
-			#ifdef _HAVE_GIA_
-			gia_core(femmodel);
-			#else
-			_error_("ISSM was not compiled with gia capabilities. Exiting");
-			#endif
-		}
-
-		/*esa: */
-		if(isesa) esa_core(femmodel);
-
-		/*Sea level rise: */
-		if(isslr) sealevelchange_core(femmodel);
+		/*Run transient step!*/
+		transient_step(femmodel);
 
 		/*unload results*/
@@ -263,5 +154,5 @@
 		}
 
-		if (iscontrol && isautodiff) {
+		if(iscontrol && isautodiff){
 			/*Go through our dependent variables, and compute the response:*/
 			for(int i=0;i<dependent_objects->Size();i++){
@@ -278,3 +169,132 @@
 	/*Free ressources:*/
 	if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
-}
+
+}/*}}}*/
+void transient_step(FemModel* femmodel){/*{{{*/
+
+	/*parameters: */
+	bool isstressbalance,ismasstransport,issmb,isthermal,isgroundingline,isgia,isesa;
+	bool isslr,ismovingfront,isdamageevolution,ishydrology,isoceancoupling;
+	bool save_results;
+	int  step,sb_coupling_frequency;
+	int  domaintype,smb_model;
+
+	/*then recover parameters common to all solutions*/
+	femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
+	femmodel->parameters->FindParam(&step,StepEnum);
+	femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
+	femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
+	femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
+	femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
+	femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
+	femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
+	femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
+	femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
+	femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
+	femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
+	femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
+	femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
+	femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
+
+	#if defined(_HAVE_OCEAN_ )
+	if(isoceancoupling) OceanExchangeDatax(femmodel,false);
+	#endif
+
+	if(isthermal && domaintype==Domain3DEnum){
+		if(issmb){
+			bool isenthalpy;
+			femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+			femmodel->parameters->FindParam(&smb_model,SmbEnum);
+			if(isenthalpy){
+				if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
+					femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
+					ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
+				}
+			}
+			else{
+				if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
+					femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+					ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
+				}
+			}
+		}
+		thermal_core(femmodel);
+	}
+	/* Using Hydrology dc  coupled we need to compute smb in the hydrology inner time loop*/
+	if(issmb) {
+		if(VerboseSolution()) _printf0_("   computing smb\n");
+		smb_core(femmodel);
+	}
+
+	if(ishydrology){
+		if(VerboseSolution()) _printf0_("   computing hydrology\n");
+		int hydrology_model;
+		hydrology_core(femmodel);
+		femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
+		if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
+	}
+
+	if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
+		if(VerboseSolution()) _printf0_("   computing stress balance\n");
+		stressbalance_core(femmodel);
+	}
+
+	if(isdamageevolution) {
+		if(VerboseSolution()) _printf0_("   computing damage\n");
+		damage_core(femmodel);
+	}
+
+	if(ismovingfront)	{
+		if(VerboseSolution()) _printf0_("   computing moving front\n");
+		movingfront_core(femmodel);
+	}
+
+	/* from here on, prepare geometry for next time step*/
+
+	if(ismasstransport){
+		if(VerboseSolution()) _printf0_("   computing mass transport\n");
+		bmb_core(femmodel);
+		masstransport_core(femmodel);
+		femmodel->UpdateVertexPositionsx();
+	}
+
+	if(isgroundingline){
+
+		/*Start profiler*/
+		femmodel->profiler->Start(GROUNDINGLINECORE);
+
+		if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
+		GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+		femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
+		extrudefrombase_core(femmodel);
+		femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
+		extrudefrombase_core(femmodel);
+		femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
+		extrudefrombase_core(femmodel);
+
+		/*Stop profiler*/
+		femmodel->profiler->Stop(GROUNDINGLINECORE);
+
+		if(save_results){
+			int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
+		}
+	}
+
+	if(isgia){
+		if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
+		#ifdef _HAVE_GIA_
+		gia_core(femmodel);
+		#else
+		_error_("ISSM was not compiled with gia capabilities. Exiting");
+		#endif
+	}
+
+	/*esa: */
+	if(isesa) esa_core(femmodel);
+
+	/*Sea level rise: */
+	if(isslr) sealevelchange_core(femmodel);
+
+}/*}}}*/