Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 815)
+++ /issm/trunk/src/c/Makefile.am	(revision 816)
@@ -222,4 +222,6 @@
 					./UpdateFromInputsx/UpdateFromInputsx.h\
 					./UpdateFromInputsx/UpdateFromInputsx.cpp\
+					./UpdateGeometryx/UpdateGeometryx.h\
+					./UpdateGeometryx/UpdateGeometryx.cpp\
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
@@ -492,4 +494,6 @@
 					./UpdateFromInputsx/UpdateFromInputsx.h\
 					./UpdateFromInputsx/UpdateFromInputsx.cpp\
+					./UpdateGeometryx/UpdateGeometryx.h\
+					./UpdateGeometryx/UpdateGeometryx.cpp\
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 816)
@@ -122,4 +122,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onsheet;
+	int node_onshelf;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -579,4 +581,6 @@
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 	/*Get number of dofs per node: */
@@ -612,4 +616,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -626,5 +633,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -665,5 +672,7 @@
 	xfree((void**)&model->gridonhutter);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 816)
@@ -94,4 +94,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -364,4 +366,6 @@
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 	
@@ -397,4 +401,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -411,5 +418,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -443,5 +450,7 @@
 	xfree((void**)&model->gridonhutter);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 816)
@@ -102,4 +102,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -461,4 +463,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -475,5 +480,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -516,4 +521,6 @@
 	xfree((void**)&model->gridonstokes);
 	xfree((void**)&model->borderstokes);
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
 
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 816)
@@ -101,4 +101,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -347,4 +349,6 @@
 	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 	
@@ -380,5 +384,8 @@
 		
 		node_onbed=(int)model->gridonbed[i];
-		node_onsurface=(int)model->gridonsurface[i];	
+		node_onsurface=(int)model->gridonsurface[i];
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+	
 		if (isnan(model->uppernodes[i])){
 			node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
@@ -389,5 +396,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*Add node to nodes dataset: */
@@ -413,5 +420,7 @@
 	xfree((void**)&model->gridonsurface);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp	(revision 816)
@@ -101,4 +101,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -392,4 +394,6 @@
 	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 
@@ -426,4 +430,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -440,5 +447,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -470,4 +477,6 @@
 	xfree((void**)&model->gridonsurface);
 	xfree((void**)&model->uppernodes);
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 816)
@@ -76,4 +76,5 @@
 	model->elementoniceshelf=NULL;
 	model->gridonicesheet=NULL;
+	model->gridoniceshelf=NULL;
 
 	model->drag_type=0;
@@ -236,4 +237,5 @@
 	xfree((void**)&model->elementoniceshelf);
 	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
 	xfree((void**)&model->segmentonneumann_diag);
 	xfree((void**)&model->segmentonneumann_diag_stokes);
Index: /issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 816)
@@ -73,4 +73,5 @@
 	double* elementoniceshelf;
 	double* gridonicesheet;
+	double* gridoniceshelf;
 
 	/*friction: */
Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 816)
@@ -123,4 +123,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -442,4 +444,6 @@
 	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 
@@ -476,4 +480,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -490,5 +497,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -523,5 +530,7 @@
 	xfree((void**)&model->gridonsurface);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 816)
@@ -101,4 +101,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -398,4 +400,6 @@
 	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 
@@ -432,4 +436,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -446,5 +453,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*set single point constraints.: */
@@ -479,5 +486,7 @@
 	xfree((void**)&model->gridonsurface);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 
 	/*Keep partitioning information into model*/
Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 815)
+++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 816)
@@ -102,4 +102,6 @@
 	int node_onbed;
 	int node_onsurface;
+	int node_onshelf;
+	int node_onsheet;
 	int node_upper_node_id;
 	int node_numdofs;
@@ -387,4 +389,6 @@
 	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
 	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
 
 	/*Get number of dofs per node: */
@@ -416,4 +420,7 @@
 		node_onbed=(int)model->gridonbed[i];
 		node_onsurface=(int)model->gridonsurface[i];	
+		node_onshelf=(int)model->gridoniceshelf[i];	
+		node_onsheet=(int)model->gridonicesheet[i];	
+
 		if (strcmp(model->meshtype,"3d")==0){
 			if (isnan(model->uppernodes[i])){
@@ -430,5 +437,5 @@
 
 		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
 
 		/*Add node to nodes dataset: */
@@ -454,5 +461,7 @@
 	xfree((void**)&model->gridonsurface);
 	xfree((void**)&model->uppernodes);
-		
+	xfree((void**)&model->gridonicesheet);
+	xfree((void**)&model->gridoniceshelf);
+	
 	cleanup_and_return:
 
Index: /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.cpp
===================================================================
--- /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.cpp	(revision 816)
+++ /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.cpp	(revision 816)
@@ -0,0 +1,126 @@
+/*!\file UpdateGeometryx
+ * \brief: update geometry node by node
+ */
+
+#include "./UpdateGeometryx.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "UpdateGeometryx"
+
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../DataSet/DataSet.h"
+
+void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, 
+		DataSet* elements, DataSet* nodes,DataSet* loads, DataSet* materials, DataSet* parameters,
+		Vec vec_newthickness,Vec vec_bed,Vec vec_surface){
+
+	int i;
+	int dof;
+
+	/*serialized input: */
+	double* newthickness=NULL;
+	double* bed=NULL;
+	double* surface=NULL;
+
+	/*objects: */
+	Node* node=NULL;
+	Object* object=NULL;
+	Matpar* matpar=NULL;
+
+	/*parameters: */
+	double h,b,s;
+	double rho_ice,rho_water;
+
+	/*output: */
+	Vec outthickness=NULL;
+	Vec outbed=NULL;
+	Vec outsurface=NULL;
+
+	/*Duplicate inputs to outputs, do not copy values!: */
+	VecDuplicate(vec_newthickness,&outthickness);
+	VecDuplicate(vec_bed,&outbed);
+	VecDuplicate(vec_surface,&outsurface);
+
+	/*First, get elements and loads configured: */
+	elements->Configure(elements,loads, nodes, materials);
+	materials->Configure(elements, loads, nodes, materials);
+	loads->Configure(elements, loads, nodes, materials);
+
+	/*Serialize inputs :*/
+	VecToMPISerial(&newthickness,vec_newthickness);
+	VecToMPISerial(&bed,vec_bed);
+	VecToMPISerial(&surface,vec_surface);
+
+	/*First, find the matpar object in materials: */
+	for(i=0;i<materials->Size();i++){
+		Object* object=materials->GetObjectByOffset(i);
+		if  (object->Enum()==MatparEnum()){
+			matpar=(Matpar*)object;
+			break;
+		}
+	}
+	/*recover material parameters: */
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+
+	/*Go through nodes and for each node, update the thickness, bed and surface, using the 
+	 * new thickness computed in the prognostic_core part of the transient solution sequence: */
+
+	for(i=0;i<nodes->Size();i++){
+
+		/*get i'th node: */
+		node=(Node*)nodes->GetObjectByOffset(i);
+
+		/*first, recover thickness, surface and bed for this node: */
+		dof=node->GetDofList1();
+		h=newthickness[dof];
+		s=surface[dof];
+		b=bed[dof];
+
+		//First, check that thickness is >0
+		if(h<0)h=10; 
+
+		//For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. 
+		if (node->IsOnSheet()){
+			s=b+h;
+		}
+
+		//For grids on ice shelt, we have hydrostatic equilibrium (for now)
+		if (node->IsOnShelf()){
+
+			b=-rho_ice/rho_water*h;
+			s=(1-rho_ice/rho_water)*h;
+
+		}
+
+		/*Ok, plug values of thickness, surafce and bed into our outputs: */
+		VecSetValues(outthickness,1,&dof,&h,INSERT_VALUES);
+		VecSetValues(outbed,1,&dof,&b,INSERT_VALUES);
+		VecSetValues(outsurface,1,&dof,&s,INSERT_VALUES);
+
+	}
+
+	/*Assemble vectors: */
+	VecAssemblyBegin(outthickness);
+	VecAssemblyEnd(outthickness);
+
+	VecAssemblyBegin(outbed);
+	VecAssemblyEnd(outbed);
+
+	VecAssemblyBegin(outsurface);
+	VecAssemblyEnd(outsurface);
+
+	/*Free ressources:*/
+	xfree((void**)&newthickness);
+	xfree((void**)&bed);
+	xfree((void**)&surface);
+
+	
+	/*Assign output pointers: */
+	*poutthickness=outthickness;
+	*poutbed=outbed;
+	*poutsurface=outsurface;
+}
Index: /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.h
===================================================================
--- /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.h	(revision 816)
+++ /issm/trunk/src/c/UpdateGeometryx/UpdateGeometryx.h	(revision 816)
@@ -0,0 +1,17 @@
+/*!\file:  UpdateGeometryx.h
+ * \brief header file for updating geometry
+ */ 
+
+#ifndef _UPDATEGEOMETRYX_H
+#define _UPDATEGEOMETRYX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, 
+		DataSet* elements, DataSet* nodes,DataSet* loads, DataSet* materials, DataSet* parameters,
+		Vec newthickness,Vec bed,Vec surface);
+
+#endif  /* _UPDATEGEOMETRYX_H */
+
Index: /issm/trunk/src/c/issm.h
===================================================================
--- /issm/trunk/src/c/issm.h	(revision 815)
+++ /issm/trunk/src/c/issm.h	(revision 816)
@@ -33,4 +33,5 @@
 #include "./SystemMatricesx/SystemMatricesx.h"
 #include "./UpdateFromInputsx/UpdateFromInputsx.h"
+#include "./UpdateGeometryx/UpdateGeometryx.h"
 #include "./PenaltySystemMatricesx/PenaltySystemMatricesx.h"
 #include "./Reducematrixfromgtofx/Reducematrixfromgtofx.h"
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 815)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 816)
@@ -21,5 +21,5 @@
 	return;
 }
-Node::Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],int node_onbed,int node_onsurface,int node_upper_node_id){
+Node::Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],int node_onbed,int node_onsurface,int node_upper_node_id,int node_onshelf,int node_onsheet){
 
 	int i;
@@ -33,4 +33,6 @@
 	onbed=node_onbed;
 	onsurface=node_onsurface;
+	onshelf=node_onshelf;
+	onsheet=node_onsheet;
 
 	/*Initialize sets: */
@@ -66,4 +68,6 @@
 	printf("   onbed: %i\n",onbed);
 	printf("   onsurface: %i\n",onsurface);
+	printf("   onshelf: %i\n",onshelf);
+	printf("   onsheet: %i\n",onsheet);
 	printf("   upper_node_id=%i\n",upper_node_id);
 	printf("   upper_node_offset=%i\n",upper_node_offset);
@@ -100,4 +104,6 @@
 	printf("   onbed: %i\n",onbed);
 	printf("   onsurface: %i\n",onsurface);
+	printf("   onshelf: %i\n",onshelf);
+	printf("   onsheet: %i\n",onsheet);
 	printf("   upper_node_id=%i\n",upper_node_id);
 	printf("   upper_node_offset=%i\n",upper_node_offset);
@@ -143,4 +149,6 @@
 	memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
 	memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
+	memcpy(marshalled_dataset,&onshelf,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf);
+	memcpy(marshalled_dataset,&onsheet,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet);
 	memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
 	memcpy(marshalled_dataset,&doflist1,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
@@ -166,4 +174,6 @@
 		sizeof(onbed)+
 		sizeof(onsurface)+
+		sizeof(onshelf)+
+		sizeof(onsheet)+
 		sizeof(doflist)+
 		sizeof(doflist1)+
@@ -200,4 +210,6 @@
 	memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
 	memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
+	memcpy(&onshelf,marshalled_dataset,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf);
+	memcpy(&onsheet,marshalled_dataset,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet);
 	memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
 	memcpy(&doflist1,marshalled_dataset,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
@@ -517,5 +529,12 @@
 	return onsurface;
 }
-		
+
+int   Node::IsOnShelf(){
+	return onshelf;
+}
+
+int   Node::IsOnSheet(){
+	return onsheet;
+}		
 void  Node::FreezeDof(int dof){
 	
Index: /issm/trunk/src/c/objects/Node.h
===================================================================
--- /issm/trunk/src/c/objects/Node.h	(revision 815)
+++ /issm/trunk/src/c/objects/Node.h	(revision 816)
@@ -23,4 +23,6 @@
 		int	    onbed; /*! for 3d, on bedrock*/
 		int	    onsurface; /*! for 3d, on surface*/
+		int	    onshelf; 
+		int	    onsheet;
 
 		/*for dof constraining: */
@@ -41,5 +43,5 @@
 
 		Node();
-		Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],int node_onbed,int node_onsurface,int upper_node_id);
+		Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],int node_onbed,int node_onsurface,int upper_node_id,int onshelf,int onsheet);
 		~Node();
 
@@ -79,5 +81,6 @@
 		void  FreezeDof(int dof);
 		void  VelocityDepthAverageAtBase(Vec ug,double* ug_serial);
-
+		int   IsOnShelf();
+		int   IsOnSheet();
 };
 
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 815)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 816)
@@ -74,4 +74,5 @@
 WriteData(fid,md.elementoniceshelf,'Mat','elementoniceshelf');
 WriteData(fid,md.gridonicesheet,'Mat','gridonicesheet');
+WriteData(fid,md.gridoniceshelf,'Mat','gridoniceshelf');
 
 WriteData(fid,md.segmentonneumann_diag,'Mat','segmentonneumann_diag');
Index: sm/trunk/src/m/solutions/cielo/UpdateGeometry.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/UpdateGeometry.m	(revision 815)
+++ 	(revision )
@@ -1,27 +1,0 @@
-function [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)
-%UPDATEGEOMETRY - update the geometry during after a time step
-%
-%   Update geometry once given a new thickness, in the transient solutions
-%
-%   Usage:
-%      [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)
-
-%initialize bed
-new_bed=zeros(md.numberofgrids,1);
-new_surface=zeros(md.numberofgrids,1);
-
-%First, check that thickness is >0
-pos=find(new_thickness<0);
-new_thickness(pos)=10; %minimum thickness
-
-%For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. 
-pos=find(md.gridonicesheet);
-new_bed(pos)=bed(pos);
-new_surface(pos)=bed(pos)+new_thickness(pos);
-
-%For grids on ice shelt, we have hydrostatic equilibrium (for now)
-pos=find(md.gridoniceshelf);
-new_bed(pos)=-md.rho_ice/md.rho_water*new_thickness(pos);
-new_surface(pos)=(1-md.rho_ice/md.rho_water)*new_thickness(pos);
-
-%Deal with grounding line migration: later on.
Index: /issm/trunk/src/m/solutions/cielo/transient2d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 815)
+++ /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 816)
@@ -78,5 +78,5 @@
 	%update surface and bed using the new thickness
 	disp(sprintf('%s','   updating geometry...'));
-	[new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,solution(n).b_g,solution(n).s_g);
+	[new_thickness,new_bed,new_surface]=UpdateGeometry(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,m_p.parameters,new_thickness,solution(n).b_g,solution(n).s_g);
 
 	%Record bed surface and thickness in the solution
@@ -95,5 +95,4 @@
 	time=time+dt;
 	n=n+1;
-
 end
 
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 815)
+++ /issm/trunk/src/mex/Makefile.am	(revision 816)
@@ -43,4 +43,5 @@
 				TriMeshRefine\
 				UpdateFromInputs\
+				UpdateGeometry\
 				VelocityExtrude\
 				VelocityDepthAverage
@@ -184,4 +185,7 @@
 			  UpdateFromInputs/UpdateFromInputs.h
 
+UpdateGeometry_SOURCES = UpdateGeometry/UpdateGeometry.cpp\
+			  UpdateGeometry/UpdateGeometry.h
+
 VelocityExtrude_SOURCES = VelocityExtrude/VelocityExtrude.cpp\
 			  VelocityExtrude/VelocityExtrude.h
Index: /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.cpp
===================================================================
--- /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.cpp	(revision 816)
+++ /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.cpp	(revision 816)
@@ -0,0 +1,73 @@
+/*\file UpdateGeometry.c
+ *\brief: update geometry from new thickness computed by prognostic_core.
+ */
+
+#include "./UpdateGeometry.h"
+
+void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+	/*diverse: */
+	int   noerr=1;
+
+	/*input datasets: */
+	DataSet* elements=NULL;
+	DataSet* nodes=NULL;
+	DataSet* loads=NULL;
+	DataSet* materials=NULL;
+	DataSet* parameters=NULL;
+	Vec      newthickness=NULL;
+	Vec      bed=NULL;
+	Vec      surface=NULL;
+
+	/* output datasets: */
+	Vec outbed=NULL;
+	Vec outsurface=NULL;
+	Vec outthickness=NULL;
+
+	/*Boot module: */
+	MODULEBOOT();
+
+	/*checks on arguments on the matlab side: */
+	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&UpdateGeometryUsage);
+        
+	/*Input datasets: */
+	FetchData((void**)&elements,NULL,NULL,ELEMENTS,"DataSet",NULL);
+	FetchData((void**)&nodes,NULL,NULL,NODES,"DataSet",NULL);
+	FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
+	FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
+	FetchData((void**)&parameters,NULL,NULL,PARAMETERS,"DataSet",NULL);
+	
+	FetchData((void**)&newthickness,NULL,NULL,NEWTHICKNESS,"Vector","Vec");
+	FetchData((void**)&bed,NULL,NULL,BED,"Vector","Vec");
+	FetchData((void**)&surface,NULL,NULL,SURFACE,"Vector","Vec");
+
+	/*!Generate internal degree of freedom numbers: */
+	UpdateGeometryx(&outthickness,&outbed,&outsurface, elements, nodes,loads, materials, parameters,newthickness,bed,surface);
+
+	/*write output data: */
+	WriteData(OUTTHICKNESS,outthickness,0,0,"Vector",NULL);
+	WriteData(OUTBED,outbed,0,0,"Vector",NULL);
+	WriteData(OUTSURFACE,outsurface,0,0,"Vector",NULL);
+
+	/*Free ressources: */
+	delete elements;
+	delete nodes;
+	delete loads;
+	delete materials;
+	delete parameters;
+	VecFree(&newthickness);
+	VecFree(&bed);
+	VecFree(&surface);
+	VecFree(&outbed);
+	VecFree(&outthickness);
+	VecFree(&outsurface);
+
+	/*end module: */
+	MODULEEND();
+}
+
+void UpdateGeometryUsage(void) {
+	_printf_("\n");
+	_printf_("   usage: [outthickness,outbed,outsrface] = %s(elements, nodes, materials, params,newthickness,bed,surface);\n",__FUNCT__);
+	_printf_("\n");
+}
Index: /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.h
===================================================================
--- /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.h	(revision 816)
+++ /issm/trunk/src/mex/UpdateGeometry/UpdateGeometry.h	(revision 816)
@@ -0,0 +1,41 @@
+
+/*
+	UpdateGeometry.h
+*/
+
+
+#ifndef _UPDATEGEOMETRY_H
+#define _UPDATEGEOMETRY_H
+
+/* local prototypes: */
+void UpdateGeometryUsage(void);
+
+#include "../../c/issm.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "UpdateGeometry"
+
+/* serial input macros: */
+#define ELEMENTS (mxArray*)prhs[0]
+#define NODES (mxArray*)prhs[1]
+#define LOADS (mxArray*)prhs[2]
+#define MATERIALS (mxArray*)prhs[3]
+#define PARAMETERS (mxArray*)prhs[4]
+#define NEWTHICKNESS (mxArray*)prhs[5]
+#define BED (mxArray*)prhs[6]
+#define SURFACE (mxArray*)prhs[7]
+
+/* serial output macros: */
+#define OUTTHICKNESS (mxArray**)&plhs[0]
+#define OUTBED (mxArray**)&plhs[1]
+#define OUTSURFACE (mxArray**)&plhs[2]
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  3
+#undef NRHS
+#define NRHS  8
+
+
+#endif  /* _UPDATEGEOMETRY_H */
+