Index: /issm/trunk/src/m/solutions/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/control_core.m	(revision 5462)
+++ /issm/trunk/src/m/solutions/control_core.m	(revision 5463)
@@ -19,4 +19,5 @@
 	cm_gradient=femmodel.parameters.CmGradient;
 	dim=femmodel.parameters.Dim;
+	isstokes=femmodel.parameters.IsStokes;
 
 	%Initialise options with tolerance and maxiter
@@ -31,8 +32,10 @@
 	adjointcore=AdjointCorePointerFromSolutionEnum(solution_type);
 
-	%Preprocess models
-%	if(solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum)
-%		femmodel=stokescontrolinit(femmodel);
-%	end
+	%Launch once a complete solution to set up all inputs
+	disp('      preparing initial solution');
+	isstokes
+	if isstokes,
+		eval(['femmodel=' solutioncore '(femmodel);']);
+	end
 
 	for n=1:nsteps,
Index: sm/trunk/src/m/solutions/stokescontrolinit.m
===================================================================
--- /issm/trunk/src/m/solutions/stokescontrolinit.m	(revision 5462)
+++ 	(revision )
@@ -1,31 +1,0 @@
-function femmodel=stokescontrolinit(femmodel),
-%STOKESCONTROLINIT - initialize femmodel if control method on Stokes
-
-	conserve_loads=true;
-	%recover parameters common to all solutions
-	verbose=femmodel.parameters.Verbose;
-	isstokes=femmodel.parameters.IsStokes;
-	stokesreconditioning=femmodel.parameters.StokesReconditioning;
-
-	%if no Stokes analysis carried out, just return
-	if (isstokes==0),
-		return;
-	end
-
-	% For Stokes inverse control method, we are going to carry out the inversion only on the Stokes part. So we need to solve the Hutter or MacAyeal/Pattyn femmodel here, and constrain the Stokes femmodel using the Hutter or MacAyeal/Pattyn at the boundary. We don't want to have to do that at every inversion step, as it needs be done only once:
-
-	%Compute slopes:
-	femmodel=bedslope_core(femmodel);
-
-	%Run a complete diagnostic to update the Stokes spcs:
-	femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-	femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);
-	femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-	femmodel=solver_linear(femmodel);
-
-	displaystring(verbose,'\n%s',['      update boundary conditions for stokes using velocities previously computed...']);
-	femmodel=ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum);
-
-	displaystring(verbose,'\n%s',['      computing stokes velocity and pressure...']);
-	femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum);
-	femmode=solver_diagnostic_nonlinear(femmodel,conserve_loads);