Index: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 2721)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 2722)
@@ -42,4 +42,5 @@
 int PrognosticAnalysisEnum(void){           return          250; }
 int BalancedthicknessAnalysisEnum(void){    return          251; }
+int BalancedvelocitiesAnalysisEnum(void){   return          252; }
 //melting
 int MeltingAnalysisEnum(void){              return          260; }
Index: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 2721)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 2722)
@@ -43,4 +43,5 @@
 int PrognosticAnalysisEnum(void);
 int BalancedthicknessAnalysisEnum(void);
+int BalancedvelocitiesAnalysisEnum(void);
 //melting
 int MeltingAnalysisEnum(void);
Index: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am	(revision 2721)
+++ issm/trunk/src/c/Makefile.am	(revision 2722)
@@ -222,4 +222,8 @@
 					./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
 					./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\
 					./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
 					./Dofx/Dofx.h\
@@ -525,4 +529,8 @@
 					./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
 					./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\
+					./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\
 					./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
 					./Dofx/Dofx.h\
@@ -625,4 +633,5 @@
 					./parallel/prognostic_core.cpp\
 					./parallel/balancedthickness_core.cpp\
+					./parallel/balancedvelocities_core.cpp\
 					./parallel/transient_core.cpp\
 					./parallel/transient_core_2d.cpp\
@@ -636,5 +645,5 @@
 bin_PROGRAMS = 
 else 
-bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe transient.exe steadystate.exe
+bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe
 endif
 
@@ -656,4 +665,7 @@
 balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
 
+balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp
+balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
+
 transient_exe_SOURCES = parallel/transient.cpp
 transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
Index: issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp	(revision 2722)
+++ issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp	(revision 2722)
@@ -0,0 +1,77 @@
+/*
+ * CreateConstraintsBalancedvelocities.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateConstraintsBalancedvelocities"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+
+void	CreateConstraintsBalancedvelocities(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+
+	int i;
+	int count;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*spc intermediary data: */
+	int spc_sid;
+	int spc_node;
+	int spc_dof;
+	double spc_value;
+	
+	double* spcvelocity=NULL;
+	
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum());
+
+	/*Fetch data: */
+	IoModelFetchData(&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
+
+	count=0;
+
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((iomodel->my_grids[i]==1)){
+	#endif
+
+		if ((int)spcvelocity[6*i+0] && (int)spcvelocity[6*i+1]){ //spc if vx and vy are constrained
+	
+			/*This grid needs to be spc'd: */
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=1; //we enforce first translation degree of freedom, for temperature
+			spc_value=pow( pow(*(spcvelocity+6*i+4),2.0) + pow(*(spcvelocity+6*i+5),2.0) ,0.5);
+
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+		}
+
+	#ifdef _PARALLEL_
+	} //if((my_grids[i]==1))
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	constraints->Presort();
+
+	/*Free data: */
+	xfree((void**)&spcvelocity);
+	
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 2722)
+++ issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 2722)
@@ -0,0 +1,541 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedvelocities.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateElementsNodesAndMaterialsBalancedvelocities"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../IoModel.h"
+
+
+void	CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+
+	/*output: int* epart, int* my_grids, double* my_bordergrids*/
+
+
+	int i,j,k,n;
+	extern int my_rank;
+	extern int num_procs;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    materials = NULL;
+	
+	/*Objects: */
+	Node*       node   = NULL;
+	Tria*       tria = NULL;
+	Penta*      penta = NULL;
+	Matice*     matice  = NULL;
+	Matpar*     matpar  = NULL;
+
+	/*output: */
+	int* epart=NULL; //element partitioning.
+	int* npart=NULL; //node partitioning.
+	int* my_grids=NULL;
+	double* my_bordergrids=NULL;
+
+
+	/*intermediary: */
+	int elements_width; //size of elements
+	double B_avg;
+			
+	/*tria constructor input: */
+	int tria_id;
+	int tria_mid;
+	int tria_mparid;
+	int tria_numparid;
+	int tria_g[3];
+	double tria_h[3];
+	double tria_s[3];
+	double tria_b[3];
+	double tria_k[3];
+	double tria_melting[3];
+	double tria_accumulation[3];
+	double tria_geothermalflux[3];
+	int    tria_friction_type;
+	double tria_p;
+	double tria_q;
+	int    tria_shelf;
+	bool   tria_onwater; 
+
+	/*matice constructor input: */
+	int    matice_mid;
+	double matice_B;
+	double matice_n;
+	
+	/*penta constructor input: */
+
+	int penta_id;
+	int penta_mid;
+	int penta_mparid;
+	int penta_numparid;
+	int penta_g[6];
+	double penta_h[6];
+	double penta_s[6];
+	double penta_b[6];
+	double penta_k[6];
+	int penta_friction_type;
+	double penta_p;
+	double penta_q;
+	int penta_shelf;
+	int penta_onbed;
+	int penta_onsurface;
+	int penta_collapse;
+	double penta_melting[6];
+	double penta_accumulation[6];
+	double penta_geothermalflux[6];
+	int penta_thermal_steadystate;
+	bool   penta_onwater;
+
+	/*matpar constructor input: */
+	int	matpar_mid;
+	double matpar_rho_ice; 
+	double matpar_rho_water;
+	double matpar_heatcapacity;
+	double matpar_thermalconductivity;
+	double matpar_latentheat;
+	double matpar_beta;
+	double matpar_meltingpoint;
+	double matpar_mixed_layer_capacity;
+	double matpar_thermal_exchange_velocity;
+	double matpar_g;
+
+	/* node constructor input: */
+	int node_id;
+	int node_partitionborder=0;
+	double node_x[3];
+	double node_sigma;
+	int node_onbed;
+	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
+	int node_upper_node_id;
+	int node_numdofs;
+
+		
+	#ifdef _PARALLEL_
+	/*Metis partitioning: */
+	int  range;
+	Vec  gridborder=NULL;
+	int  my_numgrids;
+	int* all_numgrids=NULL;
+	int  gridcount;
+	int  count;
+	#endif
+	int  first_grid_index;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum());
+	nodes     = new DataSet(NodesEnum());
+	materials = new DataSet(MaterialsEnum());
+
+	/*Width of elements: */
+	if(strcmp(iomodel->meshtype,"2d")==0){
+		elements_width=3; //tria elements
+	}
+	else{
+		elements_width=6; //penta elements
+	}
+
+	#ifdef _PARALLEL_
+	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
+	if(strcmp(iomodel->meshtype,"2d")==0){
+		/*load elements: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	}
+	else{
+		/*load elements2d: */
+		IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
+	}
+
+
+	MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
+
+	/*Free elements and elements2d: */
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->elements2d);
+
+	/*Used later on: */
+	my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
+	#endif
+
+
+
+	/*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
+
+	/*2d mesh: */
+	if (strcmp(iomodel->meshtype,"2d")==0){
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+		#ifdef _PARALLEL_
+		/*!All elements have been partitioned above, only create elements for this CPU: */
+		if(my_rank==epart[i]){ 
+		#endif
+			
+			
+			/*ids: */
+			tria_id=i+1; //matlab indexing.
+			tria_mid=-1; //no need for materials
+			tria_mparid=-1; //no need for materials
+			tria_numparid=1;
+
+			/*vertices offsets: */
+			tria_g[0]=(int)*(iomodel->elements+elements_width*i+0);
+			tria_g[1]=(int)*(iomodel->elements+elements_width*i+1);
+			tria_g[2]=(int)*(iomodel->elements+elements_width*i+2);
+
+			/*thickness,surface and bed:*/
+			tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
+			tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+			tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
+
+			tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1)); 
+			tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+			tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1)); 
+
+			tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
+			tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+			tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
+
+			/*element on iceshelf?:*/
+			tria_shelf=(int)*(iomodel->elementoniceshelf+i);
+			tria_onwater=(bool)*(iomodel->elementonwater+i);
+
+			/*Create tria element using its constructor:*/
+			tria=new Tria(tria_id, tria_mid, tria_mparid, tria_numparid,tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_onwater);
+
+			/*Add tria element to elements dataset: */
+			elements->AddObject(tria);
+
+			#ifdef _PARALLEL_
+			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
+			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
+			 will hold which grids belong to this partition*/
+			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
+			#endif
+
+		#ifdef _PARALLEL_
+		}//if(my_rank==epart[i])
+		#endif
+
+		}//for (i=0;i<numberofelements;i++)
+
+	
+		/*Free data : */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
+		IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
+		IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	
+		for (i=0;i<iomodel->numberofelements;i++){
+		#ifdef _PARALLEL_
+		/*We are using our element partition to decide which elements will be created on this node: */
+		if(my_rank==epart[i]){
+		#endif
+
+			
+			/*name and id: */
+			penta_id=i+1; //matlab indexing.
+			penta_mid=-1; 
+			penta_mparid=-1; //no need for materials
+			penta_numparid=1;
+
+			/*vertices,thickness,surface,bed and drag: */
+			for(j=0;j<6;j++){
+				penta_g[j]=(int)*(iomodel->elements+elements_width*i+j);
+				penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+				penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+				penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+			}
+
+			/*diverse: */
+			penta_shelf=(int)*(iomodel->elementoniceshelf+i);
+			penta_onbed=(int)*(iomodel->elementonbed+i);
+			penta_onsurface=(int)*(iomodel->elementonsurface+i);
+			penta_collapse=1;
+			penta_onwater=(bool)*(iomodel->elementonwater+i);
+	
+
+			/*Create Penta using its constructor:*/
+			penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
+					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,
+					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,
+					penta_thermal_steadystate,penta_onwater); 
+
+			/*Add penta element to elements dataset: */
+			elements->AddObject(penta);
+	
+			#ifdef _PARALLEL_
+			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
+			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
+			 will hold which grids belong to this partition*/
+			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
+			#endif
+
+		#ifdef _PARALLEL_
+		}//if(my_rank==epart[i])
+		#endif
+
+		}//for (i=0;i<numberofelements;i++)
+
+		/*Free data: */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonbed);
+		xfree((void**)&iomodel->elementonsurface);
+		xfree((void**)&iomodel->elementonwater);
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	#ifdef _PARALLEL_
+		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
+		 *element partition, and which are on its border with other nodes:*/
+		gridborder=NewVec(iomodel->numberofnodes);
+
+		for (i=0;i<iomodel->numberofnodes;i++){
+			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
+		}
+		VecAssemblyBegin(gridborder);
+		VecAssemblyEnd(gridborder);
+
+		#ifdef _ISSM_DEBUG_
+		VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
+		#endif
+		
+		VecToMPISerial(&my_bordergrids,gridborder);
+
+		#ifdef _ISSM_DEBUG_
+		if(my_rank==0){
+			for (i=0;i<iomodel->numberofnodes;i++){
+				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
+			}
+		}
+		#endif
+	#endif
+
+	/*Add one constant material property to materials: */
+	matpar_mid=1; //put it at the end of the materials
+	matpar_g=iomodel->g; 
+	matpar_rho_ice=iomodel->rho_ice; 
+	matpar_rho_water=iomodel->rho_water; 
+	matpar_thermalconductivity=iomodel->thermalconductivity; 
+	matpar_heatcapacity=iomodel->heatcapacity; 
+	matpar_latentheat=iomodel->latentheat; 
+	matpar_beta=iomodel->beta; 
+	matpar_meltingpoint=iomodel->meltingpoint; 
+	matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 
+	matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 
+
+	/*Create matpar object using its constructor: */
+	matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
+			matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
+			matpar_thermal_exchange_velocity,matpar_g);
+		
+	/*Add to materials datset: */
+	materials->AddObject(matpar);
+
+	/*Partition penalties in 3d: */
+	if(strcmp(iomodel->meshtype,"3d")==0){
+	
+		/*Get penalties: */
+		IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
+
+		if(iomodel->numpenalties){
+
+			iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
+			#ifdef _SERIAL_
+			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
+			#else
+			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
+
+			for(i=0;i<iomodel->numpenalties;i++){
+				first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
+				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
+					/*All grids that are being penalised belong to this node's internal grid partition.:*/
+					iomodel->penaltypartitioning[i]=1;
+				}
+				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
+					iomodel->penaltypartitioning[i]=0;
+				}
+			}
+			#endif
+		}
+
+		/*Free penalties: */
+		xfree((void**)&iomodel->penalties);
+	}
+
+	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
+	 We can therefore determine  which grids are internal to this node's partition 
+	 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
+	 that, go and create the grids*/
+
+	/*Create nodes from x,y,z, as well as the spc values on those grids: */
+		
+	/*First fetch data: */
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+	IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
+	IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
+	IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+	IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
+	IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+
+
+	/*Get number of dofs per node: */
+	DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
+
+	for (i=0;i<iomodel->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((my_grids[i]==1)){
+	#endif
+
+		node_id=i+1; //matlab indexing
+			
+		
+		
+		#ifdef _PARALLEL_
+		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+			node_partitionborder=1;
+		}
+		else{
+			node_partitionborder=0;
+		}
+		#else
+			node_partitionborder=0;
+		#endif
+
+		node_x[0]=iomodel->x[i];
+		node_x[1]=iomodel->y[i];
+		node_x[2]=iomodel->z[i];
+		node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
+		
+		node_onbed=(int)iomodel->gridonbed[i];
+		node_onsurface=(int)iomodel->gridonsurface[i];	
+		node_onshelf=(int)iomodel->gridoniceshelf[i];	
+		node_onsheet=(int)iomodel->gridonicesheet[i];	
+
+		if (strcmp(iomodel->meshtype,"3d")==0){
+			if (isnan(iomodel->uppernodes[i])){
+				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+			}
+			else{
+				node_upper_node_id=(int)iomodel->uppernodes[i];
+			}
+		}
+		else{
+			/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
+			node_upper_node_id=node_id;
+		}
+
+		/*Create node using its constructor: */
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+
+		/*set single point constraints.: */
+		if (strcmp(iomodel->meshtype,"3d")==0){
+			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
+			if (!iomodel->gridonbed[i]){
+				for(k=1;k<=node_numdofs;k++){
+					node->FreezeDof(k);
+				}
+			}
+		}
+		/*Add node to nodes dataset: */
+		nodes->AddObject(node);
+
+	#ifdef _PARALLEL_
+	} //if((my_grids[i]==1))
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	elements->Presort();
+	nodes->Presort();
+	materials->Presort();
+
+	/*Clean fetched data: */
+	xfree((void**)&iomodel->deadgrids);
+	xfree((void**)&iomodel->x);
+	xfree((void**)&iomodel->y);
+	xfree((void**)&iomodel->z);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->gridonbed);
+	xfree((void**)&iomodel->gridonsurface);
+	xfree((void**)&iomodel->uppernodes);
+	xfree((void**)&iomodel->gridonicesheet);
+	xfree((void**)&iomodel->gridoniceshelf);
+	
+
+	/*Keep partitioning information into iomodel*/
+	iomodel->epart=epart;
+	iomodel->my_grids=my_grids;
+	iomodel->my_bordergrids=my_bordergrids;
+
+	/*Free ressources:*/
+	#ifdef _PARALLEL_
+	xfree((void**)&all_numgrids);
+	xfree((void**)&npart);
+	VecFree(&gridborder);
+	#endif
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pmaterials=materials;
+
+}
Index: issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp	(revision 2722)
+++ issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp	(revision 2722)
@@ -0,0 +1,35 @@
+/*! \file CreateLoadsBalancedvelocities.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateLoadsBalancedvelocities"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../IoModel.h"
+
+
+void	CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum());
+	
+	
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	loads->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
+
+
Index: issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp	(revision 2722)
+++ issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp	(revision 2722)
@@ -0,0 +1,103 @@
+/*!\file: CreateParametersBalancedvelocities.cpp
+ * \brief driver for creating parameters dataset, for prognostic analysis.
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateParametersBalancedvelocities"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+void CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	Param*   param = NULL;
+	DataSet* parameters=NULL;
+	int      count;
+	int      i;
+	int      dim;
+
+	double* vx=NULL;
+	double* vy=NULL;
+	double* vz=NULL;
+	double* u_g=NULL;
+	double* pressure=NULL;
+	double* temperature=NULL;
+	double* thickness=NULL;
+	double* surface=NULL;
+	double* bed=NULL;
+	double* accumulation=NULL;
+	double* melting=NULL;
+
+	/*recover parameters : */
+	parameters=*pparameters;
+
+	count=parameters->Size();
+
+	/*Get vx and vy: */
+	IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz");
+
+	u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
+
+	if(vx) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
+	if(vy) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
+	if(vz) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
+
+	count++;
+	param= new Param(count,"u_g",DOUBLEVEC);
+	param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
+	parameters->AddObject(param);
+
+
+	xfree((void**)&vx);
+	xfree((void**)&vy);
+	xfree((void**)&vz);
+	xfree((void**)&u_g);
+
+	/*Get thickness: */
+	IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness");
+	
+	count++;
+	param= new Param(count,"h_g",DOUBLEVEC);
+	if(thickness) param->SetDoubleVec(thickness,iomodel->numberofnodes,1);
+	else param->SetDoubleVec(thickness,0,0);
+	parameters->AddObject(param);
+
+	/*Free thickness: */
+	xfree((void**)&thickness);
+
+	/*Get melting: */
+	IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting");
+	if(melting) for(i=0;i<iomodel->numberofnodes;i++)melting[i]=melting[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"m_g",DOUBLEVEC);
+	if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1);
+	else param->SetDoubleVec(melting,0,1);
+	parameters->AddObject(param);
+
+	/*Free melting: */
+	xfree((void**)&melting);
+
+	/*Get accumulation: */
+	IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
+	if(accumulation) for(i=0;i<iomodel->numberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"a_g",DOUBLEVEC);
+	if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1);
+	else param->SetDoubleVec(accumulation,0,0);
+	parameters->AddObject(param);
+
+	/*Free accumulation: */
+	xfree((void**)&accumulation);
+
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2721)
+++ issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2722)
@@ -94,4 +94,12 @@
 
 	}
+	else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum()){
+
+		CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle);
+		CreateParametersBalancedvelocities(pparameters,iomodel,iomodel_handle);
+
+	}
 	else{
 
Index: issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2721)
+++ issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2722)
@@ -177,23 +177,22 @@
 	int* my_grids; /*! grids that belong to this cpu*/
 	double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/
-	int*     penaltypartitioning;
+	int*    penaltypartitioning;
 
 };
-
 
 
 	/*constructor and destructor: */
 	IoModel*	NewIoModel(void);
-	void    DeleteIoModel( IoModel** pthis);
+	void     DeleteIoModel( IoModel** pthis);
 
 	/*Echo: */
-	void    IoModelEcho(IoModel* iomodel,int which_part,int rank);
+	void  IoModelEcho(IoModel* iomodel,int which_part,int rank);
 
 	/*Initialization with matlab workspace data, or marshall binary data: */
-	int	    IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle);
+	int	IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle);
 
 	/*Creation of fem datasets: general drivers*/
-	void    CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	
@@ -203,57 +202,62 @@
 	void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
-	void    CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*diagnostic vertical*/
 	void	CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
 	/*diagnostic hutter*/
 	void	CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
 	/*diagnostic stokes*/
 	void	CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
 	/*slope compute*/
 	void	CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
 
 	/*control:*/
-	void    CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*thermal:*/
 	void	CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
-	void    CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*melting:*/
 	void	CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
-	void    CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*prognostic:*/
 	void	CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
-	void    CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*balancedthickness:*/
 	void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
 	void	CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
-	void    CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
-	void    CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+
+	/*balancedvelocities:*/
+	void	CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void	CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
 	/*qmu: */
 	void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
-
 #endif  /* _IOMODEL_H */
Index: issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Penta.cpp	(revision 2721)
+++ issm/trunk/src/c/objects/Penta.cpp	(revision 2722)
@@ -325,4 +325,8 @@
 		CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
+	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
+
+		CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
+	}
 	else if (analysis_type==ThermalAnalysisEnum()){
 
@@ -344,4 +348,27 @@
 
 void  Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
+
+	/*Collapsed formulation: */
+	Tria*  tria=NULL;
+
+	/*If on water, skip: */
+	if(onwater)return;
+
+	/*Is this element on the bed? :*/
+	if(!onbed)return;
+
+	/*Spawn Tria element from the base of the Penta: */
+	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+	tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
+	delete tria;
+	return;
+
+}
+/*}}}*/
+/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreateKMatrixBalancedvelocities"
+
+void  Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Collapsed formulation: */
@@ -1410,4 +1437,8 @@
 		CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
 	}
+	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
+
+		CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type);
+	}
 	else if (analysis_type==ThermalAnalysisEnum()){
 
@@ -1429,4 +1460,26 @@
 
 void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
+
+	/*Collapsed formulation: */
+	Tria*  tria=NULL;
+
+	/*If on water, skip: */
+	if(onwater)return;
+
+	/*Is this element on the bed? :*/
+	if(!onbed)return;
+
+	/*Spawn Tria element from the base of the Penta: */
+	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+	tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
+	delete tria;
+	return;
+}
+/*}}}*/
+/*FUNCTION CreatePVectorBalancedvelocities {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::CreatePVectorBalancedvelocities"
+
+void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*Collapsed formulation: */
Index: issm/trunk/src/c/objects/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Penta.h	(revision 2721)
+++ issm/trunk/src/c/objects/Penta.h	(revision 2722)
@@ -118,4 +118,6 @@
 		void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
 
 		void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);
Index: issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Tria.cpp	(revision 2721)
+++ issm/trunk/src/c/objects/Tria.cpp	(revision 2722)
@@ -303,4 +303,9 @@
 
 	}
+	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
+
+		CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
+
+	}
 	else{
 		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
@@ -447,4 +452,216 @@
 		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
 		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
+
+		if(numpar->artdiff){
+
+			/* Compute artificial diffusivity */
+			KDL[0][0]=DL_scalar*K[0][0];
+			KDL[1][1]=DL_scalar*K[1][1];
+
+			TripleMultiply( &Bprime[0][0],2,numdof,1,
+						&KDL[0][0],2,2,0,
+						&Bprime[0][0],2,numdof,0,
+						&Ke_gg_gaussian[0][0],0);
+
+			/* Add artificial diffusivity matrix */
+			for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+
+		}
+
+#ifdef _DEBUGELEMENTS_
+		if(my_rank==RANK && id==ELID){ 
+			printf("      B:\n");
+			for(i=0;i<3;i++){
+				for(j=0;j<numdof;j++){
+					printf("%g ",B[i][j]);
+				}
+				printf("\n");
+			}
+			printf("      Bprime:\n");
+			for(i=0;i<3;i++){
+				for(j=0;j<numdof;j++){
+					printf("%g ",Bprime[i][j]);
+				}
+				printf("\n");
+			}
+		}
+#endif
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+#ifdef _DEBUGELEMENTS_
+	if(my_rank==RANK && id==ELID){ 
+		printf("      Ke_gg erms:\n");
+		for( i=0; i<numdof; i++){
+			for (j=0;j<numdof;j++){
+				printf("%g ",Ke_gg[i][j]);
+			}
+			printf("\n");
+		}
+		printf("      Ke_gg row_indices:\n");
+		for( i=0; i<numdof; i++){
+			printf("%i ",doflist[i]);
+		}
+
+	}
+#endif
+
+cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
+/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixBalancedvelocities"
+void  Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdof=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrices: */
+	double L[numgrids];
+	double B[2][numgrids];
+	double Bprime[2][numgrids];
+	double DL[2][2]={0.0};
+	double DLprime[2][2]={0.0};
+	double DL_scalar;
+	double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 
+	double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Jdettria;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double  surface_normal[3];
+	double  nx,ny,norm;
+	double  vxvy_list[numgrids][2]={0.0};
+	double  vx_list[numgrids]={0.0};
+	double  vy_list[numgrids]={0.0};
+	double  dvx[2];
+	double  dvy[2];
+	double  vx,vy;
+	double  dvxdx,dvydy;
+	double  v_gauss[2]={0.0};
+	double  K[2][2]={0.0};
+	double  KDL[2][2]={0.0};
+	int     dofs[2]={0,1};
+	int     found=0;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
+
+	for(i=0;i<numgrids;i++){
+		vx_list[i]=vxvy_list[i][0];
+		vy_list[i]=vxvy_list[i][1];
+	}
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/*Modify z so that it reflects the surface*/
+	for(i=0;i<numgrids;i++) xyz_list[i][2]=s[i];
+
+	/*Get normal vector to the surface*/
+	nx=(vx_list[0]+vx_list[1]+vx_list[2])/3;
+	ny=(vy_list[0]+vy_list[1]+vy_list[2])/3;
+	if(nx==0 && ny==0){
+		SurfaceNormal(&surface_normal[0],xyz_list);
+		nx=surface_normal[0];
+		ny=surface_normal[1];
+	}
+	if(nx==0 && ny==0){
+		nx=0;
+		ny=1;
+	}
+	norm=pow( pow(nx,2)+pow(ny,2) , (double).5);
+	nx=nx/norm;
+	ny=ny/norm;
+
+	//Create Artificial diffusivity once for all if requested
+	if(numpar->artdiff){
+		//Get the Jacobian determinant
+		gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		//Build K matrix (artificial diffusivity matrix)
+		v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
+		v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
+
+		K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
+		K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get B  and B prime matrix: */
+		GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+
+		//Get vx, vy and their derivatives at gauss point
+		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
+		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+
+		GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
+		GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
+
+		dvxdx=dvx[0];
+		dvydy=dvy[1];
+
+		DL_scalar=gauss_weight*Jdettria;
+
+		DLprime[0][0]=DL_scalar*nx;
+		DLprime[1][1]=DL_scalar*ny;
+
+		//Do the triple product tL*D*L. 
+		//Ke_gg_velocities=B'*DLprime*Bprime;
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DLprime[0][0],2,2,0,
+					&Bprime[0][0],2,numdof,0,
+					&Ke_gg_velocities2[0][0],0);
+
+		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
 
 		if(numpar->artdiff){
@@ -1550,4 +1767,8 @@
 		CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
 	}
+	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
+
+		CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type);
+	}
 	else{
 		throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
@@ -1630,4 +1851,95 @@
 
 		/* Get accumulation, melting and thickness at gauss point */
+		GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
+		GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
+
+		/* Add value into pe_g: */
+		for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
+
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add pe_g to global matrix Kgg: */
+	VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
+/*FUNCTION CreatePVectorBalancedvelocities {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorBalancedvelocities"
+void  Tria::CreatePVectorBalancedvelocities(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
+
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdof=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrix */
+	double pe_g[numgrids]={0.0};
+	double L[numgrids];
+	double Jdettria;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double  accumulation_list[numgrids]={0.0};
+	double  accumulation_g;
+	double  melting_list[numgrids]={0.0};
+	double  melting_g;
+	int     dofs[1]={0};
+	int     found=0;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!");
+	found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!");
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		/* Get accumulation, melting at gauss point */
 		GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
 		GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
Index: issm/trunk/src/c/objects/Tria.h
===================================================================
--- issm/trunk/src/c/objects/Tria.h	(revision 2721)
+++ issm/trunk/src/c/objects/Tria.h	(revision 2722)
@@ -109,6 +109,6 @@
 		int   GetOnBed();
 
-		void          GetThicknessList(double* thickness_list);
-		void          GetBedList(double* bed_list);
+		void  GetThicknessList(double* thickness_list);
+		void  GetBedList(double* bed_list);
 		Object* copy();
 		void  NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets);
@@ -125,4 +125,6 @@
 		void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorBalancedvelocities(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
 		double MassFlux(double* segment,double* ug);
 		double GetArea(void);
Index: issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 2721)
+++ issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 2722)
@@ -90,4 +90,7 @@
 	double* m_g_serial=NULL;
 	double* melting=NULL;
+
+	Vec     v_g=NULL;
+	double* v_g_serial=NULL;
 
 	int numberofnodes;
@@ -117,4 +120,7 @@
 	if(analysis_type==BalancedthicknessAnalysisEnum()){
 		fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum());
+	}
+	if(analysis_type==BalancedvelocitiesAnalysisEnum()){
+		fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum());
 	}
 	fem_t=model->GetFormulation(ThermalAnalysisEnum());
@@ -327,4 +333,27 @@
 			VecFree(&h_g);
 		}
+		else if(strcmp(result->GetFieldName(),"v_g")==0){
+			/*easy, v_g is of size numberofnodes, on 1 dof, just repartition: */
+			result->GetField(&v_g);
+			VecToMPISerial(&v_g_serial,v_g);
+			fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
+			VecToMPISerial(&partition,fem_p->partition->vector);
+
+			vel=(double*)xmalloc(numberofnodes*sizeof(double));
+
+			for(i=0;i<numberofnodes;i++){
+				vel[i]=v_g_serial[(int)partition[i]];
+			}
+
+			/*Ok, add pressure to newresults: */
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
+			newresults->AddObject(newresult);
+
+			/*do some cleanup: */
+			xfree((void**)&v_g_serial);
+			xfree((void**)&vel);
+			xfree((void**)&partition);
+			VecFree(&v_g);
+		}
 		else if(strcmp(result->GetFieldName(),"s_g")==0){
 			/*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
Index: issm/trunk/src/c/parallel/balancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/parallel/balancedvelocities.cpp	(revision 2722)
+++ issm/trunk/src/c/parallel/balancedvelocities.cpp	(revision 2722)
@@ -0,0 +1,157 @@
+/*!\file:  balancedvelocities.cpp
+ * \brief: balancedvelocities solution
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "balancedvelocities"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+int main(int argc,char* *argv){
+	
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	int   numberofnodes;
+	double waitonlock=0;
+
+	Model* model=NULL;
+
+	double* u_g_serial=NULL;
+	double* melting_g=NULL;
+	double* accumulation_g=NULL;
+	int     qmu_analysis;
+
+	/*Results: */
+	DataSet* results=NULL;
+	DataSet* processedresults=NULL;
+	Result*  result=NULL;
+
+	ParameterInputs* inputs=NULL;
+	Param*   param=NULL;
+
+	/*time*/
+	double   start, finish;
+	double   start_core, finish_core;
+	double   start_init, finish_init;
+
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	/*Initialize Petsc and get start time*/
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Initialize model structure: */
+	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
+	model=new Model();
+
+	/*Open handle to data on disk: */
+	fid=pfopen(inputfilename,"rb");
+
+	/*Initialize model structure: */
+	model=new Model();
+
+	_printf_("read and create finite element model:\n");
+	model->AddFormulation(fid,BalancedvelocitiesAnalysisEnum());
+
+	/*recover parameters: */
+	model->FindParam(&waitonlock,"waitonlock");
+	model->FindParam(&qmu_analysis,"qmu_analysis");
+
+	_printf_("initialize inputs:\n");
+	
+	model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedvelocitiesAnalysisEnum());
+	model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedvelocitiesAnalysisEnum());
+	model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedvelocitiesAnalysisEnum());
+	model->FindParam(&numberofnodes,"numberofnodes");
+	
+	inputs=new ParameterInputs;
+	inputs->Add("velocity",u_g_serial,3,numberofnodes);
+	inputs->Add("melting",melting_g,1,numberofnodes);
+	inputs->Add("accumulation",accumulation_g,1,numberofnodes);
+
+	_printf_("initialize results:\n");
+	results=new DataSet(ResultsEnum());
+	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
+
+	/*are we running the solutoin sequence, or a qmu wrapper around it? : */
+	if(!qmu_analysis){
+
+		/*run balancedvelocities analysis: */
+		_printf_("call computational core:\n");
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		balancedvelocities_core(results,model,inputs);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+	}
+	else{
+
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on balancedvelocities core:\n");
+	
+		#ifdef _HAVE_DAKOTA_ 
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		Qmux(model,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum());
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	 	#else
+		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	/*Add analysis_type to results: */
+	result=new Result(results->Size()+1,0,1,"analysis_type","balancedvelocities");
+	results->AddObject(result);
+	
+	_printf_("process results:\n");
+	ProcessResults(&processedresults,results,model,BalancedvelocitiesAnalysisEnum());
+	
+	_printf_("write results to disk:\n");
+	OutputResults(processedresults,outputfilename);
+
+	if (waitonlock>0){
+		_printf_("write lock file:\n");
+		WriteLockFile(lockname);
+	}
+
+	/*Free ressources:*/
+	delete processedresults;
+	delete results;
+	delete model;
+	delete inputs;
+
+	/*Get finish time and close*/
+	MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
+	_printf_("\n   %-34s %f seconds  \n","Model initialization elapsed time:",finish_init-start_init);
+	_printf_("   %-34s %f seconds  \n","Core solution elapsed time:",finish_core-start_core);
+	_printf_("   %-34s %f seconds\n\n","Total elapsed time:",finish-start);
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+
+	/*end module: */
+	MODULEEND();
+
+	return 0; //unix success return;
+}
Index: issm/trunk/src/c/parallel/balancedvelocities_core.cpp
===================================================================
--- issm/trunk/src/c/parallel/balancedvelocities_core.cpp	(revision 2722)
+++ issm/trunk/src/c/parallel/balancedvelocities_core.cpp	(revision 2722)
@@ -0,0 +1,64 @@
+/*!\file: balancedvelocities_core.cpp
+ * \brief: core of the balancedvelocities solution 
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "balancedvelocities_core"
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./parallel.h"
+#include "../issm.h"
+
+void balancedvelocities_core(DataSet* results,Model* model,ParameterInputs* inputs){
+
+	extern int my_rank;
+
+	/*output: */
+	Result* result=NULL;
+
+	/*intermediary: */
+	Vec u_g=NULL;
+
+	/*solutions: */
+	Vec v_g=NULL;
+
+	/*flags: */
+	int verbose=0;
+	int numberofdofspernode;
+	int numberofnodes;
+	int dofs[2]={1,1};
+
+	/*fem balancedvelocities model: */
+	FemModel* fem_p=NULL;
+
+
+	/*recover fem model: */
+	fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum());
+
+	//first recover parameters common to all solutions
+	model->FindParam(&verbose,"verbose");
+	model->FindParam(&numberofnodes,"numberofnodes");
+	model->FindParam(&numberofdofspernode,"numberofdofspernode");
+
+	_printf_("depth averaging velocity...\n");
+	u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
+	FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
+	inputs->Add("velocity_average",u_g,2,numberofnodes);
+	
+	_printf_("call computational core:\n");
+	diagnostic_core_linear(&v_g,fem_p,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum());
+
+	_printf_("extrude computed thickness on all layers:\n");
+	FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"v_g",v_g);
+	results->AddObject(result);
+
+	/*Free ressources:*/
+	VecFree(&u_g);
+	VecFree(&v_g);
+}
Index: issm/trunk/src/c/parallel/parallel.h
===================================================================
--- issm/trunk/src/c/parallel/parallel.h	(revision 2721)
+++ issm/trunk/src/c/parallel/parallel.h	(revision 2722)
@@ -18,4 +18,5 @@
 void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs);
 void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs);
+void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs);
 void control_core(DataSet* results,Model* model, ParameterInputs* inputs);
 
Index: issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp
===================================================================
--- issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 2721)
+++ issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 2722)
@@ -59,4 +59,7 @@
 		numdofs=1;
 	}
+	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
+		numdofs=1;
+	}
 	else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type,"  not implemented yet!"));
 
Index: issm/trunk/src/m/classes/public/process_solve_options.m
===================================================================
--- issm/trunk/src/m/classes/public/process_solve_options.m	(revision 2721)
+++ issm/trunk/src/m/classes/public/process_solve_options.m	(revision 2722)
@@ -17,9 +17,10 @@
 
 %check solution type is supported
-if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient','balancedthickness'}),
+if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient',...
+		'balancedthickness','balancedvelocities'}),
 	error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
 else
 	%convert to enum
-	outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum']);
+	outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum()']);
 end
 if ~ismemberi(sub_analysis_type,{'steady','transient','none','horiz','adjoint','gradient','inverse','vert',''}),
@@ -27,5 +28,5 @@
 else
 	%convert to enum
-	outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum']);
+	outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum()']);
 end
 
Index: issm/trunk/src/m/classes/public/solve.m
===================================================================
--- issm/trunk/src/m/classes/public/solve.m	(revision 2721)
+++ issm/trunk/src/m/classes/public/solve.m	(revision 2722)
@@ -73,4 +73,8 @@
 elseif md.analysis_type==BalancedthicknessAnalysisEnum,
 	md=balancedthickness(md);
+
+elseif md.analysis_type==BalancedvelocitiesAnalysisEnum,
+	md=balancedvelocities(md);
+
 else
 	error('solution type not supported for this model configuration.');
Index: issm/trunk/src/m/enum/AnalysisTypeFromEnum.m
===================================================================
--- issm/trunk/src/m/enum/AnalysisTypeFromEnum.m	(revision 2721)
+++ issm/trunk/src/m/enum/AnalysisTypeFromEnum.m	(revision 2722)
@@ -89,4 +89,8 @@
 end
 
+if enum==BalancedvelocitiesAnalysisEnum(),
+	string='balancedvelocities';
+end
+
 if enum==MeltingAnalysisEnum(),
 	string='melting';
Index: issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m
===================================================================
--- issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m	(revision 2722)
+++ issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m	(revision 2722)
@@ -0,0 +1,9 @@
+function macro=BalancedvelocitiesAnalysisEnum()
+%BALANCEDVELOCITIESANALYSISENUM - Enum of BalancedvelocitiesAnalysis
+%
+%   file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum
+%
+%   Usage:
+%      macro=BalancedvelocitiesAnalysisEnum()
+
+macro=252;
Index: issm/trunk/src/m/solutions/jpl/balancedvelocities.m
===================================================================
--- issm/trunk/src/m/solutions/jpl/balancedvelocities.m	(revision 2722)
+++ issm/trunk/src/m/solutions/jpl/balancedvelocities.m	(revision 2722)
@@ -0,0 +1,37 @@
+function md=balancedvelocities(md)
+%BALANCEDVELOCITIES - balancedvelocities solution sequence.
+%
+%   Usage:
+%      md=balancedvelocities(md)
+
+	%timing
+	t1=clock;
+
+	%Build all models requested for diagnostic simulation
+	models.analysis_type=BalancedvelocitiesAnalysisEnum; %needed for processresults
+	
+	displaystring(md.verbose,'%s',['reading balancedvelocities model data']);
+	md.analysis_type=BalancedvelocitiesAnalysisEnum; models.p=CreateFemModel(md);
+
+	% figure out number of dof: just for information purposes.
+	md.dof=modelsize(models);
+
+	%initialize inputs
+	displaystring(md.verbose,'\n%s',['setup inputs...']);
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double');
+
+	displaystring(md.verbose,'\n%s',['call computational core:']);
+	results=balancedvelocities_core(models,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum());
+
+	displaystring(md.verbose,'\n%s',['load results...']);
+	if ~isstruct(md.results), md.results=struct(); end
+	md.results.balancedvelocities=processresults(models, results);
+
+	%stop timing
+	t2=clock;
+	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);	
Index: issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m
===================================================================
--- issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m	(revision 2722)
+++ issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m	(revision 2722)
@@ -0,0 +1,24 @@
+function results=balancedvelocities_core(models,inputs,analysis_type,sub_analysis_type)
+%BALANCEDVELOCITIES_CORE - linear solution sequence
+%
+%   Usage:
+%      v_g=balancedvelocities_core(m,inputs,analysis_type,sub_analysis_type)
+
+	%get FE model
+	m=models.p;
+	results.time=0;
+	results.step=1;
+
+	displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']);
+	%Take only the first two dofs of m.parameters.u_g
+	u_g=get(inputs,'velocity',[1 1 0 0]);
+	u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
+	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
+
+	displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
+	results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
+
+	displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
+	results.v_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0);
+
+end %end function
Index: issm/trunk/src/m/solutions/jpl/processresults.m
===================================================================
--- issm/trunk/src/m/solutions/jpl/processresults.m	(revision 2721)
+++ issm/trunk/src/m/solutions/jpl/processresults.m	(revision 2722)
@@ -126,4 +126,11 @@
 			newresults(i).parameter=param_g;
 
+		elseif strcmpi(resultsname{j},'v_g'),
+
+			v_g=results(i).v_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).vel=v_g;
+
 		else
 
