Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp
===================================================================
--- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3587)
+++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3588)
@@ -18,11 +18,11 @@
 	extern int my_rank;
 	
-	//_printf_("      Configuring elements...\n");
+	_printf_("      Configuring elements...\n");
 	elements->Configure(elements,loads,nodes,vertices,materials,parameters);
-	//_printf_("      Configuring loads...\n");
+	_printf_("      Configuring loads...\n");
 	loads->Configure(elements,loads,nodes,vertices,materials,parameters);
-	//_printf_("      Configuring nodes...\n");
+	_printf_("      Configuring nodes...\n");
 	nodes->Configure(elements,loads,nodes,vertices,materials,parameters);
-	//_printf_("      Configuring parameters...\n");
+	_printf_("      Configuring parameters...\n");
 	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
 
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 3587)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 3588)
@@ -48,4 +48,5 @@
 	Prognostic2AnalysisEnum,
 	BalancedthicknessAnalysisEnum,
+	Balancedthickness2AnalysisEnum,
 	BalancedvelocitiesAnalysisEnum,
 	//melting
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 3587)
+++ /issm/trunk/src/c/Makefile.am	(revision 3588)
@@ -245,4 +245,8 @@
 					./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
 					./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\
 					./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
 					./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
@@ -645,4 +649,8 @@
 					./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
 					./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\
+					./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\
 					./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
 					./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
@@ -761,4 +769,5 @@
 					./parallel/prognostic2_core.cpp\
 					./parallel/balancedthickness_core.cpp\
+					./parallel/balancedthickness2_core.cpp\
 					./parallel/balancedvelocities_core.cpp\
 					./parallel/slopecompute_core.cpp\
@@ -835,5 +844,5 @@
 bin_PROGRAMS = 
 else 
-bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe
+bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedthickness2.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe
 endif
 
@@ -858,4 +867,7 @@
 balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
 
+balancedthickness2_exe_SOURCES = parallel/balancedthickness2.cpp
+balancedthickness2_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
+
 balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp
 balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 3588)
@@ -11,9 +11,36 @@
 
 void	CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
-	
+
+	int i;
+	int count=0;
 	DataSet* constraints = NULL;
-	
+
 	/*Create constraints: */
 	constraints = new DataSet(ConstraintsEnum);
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
+
+	count=1; //matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		if(iomodel->my_vertices[i]){
+			
+			if ((int)iomodel->spcthickness[2*i]){
+		
+				constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));//we enforce first translation degree of freedom, for temperature
+				count++;
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->spcthickness);
+
+	/*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();
+
+	cleanup_and_return:
 	
 	/*Assign output pointer: */
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3588)
@@ -8,21 +8,18 @@
 #include "../../objects/objects.h"
 #include "../../shared/shared.h"
-#include "../../include/macros.h"
 #include "../../include/typedefs.h"
 #include "../../MeshPartitionx/MeshPartitionx.h"
 #include "../IoModel.h"
 
-void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 	/*Intermediary*/
-	int i,j;
-	int vertex_index;
-	int node_index;
+	int i,j,k,n;
 
 	/*DataSets: */
-	DataSet* elements  = NULL;
-	DataSet* nodes = NULL;
-	DataSet* vertices = NULL;
-	DataSet* materials = NULL;
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
 
 	/*First create the elements, nodes and material properties: */
@@ -35,5 +32,4 @@
 	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
 
-	/*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
 	/*2d mesh: */
 	if (strcmp(iomodel->meshtype,"2d")==0){
@@ -58,6 +54,7 @@
 			}
 		}//for (i=0;i<numberofelements;i++)
-	
+
 		/*Free data : */
+		xfree((void**)&iomodel->elements);
 		xfree((void**)&iomodel->thickness);
 		xfree((void**)&iomodel->surface);
@@ -68,11 +65,45 @@
 	}
 	else{ //	if (strcmp(meshtype,"2d")==0)
-		ISSMERROR("not implemented yet");
+
+		/*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++){
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i,iomodel,6));
+			}
+		}//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)
 
-	/*Add new constrant material property tgo materials, at the end: */
+	/*Add new constrant material property to materials, at the end: */
 	materials->AddObject(new Matpar(iomodel));
 
-	/*Create nodes and vertices: */
+	/*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");
@@ -82,13 +113,7 @@
 	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
 	IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
-	IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
 	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
 	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
-	if (strcmp(iomodel->meshtype,"3d")==0){
-		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
-		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
-	}
 
-	/*Build Vertices dataset*/
 	for (i=0;i<iomodel->numberofvertices;i++){
 
@@ -99,24 +124,7 @@
 			vertices->AddObject(new Vertex(i,iomodel));
 
-		}
-	}
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i,iomodel));
 
-	/*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
-	for (i=0;i<iomodel->numberofelements;i++){
-		for (j=0;j<3;j++){
-
-			if(iomodel->my_nodes[3*i+j]){ 
-
-				//Get index of the vertex on which the current node is located
-				vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
-				ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);
-
-				//Compute Node index (id-1)
-				node_index=3*i+j;
-
-				/*Add node to nodes dataset: */
-				nodes->AddObject(new Node(vertex_index,node_index,iomodel));
-
-			}
 		}
 	}
@@ -131,5 +139,4 @@
 	xfree((void**)&iomodel->gridonbed);
 	xfree((void**)&iomodel->gridonsurface);
-	xfree((void**)&iomodel->gridonhutter);
 	xfree((void**)&iomodel->uppernodes);
 	xfree((void**)&iomodel->gridonicesheet);
@@ -139,6 +146,6 @@
 	 * datasets, it will not be redone: */
 	elements->Presort();
+	vertices->Presort();
 	nodes->Presort();
-	vertices->Presort();
 	materials->Presort();
 
@@ -150,3 +157,4 @@
 	*pvertices=vertices;
 	*pmaterials=materials;
+
 }
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 3588)
@@ -8,108 +8,20 @@
 #include "../../shared/shared.h"
 #include "../../include/macros.h"
-#include "../../include/typedefs.h"
 #include "../IoModel.h"
+
 
 void	CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
-	/*Intermediary*/
-	int i,j;
-	int i1,i2;
-	int pos1,pos2;
-	double e1,e2;
-
-	/*Output*/
-	DataSet* loads=NULL;
-
-	/*numericalflux intermediary data: */
-	char numericalflux_type[NUMERICALFLUXSTRING];
-	int  numericalflux_id;
-	int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
-	int  numericalflux_elem_id;
-	double numericalflux_h[MAX_NUMERICALFLUX_NODES];
+	DataSet*    loads    = NULL;
 
 	/*Create loads: */
 	loads   = new DataSet(LoadsEnum);
-
-	/*Get edges and elements*/
-	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
-	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
-	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
-
-	/*First load data:*/
-	for (i=0;i<iomodel->numberofedges;i++){
-
-		/*Get left and right elements*/
-		e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
-		e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
-
-		/*Now, if this element is not in the partition, pass: */
-		if(!iomodel->my_elements[(int)e1]) continue;
-
-		/*Create load*/
-		numericalflux_id=i+1; //Matlab indexing
-		numericalflux_elem_id=(int)e1+1;//id is in matlab index
-
-		/*1: Get vertices ids*/
-		i1=(int)iomodel->edges[4*i+0];
-		i2=(int)iomodel->edges[4*i+1];
-
-		if (!isnan(e2)){
-			strcpy(numericalflux_type,"internal");
-
-			/*Now, we must get the nodes of the 4 nodes located on the edge*/
-
-			/*2: Get the column where these ids are located in the index*/
-			pos1=pos2=UNDEF;
-			for(j=0;j<3;j++){
-				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
-				if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1;
-			}
-			ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
-
-			/*3: We have the id of the elements and the position of the vertices in the index
-			 * we can compute their dofs!*/
-			numericalflux_node_ids[0]=3*(int)e1+pos1;       //ex: 1 2 3
-			numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1
-			numericalflux_node_ids[2]=3*(int)e2+pos2;           //ex: 1 2 3
-			numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2
-
-			numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1];
-			numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
-			numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1];
-			numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1];
-		}
-		else{
-			strcpy(numericalflux_type,"boundary");
-
-			/*2: Get the column where these ids are located in the index*/
-			pos1==UNDEF;
-			for(j=0;j<3;j++){
-				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
-			}
-			ISSMASSERT(pos1!=UNDEF);
-
-			/*3: We have the id of the elements and the position of the vertices in the index
-			 * we can compute their dofs!*/
-			numericalflux_node_ids[0]=3*(int)e1+pos1;
-			numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1;
-
-			numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1];
-			numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
-			numericalflux_h[2]=UNDEF;
-			numericalflux_h[3]=UNDEF;
-		}
-
-		loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h));
-	}
-
-	/*Free data: */
-	xfree((void**)&iomodel->edges);
-	xfree((void**)&iomodel->elements);
-	xfree((void**)&iomodel->thickness);
-
+	
+	
 	/*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: */
@@ -117,2 +29,4 @@
 
 }
+
+
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp	(revision 3588)
@@ -1,4 +1,3 @@
-/*!\file: CreateParametersPrognostic.cpp
- * \brief driver for creating parameters dataset, for prognostic analysis.
+ /* \brief driver for creating parameters dataset, for prognostic analysis.
  */ 
 
@@ -16,7 +15,5 @@
 	int      count;
 	int      i;
-	int      dim;
-	double*  vx_g=NULL;
-	double*  vy_g=NULL;
+	double* u_g=NULL;
 
 	/*recover parameters : */
@@ -28,49 +25,22 @@
 	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
 	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz");
 
-	vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
-	vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
+	u_g=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double));
 
-	if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts;
-	if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts;
+	if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+0]=iomodel->vx[i]/iomodel->yts;
+	if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+1]=iomodel->vy[i]/iomodel->yts;
+	if(iomodel->vz) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+2]=iomodel->vz[i]/iomodel->yts;
 
 	count++;
-	param= new Param(count,"vx_g",DOUBLEVEC);
-	param->SetDoubleVec(vx_g,iomodel->numberofvertices,1);
+	param= new Param(count,"u_g",DOUBLEVEC);
+	param->SetDoubleVec(u_g,3*iomodel->numberofvertices,3);
 	parameters->AddObject(param);
-	count++;
-	param= new Param(count,"vy_g",DOUBLEVEC);
-	param->SetDoubleVec(vy_g,iomodel->numberofvertices,1);
-	parameters->AddObject(param);
+
 
 	xfree((void**)&iomodel->vx);
 	xfree((void**)&iomodel->vy);
-	xfree((void**)&vx_g);
-	xfree((void**)&vy_g);
-
-	/*Get thickness: */
-	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
-
-	count++;
-	param= new Param(count,"h_g",DOUBLEVEC);
-	if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1);
-	else param->SetDoubleVec(iomodel->thickness,0,0);
-	parameters->AddObject(param);
-
-	/*Free thickness: */
-	xfree((void**)&iomodel->thickness);
-
-	/*Get dhdt: */
-	IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
-	if(iomodel->dhdt) for(i=0;i<iomodel->numberofvertices;i++)iomodel->dhdt[i]=iomodel->dhdt[i]/iomodel->yts;
-
-	count++;
-	param= new Param(count,"dhdt_g",DOUBLEVEC);
-	if(iomodel->dhdt) param->SetDoubleVec(iomodel->dhdt,iomodel->numberofvertices,1);
-	else param->SetDoubleVec(iomodel->dhdt,0,1);
-	parameters->AddObject(param);
-
-	/*Free dhdt: */
-	xfree((void**)&iomodel->dhdt);
+	xfree((void**)&iomodel->vz);
+	xfree((void**)&u_g);
 
 	/*Get melting: */
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp	(revision 3588)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp	(revision 3588)
@@ -0,0 +1,21 @@
+/*
+ * CreateConstraintsBalancedthickness2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+void	CreateConstraintsBalancedthickness2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	DataSet* constraints = NULL;
+	
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp	(revision 3588)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp	(revision 3588)
@@ -0,0 +1,152 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedthickness2.c:
+ */
+
+#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 "../../include/typedefs.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../IoModel.h"
+
+void	CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j;
+	int vertex_index;
+	int node_index;
+
+	/*DataSets: */
+	DataSet* elements  = NULL;
+	DataSet* nodes = NULL;
+	DataSet* vertices = NULL;
+	DataSet* materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i,iomodel,3));
+			}
+		}//for (i=0;i<numberofelements;i++)
+	
+		/*Free data : */
+		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)
+		ISSMERROR("not implemented yet");
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel));
+
+	/*Create nodes and vertices: */
+	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->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+
+	/*Build Vertices dataset*/
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i,iomodel));
+
+		}
+	}
+
+	/*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
+	for (i=0;i<iomodel->numberofelements;i++){
+		for (j=0;j<3;j++){
+
+			if(iomodel->my_nodes[3*i+j]){ 
+
+				//Get index of the vertex on which the current node is located
+				vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
+				ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);
+
+				//Compute Node index (id-1)
+				node_index=3*i+j;
+
+				/*Add node to nodes dataset: */
+				nodes->AddObject(new Node(vertex_index,node_index,iomodel));
+
+			}
+		}
+	}
+
+	/*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->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	xfree((void**)&iomodel->gridonicesheet);
+	xfree((void**)&iomodel->gridoniceshelf);
+
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+}
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp	(revision 3588)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp	(revision 3588)
@@ -0,0 +1,118 @@
+/*! \file CreateLoadsBalancedthickness2.c:
+ */
+
+#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 "../../include/typedefs.h"
+#include "../IoModel.h"
+
+void	CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j;
+	int i1,i2;
+	int pos1,pos2;
+	double e1,e2;
+
+	/*Output*/
+	DataSet* loads=NULL;
+
+	/*numericalflux intermediary data: */
+	char numericalflux_type[NUMERICALFLUXSTRING];
+	int  numericalflux_id;
+	int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
+	int  numericalflux_elem_id;
+	double numericalflux_h[MAX_NUMERICALFLUX_NODES];
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Get edges and elements*/
+	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+
+	/*First load data:*/
+	for (i=0;i<iomodel->numberofedges;i++){
+
+		/*Get left and right elements*/
+		e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
+		e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
+
+		/*Now, if this element is not in the partition, pass: */
+		if(!iomodel->my_elements[(int)e1]) continue;
+
+		/*Create load*/
+		numericalflux_id=i+1; //Matlab indexing
+		numericalflux_elem_id=(int)e1+1;//id is in matlab index
+
+		/*1: Get vertices ids*/
+		i1=(int)iomodel->edges[4*i+0];
+		i2=(int)iomodel->edges[4*i+1];
+
+		if (!isnan(e2)){
+			strcpy(numericalflux_type,"internal");
+
+			/*Now, we must get the nodes of the 4 nodes located on the edge*/
+
+			/*2: Get the column where these ids are located in the index*/
+			pos1=pos2=UNDEF;
+			for(j=0;j<3;j++){
+				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
+				if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1;
+			}
+			ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
+
+			/*3: We have the id of the elements and the position of the vertices in the index
+			 * we can compute their dofs!*/
+			numericalflux_node_ids[0]=3*(int)e1+pos1;       //ex: 1 2 3
+			numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1
+			numericalflux_node_ids[2]=3*(int)e2+pos2;           //ex: 1 2 3
+			numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2
+
+			numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1];
+			numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
+			numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1];
+			numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1];
+		}
+		else{
+			strcpy(numericalflux_type,"boundary");
+
+			/*2: Get the column where these ids are located in the index*/
+			pos1==UNDEF;
+			for(j=0;j<3;j++){
+				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
+			}
+			ISSMASSERT(pos1!=UNDEF);
+
+			/*3: We have the id of the elements and the position of the vertices in the index
+			 * we can compute their dofs!*/
+			numericalflux_node_ids[0]=3*(int)e1+pos1;
+			numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1;
+
+			numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1];
+			numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
+			numericalflux_h[2]=UNDEF;
+			numericalflux_h[3]=UNDEF;
+		}
+
+		loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h));
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->edges);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->thickness);
+
+	/*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();
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp	(revision 3588)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp	(revision 3588)
@@ -0,0 +1,104 @@
+/*!\file: CreateParametersPrognostic.cpp
+ * \brief driver for creating parameters dataset, for prognostic analysis.
+ */ 
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+void CreateParametersBalancedthickness2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	Param*   param = NULL;
+	DataSet* parameters=NULL;
+	int      count;
+	int      i;
+	int      dim;
+	double*  vx_g=NULL;
+	double*  vy_g=NULL;
+
+	/*recover parameters : */
+	parameters=*pparameters;
+
+	count=parameters->Size();
+
+	/*Get vx and vy: */
+	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+
+	vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
+	vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
+
+	if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts;
+	if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts;
+
+	count++;
+	param= new Param(count,"vx_g",DOUBLEVEC);
+	param->SetDoubleVec(vx_g,iomodel->numberofvertices,1);
+	parameters->AddObject(param);
+	count++;
+	param= new Param(count,"vy_g",DOUBLEVEC);
+	param->SetDoubleVec(vy_g,iomodel->numberofvertices,1);
+	parameters->AddObject(param);
+
+	xfree((void**)&iomodel->vx);
+	xfree((void**)&iomodel->vy);
+	xfree((void**)&vx_g);
+	xfree((void**)&vy_g);
+
+	/*Get thickness: */
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+
+	count++;
+	param= new Param(count,"h_g",DOUBLEVEC);
+	if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1);
+	else param->SetDoubleVec(iomodel->thickness,0,0);
+	parameters->AddObject(param);
+
+	/*Free thickness: */
+	xfree((void**)&iomodel->thickness);
+
+	/*Get dhdt: */
+	IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
+	if(iomodel->dhdt) for(i=0;i<iomodel->numberofvertices;i++)iomodel->dhdt[i]=iomodel->dhdt[i]/iomodel->yts;
+
+	count++;
+	param= new Param(count,"dhdt_g",DOUBLEVEC);
+	if(iomodel->dhdt) param->SetDoubleVec(iomodel->dhdt,iomodel->numberofvertices,1);
+	else param->SetDoubleVec(iomodel->dhdt,0,1);
+	parameters->AddObject(param);
+
+	/*Free dhdt: */
+	xfree((void**)&iomodel->dhdt);
+
+	/*Get melting: */
+	IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting");
+	if(iomodel->melting) for(i=0;i<iomodel->numberofvertices;i++)iomodel->melting[i]=iomodel->melting[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"m_g",DOUBLEVEC);
+	if(iomodel->melting) param->SetDoubleVec(iomodel->melting,iomodel->numberofvertices,1);
+	else param->SetDoubleVec(iomodel->melting,0,1);
+	parameters->AddObject(param);
+
+	/*Free melting: */
+	xfree((void**)&iomodel->melting);
+
+	/*Get accumulation: */
+	IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation");
+	if(iomodel->accumulation) for(i=0;i<iomodel->numberofvertices;i++)iomodel->accumulation[i]=iomodel->accumulation[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"a_g",DOUBLEVEC);
+	if(iomodel->accumulation) param->SetDoubleVec(iomodel->accumulation,iomodel->numberofvertices,1);
+	else param->SetDoubleVec(iomodel->accumulation,0,0);
+	parameters->AddObject(param);
+
+	/*Free accumulation: */
+	xfree((void**)&iomodel->accumulation);
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 3588)
@@ -100,4 +100,12 @@
 
 	}
+	else if (iomodel->analysis_type==Balancedthickness2AnalysisEnum){
+
+		CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle);
+		CreateParametersBalancedthickness2(pparameters,iomodel,iomodel_handle);
+
+	}
 	else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum){
 
Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 3588)
@@ -247,5 +247,5 @@
 	if (
 				iomodel->analysis_type==Prognostic2AnalysisEnum ||
-				iomodel->analysis_type==BalancedthicknessAnalysisEnum
+				iomodel->analysis_type==Balancedthickness2AnalysisEnum
 				)
 	 param->SetDouble(3*iomodel->numberofelements);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 3588)
@@ -263,4 +263,10 @@
 	void  CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
+	/*balancedthickness2:*/
+	void	CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void	CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+	void  CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+	void  CreateParametersBalancedthickness2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+
 	/*balancedvelocities:*/
 	void	CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
Index: /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp	(revision 3587)
+++ /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp	(revision 3588)
@@ -22,5 +22,5 @@
 void  Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){
 	
-	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum)
+	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum)
 		DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle);
 	else
Index: /issm/trunk/src/c/objects/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3587)
+++ /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3588)
@@ -225,5 +225,5 @@
 		if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
 	}
-	else if (analysis_type==BalancedthicknessAnalysisEnum){
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
 		/*No transient term is involved*/
 		dt=1;
@@ -256,4 +256,7 @@
 		GetParameterValue(&vy, &vy_list[0],gauss_coord);
 		UdotN=vx*normal[0]+vy*normal[1];
+		if (fabs(UdotN)<1.0e-9) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN);
+		//if (fabs(UdotN)>0 && fabs(UdotN)< 1.0e-7) UdotN= 1.0e-7;
+		//if (fabs(UdotN)<0 && fabs(UdotN)>-1.0e-7) UdotN=-1.0e-7;
 
 		/*Get L and B: */
@@ -341,5 +344,5 @@
 		if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
 	}
-	else if (analysis_type==BalancedthicknessAnalysisEnum){
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
 		/*No transient term is involved*/
 		dt=1;
@@ -483,5 +486,5 @@
 		if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
 	}
-	else if (analysis_type==BalancedthicknessAnalysisEnum){
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
 		/*No transient term is involved*/
 		dt=1;
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3587)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3588)
@@ -77,5 +77,5 @@
 	/*hooks: */
 	//go recover node ids, needed to initialize the node hook.
-	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum){
+	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){
 		/*Discontinuous Galerkin*/
 		tria_node_ids[0]=3*i+1;
@@ -598,4 +598,8 @@
 		CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
+
+		CreateKMatrixBalancedthickness2( Kgg,inputs,analysis_type,sub_analysis_type);
+	}
 	else if (analysis_type==BalancedvelocitiesAnalysisEnum){
 
@@ -611,4 +615,217 @@
 /*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/
 void  Tria::CreateKMatrixBalancedthickness(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_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+
+	double Jdettria;
+
+	/*input parameters for structural analysis (diagnostic): */
+	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;
+
+	/*dynamic objects pointed to by hooks: */
+	Node**  nodes=NULL;
+	Matpar* matpar=NULL;
+	Matice* matice=NULL;
+	Numpar* numpar=NULL;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover objects from hooks: */
+	nodes=(Node**)hnodes.deliverp();
+	matpar=(Matpar*)hmatpar.delivers();
+	matice=(Matice*)hmatice.delivers();
+	numpar=(Numpar*)hnumpar.delivers();
+
+	/*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)ISSMERROR(" 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: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	//Create Artificial diffusivity once for all if requested
+	if(numpar->artdiff){
+		//Get the Jacobian determinant
+		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		//Build K matrix (artificial diffusivity matrix)
+		v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
+		v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
+
+		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
+		K[1][1]=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;
+
+		//Create DL and DLprime matrix
+		DL[0][0]=DL_scalar*dvxdx;
+		DL[1][1]=DL_scalar*dvydy;
+
+		DLprime[0][0]=DL_scalar*vx;
+		DLprime[1][1]=DL_scalar*vy;
+
+		//Do the triple product tL*D*L. 
+		//Ke_gg_thickness=B'*DLprime*Bprime;
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DL[0][0],2,2,0,
+					&B[0][0],2,numdof,0,
+					&Ke_gg_thickness1[0][0],0);
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DLprime[0][0],2,2,0,
+					&Bprime[0][0],2,numdof,0,
+					&Ke_gg_thickness2[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_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 Tria::CreateKMatrixBalancedthickness2 {{{1*/
+void  Tria::CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -2109,4 +2326,8 @@
 		CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
 	}
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
+
+		CreatePVectorBalancedthickness2( pg,inputs,analysis_type,sub_analysis_type);
+	}
 	else if (analysis_type==BalancedvelocitiesAnalysisEnum){
 
@@ -2121,4 +2342,107 @@
 /*FUNCTION Tria::CreatePVectorBalancedthickness {{{1*/
 void  Tria::CreatePVectorBalancedthickness(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;
+	double  thickness_list[numgrids]={0.0};
+	double  thickness_g;
+	int     dofs[1]={0};
+	int     found=0;
+
+	/*dynamic objects pointed to by hooks: */
+	Node**  nodes=NULL;
+	Matpar* matpar=NULL;
+	Matice* matice=NULL;
+	Numpar* numpar=NULL;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover objects from hooks: */
+	nodes=(Node**)hnodes.deliverp();
+	matpar=(Matpar*)hmatpar.delivers();
+	matice=(Matice*)hmatice.delivers();
+	numpar=(Numpar*)hnumpar.delivers();
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find accumulation in inputs!");
+	found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find melting in inputs!");
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&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 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 Tria::CreatePVectorBalancedthickness2 {{{1*/
+void  Tria::CreatePVectorBalancedthickness2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
 
 
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 3587)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 3588)
@@ -117,4 +117,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  CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorBalancedthickness2(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);
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 3587)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 3588)
@@ -135,4 +135,7 @@
 	if(analysis_type==BalancedthicknessAnalysisEnum){
 		fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
+	}
+	if(analysis_type==Balancedthickness2AnalysisEnum){
+		fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
 	}
 	if(analysis_type==BalancedvelocitiesAnalysisEnum){
Index: /issm/trunk/src/c/parallel/balancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 3587)
+++ /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 3588)
@@ -25,10 +25,8 @@
 	Model* model=NULL;
 
-	double* vx_g=NULL;
-	double* vy_g=NULL;
-	double* m_g=NULL;
-	double* a_g=NULL;
-	double* h_g=NULL;
-	double* dhdt_g=NULL;
+	Vec     u_g=NULL;
+	double* u_g_serial=NULL;
+	double* melting_g=NULL;
+	double* accumulation_g=NULL;
 	double  dt;
 	double  yts;
@@ -62,5 +60,5 @@
 	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
 
-	_printf_("recover input file name and output file name:\n");
+	_printf_("recover , input file name and output file name:\n");
 	inputfilename=argv[2];
 	outputfilename=argv[3];
@@ -86,19 +84,13 @@
 	_printf_("initialize inputs:\n");
 	
-	model->FindParam(&vx_g,NULL,NULL,"vx_g",BalancedthicknessAnalysisEnum);
-	model->FindParam(&vy_g,NULL,NULL,"vy_g",BalancedthicknessAnalysisEnum);
-	model->FindParam(&m_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum);
-	model->FindParam(&a_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum);
-	model->FindParam(&h_g,NULL,NULL,"h_g",BalancedthicknessAnalysisEnum);
-	model->FindParam(&dhdt_g,NULL,NULL,"dhdt_g",BalancedthicknessAnalysisEnum);
+	model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedthicknessAnalysisEnum);
+	model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum);
+	model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum);
 	model->FindParam(&numberofnodes,"numberofnodes");
 	
 	inputs=new ParameterInputs;
-	inputs->Add("vx",vx_g,1,numberofnodes);
-	inputs->Add("vy",vy_g,1,numberofnodes);
-	inputs->Add("melting",m_g,1,numberofnodes);
-	inputs->Add("accumulation",a_g,1,numberofnodes);
-	inputs->Add("thickness",h_g,1,numberofnodes);
-	inputs->Add("dhdt",dhdt_g,1,numberofnodes);
+	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");
@@ -146,13 +138,6 @@
 
 	/*Free ressources:*/
-	xfree((void**)&vx_g);
-	xfree((void**)&vy_g);
-	xfree((void**)&m_g);
-	xfree((void**)&a_g);
-	xfree((void**)&dhdt_g);
-	xfree((void**)&h_g);
 	delete processedresults;
 	delete results;
-	delete model;
 	delete inputs;
 
Index: /issm/trunk/src/c/parallel/balancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness2.cpp	(revision 3588)
+++ /issm/trunk/src/c/parallel/balancedthickness2.cpp	(revision 3588)
@@ -0,0 +1,171 @@
+/*!\file:  balancedthickness2.cpp
+ * \brief: balancedthickness2 solution
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#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* vx_g=NULL;
+	double* vy_g=NULL;
+	double* m_g=NULL;
+	double* a_g=NULL;
+	double* h_g=NULL;
+	double* dhdt_g=NULL;
+	double  dt;
+	double  yts;
+	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_))
+	ISSMERROR(" 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,Balancedthickness2AnalysisEnum);
+
+	/*recover parameters: */
+	model->FindParam(&waitonlock,"waitonlock");
+	model->FindParam(&qmu_analysis,"qmu_analysis");
+
+	_printf_("initialize inputs:\n");
+	
+	model->FindParam(&vx_g,NULL,NULL,"vx_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&vy_g,NULL,NULL,"vy_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&m_g,NULL,NULL,"m_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&a_g,NULL,NULL,"a_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&h_g,NULL,NULL,"h_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&dhdt_g,NULL,NULL,"dhdt_g",Balancedthickness2AnalysisEnum);
+	model->FindParam(&numberofnodes,"numberofnodes");
+	
+	inputs=new ParameterInputs;
+	inputs->Add("vx",vx_g,1,numberofnodes);
+	inputs->Add("vy",vy_g,1,numberofnodes);
+	inputs->Add("melting",m_g,1,numberofnodes);
+	inputs->Add("accumulation",a_g,1,numberofnodes);
+	inputs->Add("thickness",h_g,1,numberofnodes);
+	inputs->Add("dhdt",dhdt_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 balancedthickness2 analysis: */
+		_printf_("call computational core:\n");
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		balancedthickness2_core(results,model,inputs);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+	}
+	else{
+
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on balancedthickness2 core:\n");
+	
+		#ifdef _HAVE_DAKOTA_ 
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		Qmux(model,inputs,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	 	#else
+		ISSMERROR(" Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	/*Add analysis_type to results: */
+	result=new Result(results->Size()+1,0,1,"analysis_type","balancedthickness2");
+	results->AddObject(result);
+	
+	_printf_("process results:\n");
+	ProcessResults(&processedresults,results,model,Balancedthickness2AnalysisEnum);
+	
+	_printf_("write results to disk:\n");
+	OutputResults(processedresults,outputfilename);
+
+	if (waitonlock>0){
+		_printf_("write lock file:\n");
+		WriteLockFile(lockname);
+	}
+
+	/*Free ressources:*/
+	xfree((void**)&vx_g);
+	xfree((void**)&vy_g);
+	xfree((void**)&m_g);
+	xfree((void**)&a_g);
+	xfree((void**)&dhdt_g);
+	xfree((void**)&h_g);
+	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_("\n   %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+
+	/*end module: */
+	MODULEEND();
+
+	return 0; //unix success return;
+}
Index: /issm/trunk/src/c/parallel/balancedthickness2_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness2_core.cpp	(revision 3588)
+++ /issm/trunk/src/c/parallel/balancedthickness2_core.cpp	(revision 3588)
@@ -0,0 +1,73 @@
+/*!\file: balancedthickness2_core.cpp
+ * \brief: core of the balancedthickness2 solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./parallel.h"
+#include "../issm.h"
+
+void balancedthickness2_core(DataSet* results,Model* model,ParameterInputs* inputs){
+
+	extern int my_rank;
+
+	/*output: */
+	Result* result=NULL;
+
+	/*intermediary: */
+	Vec vx_g=NULL;
+	Vec vy_g=NULL;
+
+	/*solutions: */
+	Vec h_g=NULL;
+
+	/*flags: */
+	int verbose=0;
+	int numberofdofspernode;
+	int numberofnodes;
+	int numberofvertices;
+	int dofs[1]={1};
+
+	/*fem balancedthickness2 model: */
+	FemModel* fem_p=NULL;
+
+
+	fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
+
+	//first recover parameters common to all solutions
+	model->FindParam(&verbose,"verbose");
+	model->FindParam(&numberofnodes,"numberofnodes");
+	model->FindParam(&numberofvertices,"numberofvertices");
+	model->FindParam(&numberofdofspernode,"numberofdofspernode");
+
+	_printf_("depth averaging velocity...\n");
+	vx_g=inputs->Get("vx",&dofs[0],1);
+	vy_g=inputs->Get("vy",&dofs[0],1);
+	/* NOT WORKING YET....
+	FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
+	FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
+	*/
+	inputs->Add("vx_average",vx_g,1,numberofvertices);
+	inputs->Add("vy_average",vy_g,1,numberofvertices);
+	
+	_printf_("call computational core:\n");
+	diagnostic_core_linear(&h_g,fem_p,inputs,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
+
+	_printf_("Averaging over vertices:\n");
+	FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
+
+//	_printf_("extrude computed thickness on all layers:\n");
+//	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"h_g",h_g);
+	results->AddObject(result);
+
+	/*Free ressources:*/
+	VecFree(&vx_g);
+	VecFree(&vy_g);
+	VecFree(&h_g);
+}
Index: /issm/trunk/src/c/parallel/balancedthickness_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness_core.cpp	(revision 3587)
+++ /issm/trunk/src/c/parallel/balancedthickness_core.cpp	(revision 3588)
@@ -19,6 +19,5 @@
 
 	/*intermediary: */
-	Vec vx_g=NULL;
-	Vec vy_g=NULL;
+	Vec u_g=NULL;
 
 	/*solutions: */
@@ -29,6 +28,5 @@
 	int numberofdofspernode;
 	int numberofnodes;
-	int numberofvertices;
-	int dofs[1]={1};
+	int dofs[2]={1,1};
 
 	/*fem balancedthickness model: */
@@ -36,4 +34,5 @@
 
 
+	/*recover fem model: */
 	fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
 
@@ -41,25 +40,16 @@
 	model->FindParam(&verbose,"verbose");
 	model->FindParam(&numberofnodes,"numberofnodes");
-	model->FindParam(&numberofvertices,"numberofvertices");
 	model->FindParam(&numberofdofspernode,"numberofdofspernode");
 
 	_printf_("depth averaging velocity...\n");
-	vx_g=inputs->Get("vx",&dofs[0],1);
-	vy_g=inputs->Get("vy",&dofs[0],1);
-	/* NOT WORKING YET....
-	FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
-	FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
-	*/
-	inputs->Add("vx_average",vx_g,1,numberofvertices);
-	inputs->Add("vy_average",vy_g,1,numberofvertices);
+	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->vertices,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(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
 
-	_printf_("Averaging over vertices:\n");
-	FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
-
-//	_printf_("extrude computed thickness on all layers:\n");
-//	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+	_printf_("extrude computed thickness on all layers:\n");
+	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
 
 	/*Plug results into output dataset: */
@@ -68,6 +58,5 @@
 
 	/*Free ressources:*/
-	VecFree(&vx_g);
-	VecFree(&vy_g);
+	VecFree(&u_g);
 	VecFree(&h_g);
 }
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 3587)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 3588)
@@ -19,4 +19,5 @@
 void prognostic2_core(DataSet* results,Model* model, ParameterInputs* inputs);
 void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs);
+void balancedthickness2_core(DataSet* results,Model* model, ParameterInputs* inputs);
 void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs);
 void slopecompute_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 3587)
+++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 3588)
@@ -59,4 +59,7 @@
 		numdofs=1;
 	}
+	else if (analysis_type==Balancedthickness2AnalysisEnum){
+		numdofs=1;
+	}
 	else if (analysis_type==BalancedvelocitiesAnalysisEnum){
 		numdofs=1;