Index: sm/trunk/src/c/solutions/ControlInitialization.cpp
===================================================================
--- /issm/trunk/src/c/solutions/ControlInitialization.cpp	(revision 4051)
+++ 	(revision )
@@ -1,175 +1,0 @@
-/*!\file: ControlInitialization.cpp
- * \brief: ...
- */ 
-
-#include "../toolkits/toolkits.h"
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "./solutions.h"
-#include "../modules/modules.h"
-
-void ControlInitialization(Model* model){
-
-	extern int my_rank;
-
-	/*fem models: */
-	FemModel* fem_dh=NULL;
-	FemModel* fem_dv=NULL;
-	FemModel* fem_dhu=NULL;
-	FemModel* fem_ds=NULL;
-	FemModel* fem_sl=NULL;
-
-	/*solutions: */
-	Vec ug=NULL;
-	Vec ug_horiz=NULL;
-	Vec ug_vert=NULL;
-	Vec ug_stokes=NULL;
-	Vec pg=NULL;
-	Vec slopex=NULL;
-	Vec slopey=NULL;
-	double* vx=NULL;
-	double* vy=NULL;
-	double* vz=NULL;
-	double* pressure=NULL;
-
-	/*flags: */
-	int verbose=0;
-	int dim=-1;
-	int ishutter=0;
-	int ismacayealpattyn=0;
-	int isstokes=0;
-	int numberofdofspernode_sl;
-	int numberofdofspernode_dh;
-	int numberofdofspernode_ds;
-	int numberofnodes;
-
-	double stokesreconditioning;
-
-	/*dof recovery: */
-	int dof01[2]={0,1};
-	int dof2[1]={2};
-	int dof012[3]={0,1,2};
-	int dof3[1]={3};
-	double* dofset=NULL;
-
-	/*first recover parameters common to all solutions:*/
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&dim,DimEnum);
-	model->FindParam(&ishutter,IsHutterEnum);
-	model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum);
-	model->FindParam(&isstokes,IsStokesEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&stokesreconditioning,StokesReconditioningEnum);
-
-	/*recover fem models: */
-	fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum);
-	fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum);
-
-	//specific parameters for specific models
-	fem_dh->parameters->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum);
-	fem_sl->parameters->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum);
-	fem_ds->parameters->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum);
-
-	/*if no Stokes, assign output and return*/
-	if (!isstokes){
-		model->SetActiveFormulation(fem_dh);
-		return;
-	}
-
-	/*1: compute slopes once for all*/
-
-	//compute slopes
-	if(verbose)_printf_("%s\n","computing bed slope (x and y derivatives)...");
-	diagnostic_core_linear(&slopex,fem_sl,SlopecomputeAnalysisEnum,BedSlopeXAnalysisEnum);
-	diagnostic_core_linear(&slopey,fem_sl,SlopecomputeAnalysisEnum,BedSlopeYAnalysisEnum);
-	FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->vertices,fem_sl->loads,fem_sl->materials,fem_sl->parameters,"slopex",0);
-	FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->vertices,fem_sl->loads,fem_sl->materials,fem_sl->parameters,"slopey",0);
-
-	//Add in inputs
-	model->UpdateInputsFromVector(slopex,BedSlopexEnum,VertexEnum);
-	model->UpdateInputsFromVector(slopey,BedSlopeyEnum,VertexEnum);
-	VecFree(&slopex); VecFree(&slopey);
-
-	/*2: run a complete diagnostic to update spcs*/
-
-	//horizontal velocity
-	if(verbose)_printf_("%s\n"," computing horizontal velocities...");
-	diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_dh,DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	if(verbose)_printf_("%s\n"," extruding horizontal velocities...");
-	VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh-> materials,fem_dh->parameters,"velocity",1);
-
-	//Add to inputs:
-	SplitSolutionVectorx(ug_horiz,numberofnodes,numberofdofspernode_dh,&vx,&vy);
-	model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-	model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-	
-
-	//vertical velocity
-	if(verbose)_printf_("%s\n"," computing vertical velocities...");
-	diagnostic_core_linear(&ug_vert,fem_dv,DiagnosticAnalysisEnum,VertAnalysisEnum);
-
-	//Create 3d u_g
-	if(verbose)_printf_("%s\n"," combining horizontal and vertical velocities...");
-	VecFree(&ug); ug=NewVec(numberofnodes*3);
-	xfree((void**)&dofset);dofset=dofsetgen(2,&dof01[0],3,numberofnodes*3); VecMerge(ug,ug_horiz,dofset,numberofnodes*2);
-	xfree((void**)&dofset);dofset=dofsetgen(1,&dof2[0],3,numberofnodes*3); VecMerge(ug,ug_vert,dofset,numberofnodes*1);
-	VecFree(&ug_vert); VecFree(&ug_horiz);
-
-	//Create 4d u_g
-	if(verbose)_printf_("%s\n"," computing pressure according to Pattyn...");
-	ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials, fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	VecScale(pg,1.0/stokesreconditioning);
-	ug_stokes=NewVec(fem_ds->nodesets->GetGSize());
-	xfree((void**)&dofset);dofset=dofsetgen(3,dof012,4,numberofnodes*4); VecMerge(ug_stokes,ug,dofset,numberofnodes*3);
-	xfree((void**)&dofset);dofset=dofsetgen(1,dof3,4,numberofnodes*4); VecMerge(ug_stokes,pg,dofset,numberofnodes);
-
-	//Add in inputs
-	SplitSolutionVectorx(ug_stokes,numberofnodes,numberofdofspernode_ds,&vx,&vy,&vz,&pressure);
-	model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-	model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-	model->UpdateInputsFromVector(vz,VzEnum,VertexEnum);
-	model->UpdateInputsFromVector(pressure,PressureEnum,VertexEnum);
-	VecFree(&ug_stokes);
-
-	//update spcs
-	if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed...");
-	xfree((void**)&dofset);dofset=dofsetgen(3,dof012,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,3*numberofnodes);
-	
-	VecFree(&fem_ds->ys);
-	Reducevectorgtosx(&fem_ds->ys,fem_ds->yg->vector,fem_ds->nodesets);
-
-	//Compute Stokes velocities to speed up later runs
-	if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
-	VecFree(&ug);
-	diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum);
-
-	//Add in inputs
-	xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); xfree((void**)&pressure);
-	SplitSolutionVectorx(ug,numberofnodes,numberofdofspernode_ds,&vx,&vy,&vz,&pressure);
-	model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-	model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-	model->UpdateInputsFromVector(vz,VzEnum,VertexEnum);
-	model->UpdateInputsFromVector(pressure,PressureEnum,VertexEnum);
-
-	/*Assign output*/
-	model->SetActiveFormulation(fem_ds);
-
-	/*Free ressources:*/
-	xfree((void**)&dofset);
-	xfree((void**)&vx);
-	xfree((void**)&vy);
-	xfree((void**)&vz);
-	xfree((void**)&pressure);
-	VecFree(&ug);
-	VecFree(&ug_horiz);
-	VecFree(&ug_vert);
-	VecFree(&ug_stokes);
-	VecFree(&pg);
-	VecFree(&slopex);
-	VecFree(&slopey);
-
-}
Index: sm/trunk/src/c/solutions/ControlRestart.cpp
===================================================================
--- /issm/trunk/src/c/solutions/ControlRestart.cpp	(revision 4051)
+++ 	(revision )
@@ -1,26 +1,0 @@
-/*!\file: ControlRestart.cpp
- * \brief: save as much as possible, to be able to restart the control_core solution
- */ 
-
-#include "./solutions.h"
-#include "../modules/modules.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-
-void ControlRestart(FemModel* femmodel){
-
-	char*    outputfilename=NULL;
-	
-	/*retrieve output file name: */
-	model->FindParam(&outputfilename,OutputFileNameEnum);
-
-	/*we essentially want J and the parameter: */
-	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
-	femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
-	femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
-
-	/*write to disk: */
-	OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
-	
-	/*Free ressources:*/
-	xfree((void**)&outputfilename);
-}
Index: /issm/trunk/src/c/solutions/control_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/control_core.cpp	(revision 4051)
+++ /issm/trunk/src/c/solutions/control_core.cpp	(revision 4052)
@@ -41,8 +41,8 @@
 	double* J=NULL;
 
-	/*Process models*/
-	ControlInitialization(femmodel);
+	/*some preliminary work to be done if running full-Stokes analysis: */
+	stokescontrolinit(femmodel);
 
-	/*Recover parameters used throughout the solution:*/
+	/*Recover parameters used throughout the solution:{{{1*/
 	femmodel->parameters->FindParam(&nsteps,NStepsEnum);
 	femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
@@ -57,4 +57,5 @@
 	femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
 	femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
+	/*}}}*/
 
 	/*Initialize misfit: */
@@ -69,15 +70,15 @@
 
 		_printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-		femmodel->UpdateInputsFromConstant(fit[n],FitEnum);
+		UpdateInputsFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,fit[n],FitEnum);
 		
 		/*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */
-		if (control_steady) steadystate_core(model);
+		if (control_steady) steadystate_core(femmodel);
 	
-		if(verbose)_printf_("%s\n","      computing gradJ...");
+		_printf_("%s\n","      computing gradJ...");
 		gradient_core(femmodel);
 
 		/*Return gradient if asked: */
 		if (cm_gradient){
-			InputToResultx(femodel->elements,femodel->nodes,femodel->vertices,femodel->loads,femodel->materials,femodel->parameters,GradientEnum);
+			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum);
 			goto cleanup_and_return;
 		}
@@ -91,5 +92,5 @@
 
 		_printf_("%s","      constraining the new distribution...");    
-		ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
+		ControlConstrainInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
 		
 		_printf_("%s","      save new parameter...");
@@ -98,13 +99,15 @@
 		_printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
 		
-		if(controlconvergence)goto cleanup_and_return;
+		if(controlconvergence(J,fit,eps_cm,n))goto convergence_point;
 
 		/*Temporary saving every 5 control steps: */
 		if (((n+1)%5)==0){
 			_printf_("%s","      saving temporary results...");
-			ControlRestart(femmodel);
+			controlrestart(femmodel);
 		}
 	}
 
+
+	convergence_point:
 	_printf_("%s","      preparing final velocity solution");
 	if (control_steady) steadystate_core(femmodel);
@@ -113,9 +116,8 @@
 	/*some results not computed by steadystate_core or diagnostic_core: */
 	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
-	femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
-	femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
+	femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
+	femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumAsString(control_type),1,0));
 
 	cleanup_and_return:
-	
 	/*Free ressources: */
 	xfree((void**)&control_type);
Index: /issm/trunk/src/c/solutions/controlrestart.cpp
===================================================================
--- /issm/trunk/src/c/solutions/controlrestart.cpp	(revision 4052)
+++ /issm/trunk/src/c/solutions/controlrestart.cpp	(revision 4052)
@@ -0,0 +1,26 @@
+/*!\file: controlrestart.cpp
+ * \brief: save as much as possible, to be able to restart the control_core solution
+ */ 
+
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void controlrestart(FemModel* femmodel){
+
+	char*    outputfilename=NULL;
+	
+	/*retrieve output file name: */
+	model->FindParam(&outputfilename,OutputFileNameEnum);
+
+	/*we essentially want J and the parameter: */
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
+	femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
+	femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
+
+	/*write to disk: */
+	OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
+	
+	/*Free ressources:*/
+	xfree((void**)&outputfilename);
+}
Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4051)
+++ /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4052)
@@ -23,5 +23,5 @@
 	double stokesreconditioning;
 
-	/* recover parameters: */
+	/* recover parameters: {{{1*/
 	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 	femmodel->parameters->FindParam(&dim,DimEnum);
@@ -31,4 +31,5 @@
 	femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
 	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*}}}*/
 
 	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4051)
+++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4052)
@@ -54,5 +54,5 @@
 
 	/*Constrain:*/
-	ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
+	ControlConstrainInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
 
 	/*Run diagnostic with updated inputs: */
Index: /issm/trunk/src/c/solutions/solutions.h
===================================================================
--- /issm/trunk/src/c/solutions/solutions.h	(revision 4051)
+++ /issm/trunk/src/c/solutions/solutions.h	(revision 4052)
@@ -13,21 +13,21 @@
 
 /*cores: */
-Results* prognostic_core(FemModel* model);
-Results* prognostic2_core(FemModel* model);
-Results* balancedthickness_core(FemModel* model);
-Results* balancedthickness2_core(FemModel* model);
-Results* balancedvelocities_core(FemModel* model);
-Results* slopecompute_core(FemModel* model);
-Results* control_core(FemModel* model);
-Results* steadystate_core(FemModel* model);
-Results* transient_core(FemModel* model);
-Results* transient_core_2d(FemModel* model);
-Results* transient_core_3d(FemModel* model);
-void adjoint_core(FemModel* model);
-void gradient_core(FemModel* model,int step=0, double search_scalar=0);
-void diagnostic_core(FemModel* model);
-void thermal_core(FemModel* model);
-void surfaceslope_core(FemModel* femmodel);
-void bedslope_core(FemModel* femmodel);
+Results* prognostic_core(FemModel* femmodel);
+Results* prognostic2_core(FemModel* femmodel);
+Results* balancedthickness_core(FemModel* femmodel);
+Results* balancedthickness2_core(FemModel* femmodel);
+Results* balancedvelocities_core(FemModel* femmodel);
+Results* slopecompute_core(FemModel* femmodel);
+Results* steadystate_core(FemModel* femmodel);
+Results* transient_core(FemModel* femmodel);
+Results* transient_core_2d(FemModel* femmodel);
+Results* transient_core_3d(FemModel* femmodel);
+void adjoint_core(FemModel* femmodel);
+void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0);
+void diagnostic_core(FemModel* femmodel);
+void thermal_core(FemModel* femmodel);
+void surfaceslope_core(FemModel* femfemmodel);
+void bedslope_core(FemModel* femfemmodel);
+void control_core(FemModel* femmodel);
 
 //int GradJOrth(WorkspaceParams* workspaceparams);
@@ -35,11 +35,11 @@
 int controlconvergence(double* J, double* fit, double eps_cm, int n);
 
-int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
+int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel);
 
-int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
+int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel);
 	
 double objectivefunctionC(double search_scalar,OptArgs* optargs);
 
-int GradJSearch(double* search_vector,FemModel* femmodel,int step);
+int GradJSearch(double* search_vector,FemModel* femfemmodel,int step);
 //int GradJCheck(WorkspaceParams* workspaceparams,int step,int status);
 
@@ -47,10 +47,10 @@
 void WriteLockFile(char* filename);
 
-void ControlInitialization(FemModel* model);
-void ControlRestart(FemModel* femmodel);
+void stokescontrolinit(FemModel* femfemmodel);
+void controlrestart(FemModel* femfemmodel);
 
-void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
-//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
-void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type);
+void CreateFemModel(FemModel* femfemmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
+//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femfemmodel,char* filename);
+void ResetBoundaryConditions(FemModel* femfemmodel, int analysis_type, int sub_analysis_type);
 
 #endif
Index: /issm/trunk/src/c/solutions/stokescontrolinit.cpp
===================================================================
--- /issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4052)
+++ /issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4052)
@@ -0,0 +1,55 @@
+/*!\file: stokescontrolinit.cpp
+ * \brief: ...
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../include/include.h"
+#include "../solvers/solvers.h"
+
+void stokescontrolinit(FemModel* femmodel){
+
+	/*flags: */
+	int verbose=0;
+	int isstokes=0;
+	double stokesreconditioning;
+
+	/*first recover parameters common to all solutions:*/
+	model->FindParam(&verbose,VerboseEnum);
+	model->FindParam(&isstokes,IsStokesEnum);
+	model->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+
+	/*if no Stokes analysis carried out, assign output and return*/
+	if (!isstokes)femmodel->SetAnalysisType(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum);
+
+	/* 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 model here, and constrain the Stokes model 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: */
+	bedslope_core(femmodel);
+	
+	/*Run a complete diagnostic to update the Stokes spcs: */
+	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution
+
+	//vertical velocity
+	solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum);
+
+	//recondition" pressure computed previously:
+	DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
+	ScaleInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning);
+
+	if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed...");
+	ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
+
+	if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
+	solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
+
+	/*Assign output*/
+	model->SetActiveFormulation(fem_ds);
+}