Index: /issm/trunk/src/c/Container/Elements.cpp
===================================================================
--- /issm/trunk/src/c/Container/Elements.cpp	(revision 7088)
+++ /issm/trunk/src/c/Container/Elements.cpp	(revision 7089)
@@ -204,2 +204,18 @@
 }
 /*}}}*/
+/*FUNCTION Elements::NumberOfElements{{{1*/
+int Elements::NumberOfElements(void){
+
+	int local_nelem=0;
+	int numberofelements;
+
+	#ifdef _PARALLEL_
+	local_nelem=this->Size();
+	MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+	#else
+	numberofelements=this->Size();
+	#endif
+
+	return numberofelements;
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Elements.h
===================================================================
--- /issm/trunk/src/c/Container/Elements.h	(revision 7088)
+++ /issm/trunk/src/c/Container/Elements.h	(revision 7089)
@@ -33,4 +33,5 @@
 		void ToResults(Results* results,Parameters* parameters,int step, double time);
 		Patch* ResultsToPatch(void);
+		int   NumberOfElements(void);
 		/*}}}*/
 
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 7088)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 7089)
@@ -36,4 +36,5 @@
 	/*Fetch data needed: */
 	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	IoModelFetchData(&iomodel->elementconnectivity,NULL,NULL,iomodel_handle,"elementconnectivity");
 	IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
 	IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
@@ -48,6 +49,6 @@
 
 			/*Create and add tria element to elements dataset: */
-			if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,iomodel,nummodels));
-			else elements->AddObject(new Penta(i+1,i,iomodel,nummodels));
+			if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,i,iomodel,nummodels));
+			else elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels));
 
 			/*Create and add material property to materials dataset: */
@@ -58,4 +59,5 @@
 	/*Free data: */
 	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->elementconnectivity);
 	xfree((void**)&iomodel->upperelements);
 	xfree((void**)&iomodel->lowerelements);
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 7088)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 7089)
@@ -75,4 +75,5 @@
 	parameters->AddObject(new BoolParam(KffEnum,iomodel->kff));
 	parameters->AddObject(new BoolParam(IoGatherEnum,iomodel->io_gather));
+	parameters->AddObject(new BoolParam(GroundingLineMigrationEnum,iomodel->grounding_line_migration));
 
 	/*Deal with more complex parameters*/
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 7088)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 7089)
@@ -39,4 +39,5 @@
 	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
 	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+	if(iomodel->grounding_line_migration) IoModelFetchData(&iomodel->bathymetry,NULL,NULL,iomodel_handle,"bathymetry");
 	if (iomodel->dim==3){
 		IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 7088)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 7089)
@@ -32,5 +32,7 @@
 		virtual void   GetSolutionFromInputs(Vec solution)=0;
 		virtual int    GetNodeIndex(Node* node)=0;
+		virtual int    Sid()=0;
 		virtual bool   IsOnShelf()=0; 
+		virtual bool   IsNodeOnShelf()=0; 
 		virtual bool   IsOnBed()=0;
 		virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
@@ -79,5 +81,8 @@
 		virtual bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0;
 		virtual void   AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part)=0;
+		virtual int*   GetHorizontalNeighboorSids(void)=0;
 		virtual double TimeAdapt()=0;
+		virtual void   MigrateGroundingLine()=0;
+		virtual void   UpdateShelfStatus()=0;
 
 		/*Implementation: */
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 7089)
@@ -26,11 +26,15 @@
 /*FUNCTION Penta::Penta(){{{1*/
 Penta::Penta(){
+
+	int i;
+
 	this->nodes=NULL;
 	this->matice=NULL;
 	this->matpar=NULL;
-	this->neighbors=NULL;
+	this->verticalneighbors=NULL;
 	this->inputs=NULL;
 	this->parameters=NULL;
 	this->results=NULL;
+	for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
 }
 /*}}}*/
@@ -43,5 +47,5 @@
 /*}}}*/
 /*FUNCTION Penta::Penta(int id, int index, IoModel* iomodel,int nummodels) {{{1*/
-Penta::Penta(int penta_id, int index, IoModel* iomodel,int nummodels)
+Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)
 	:PentaRef(nummodels)
 	,PentaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id
@@ -59,4 +63,5 @@
 	/*id: */
 	this->id=penta_id;
+	this->sid=penta_sid;
 
 	/*Build neighbors list*/
@@ -66,4 +71,7 @@
 	else                                    penta_elements_ids[0]=(int)(iomodel->lowerelements[index]);
 	this->InitHookNeighbors(penta_elements_ids);
+
+	/*Build horizontalneighborsids list: */
+	for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
 
 	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
@@ -78,5 +86,5 @@
 	this->matice=NULL;
 	this->matpar=NULL;
-	this->neighbors=NULL;
+	this->verticalneighbors=NULL;
 }
 /*}}}*/
@@ -104,4 +112,5 @@
 	/*deal with Penta  copy fields: */
 	penta->id=this->id;
+	penta->sid=this->sid;
 	if(this->inputs){
 		penta->inputs=(Inputs*)this->inputs->Copy();
@@ -124,5 +133,8 @@
 	penta->matice=(Matice*)penta->hmatice->delivers();
 	penta->matpar=(Matpar*)penta->hmatpar->delivers();
-	penta->neighbors=(Penta**)penta->hneighbors->deliverp();
+	penta->verticalneighbors=(Penta**)penta->hneighbors->deliverp();
+
+	/*neighbors: */
+	for(i=0;i<3;i++)penta->horizontalneighborsids[i]=this->horizontalneighborsids[i];
 
 	return penta;
@@ -154,4 +166,5 @@
 	/*marshall Penta data: */
 	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
 	memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
 
@@ -195,4 +208,7 @@
 	xfree((void**)&marshalled_results);
 
+	/*marshall horizontal neighbors: */
+	memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
+
 	*pmarshalled_dataset=marshalled_dataset;
 	return;
@@ -211,4 +227,5 @@
 
 	return sizeof(id)
+		+sizeof(sid)
 		+hnodes_size
 		+sizeof(numanalyses)
@@ -219,4 +236,5 @@
 		+inputs->MarshallSize()
 		+results->MarshallSize()
+		+3*sizeof(int)
 		+sizeof(int); //sizeof(int) for enum type
 }
@@ -235,4 +253,5 @@
 	 *object data (thanks to DataSet::Demarshall):*/
 	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
 	memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
 
@@ -260,5 +279,5 @@
 	matice=NULL;
 	matpar=NULL;
-	neighbors=NULL;
+	verticalneighbors=NULL;
 
 	/*demarshall inputs and results: */
@@ -268,4 +287,7 @@
 	/*parameters: may not exist even yet, so let Configure handle it: */
 	this->parameters=NULL;
+
+	/*neighbors: */
+	memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
 
 	/*return: */
@@ -425,5 +447,5 @@
 	this->matice=(Matice*)this->hmatice->delivers();
 	this->matpar=(Matpar*)this->hmatpar->delivers();
-	this->neighbors=(Penta**)this->hneighbors->deliverp();
+	this->verticalneighbors=(Penta**)this->hneighbors->deliverp();
 
 	/*point parameters to real dataset: */
@@ -499,5 +521,5 @@
 
 	/*Checks in debugging {{{2*/
-	_assert_(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
+	_assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
 	/*}}}*/
 
@@ -1794,5 +1816,5 @@
 
 	/*if debugging mode, check that all pointers exist {{{2*/
-	_assert_(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
+	_assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
 	/*}}}*/
 
@@ -2878,5 +2900,5 @@
 	matice->DeepEcho();
 	matpar->DeepEcho();
-	printf("   neighbor ids: %i-%i\n",neighbors[0]->Id(),neighbors[1]->Id());
+	printf("   neighbor ids: %i-%i\n",verticalneighbors[0]->Id(),verticalneighbors[1]->Id());
 	printf("   parameters\n");
 	parameters->DeepEcho();
@@ -2885,4 +2907,7 @@
 	printf("   results\n");
 	results->DeepEcho();
+	printf("neighboor sids: \n");
+	printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
+
 	return;
 }
@@ -3063,4 +3088,12 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/
+int* Penta::GetHorizontalNeighboorSids(){
+
+	/*return PentaRef field*/
+	return &this->horizontalneighborsids[0];
+
+}
+/*}}}*/
 /*FUNCTION Penta::GetLowerElement{{{1*/
 Penta* Penta::GetLowerElement(void){
@@ -3068,5 +3101,5 @@
 	Penta* upper_penta=NULL;
 
-	upper_penta=(Penta*)neighbors[0]; //first one (0) under, second one (1) above
+	upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
 
 	return upper_penta;
@@ -3497,5 +3530,5 @@
 	Penta* upper_penta=NULL;
 
-	upper_penta=(Penta*)neighbors[1]; //first one under, second one above
+	upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
 
 	return upper_penta;
@@ -3533,4 +3566,11 @@
 
 	return z;
+}
+/*}}}*/
+/*FUNCTION Penta::Sid {{{1*/
+int    Penta::Sid(){
+	
+	return sid;
+
 }
 /*}}}*/
@@ -5191,4 +5231,19 @@
 }
 /*}}}*/
+/*FUNCTION Penta::IsNodeOnShelf {{{1*/
+bool   Penta::IsNodeOnShelf(){
+
+	int  i;
+	bool shelf=false;
+
+	for(i=0;i<6;i++){
+		if (nodes[i]->IsOnShelf()){
+			shelf=true;
+			break;
+		}
+	}
+	return shelf;
+}
+/*}}}*/
 /*FUNCTION Penta::IsOnSurface{{{1*/
 bool Penta::IsOnSurface(void){
@@ -5302,4 +5357,9 @@
 	/*Assign output pointers:*/
 	*pmaxvz=maxvz;
+}
+/*}}}*/
+/*FUNCTION Penta::MigrateGroundingLine{{{1*/
+void  Penta::MigrateGroundingLine(void){
+	_error_("not supported yet!");
 }
 /*}}}*/
@@ -5417,4 +5477,5 @@
 	int     numrows  = 0;
 	int     numnodes = 0;
+	int     temp_numnodes=0;
 
 	/*Go through all the results objects, and update the counters: */
@@ -5424,5 +5485,6 @@
 		numrows++;
 		/*now, how many vertices and how many nodal values for this result? :*/
-		numnodes=elementresult->NumberOfNodalValues(); //ask result object.
+		temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
+		if(temp_numnodes>numnodes)numnodes=temp_numnodes;
 	}
 
@@ -6071,2 +6133,7 @@
 }
 /*}}}*/
+/*FUNCTION Penta::UpdateShelfStatus{{{1*/
+void Penta::UpdateShelfStatus(void){
+	_error_("Not implemented yet");
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 7088)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 7089)
@@ -32,9 +32,11 @@
 
 		int          id;
+		int          sid;
 
 		Node       **nodes;        // 6 nodes
 		Matice      *matice;       // 1 material ice
 		Matpar      *matpar;       // 1 material parameter
-		Penta      **neighbors;   // 2 neighbors: first one under, second one above
+		Penta      **verticalneighbors;   // 2 neighbors: first one under, second one above
+		int          horizontalneighborsids[3];
 
 		Parameters  *parameters;   //pointer to solution parameters
@@ -44,5 +46,5 @@
 		/*Penta constructors and destructor: {{{1*/
 		Penta();
-		Penta(int penta_id,int i, IoModel* iomodel,int nummodels);
+		Penta(int penta_id,int penta_sid,int i, IoModel* iomodel,int nummodels);
 		~Penta();
 		/*}}}*/
@@ -88,4 +90,5 @@
 		void   GradjB(Vec gradient);
 		void   GradjDrag(Vec gradient);
+		int    Sid();
 		void   InputControlUpdate(double scalar,bool save_parameter);
 		void   InputArtificialNoise(int enum_type,double min, double max);
@@ -106,4 +109,5 @@
 		void   MaxVy(double* pmaxvy, bool process_units);
 		void   MaxVz(double* pmaxvz, bool process_units);
+		void   MigrateGroundingLine();
 		void   MinVel(double* pminvel, bool process_units);
 		void   MinVx(double* pminvx, bool process_units);
@@ -121,5 +125,7 @@
 		double SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
+		void   UpdateShelfStatus(void);
 		double TimeAdapt();
+		int*   GetHorizontalNeighboorSids(void);
 		/*}}}*/
 		/*Penta specific routines:{{{1*/
@@ -229,4 +235,5 @@
 		bool	  IsOnBed(void);
 		bool    IsOnShelf(void); 
+		bool    IsNodeOnShelf(void); 
 		bool    IsOnWater(void); 
 		double  MinEdgeLength(double xyz_list[6][3]);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 7089)
@@ -26,7 +26,10 @@
 Tria::Tria(){
 
+	int i;
+
 	this->nodes=NULL;
 	this->matice=NULL;
 	this->matpar=NULL;
+	for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
 	this->inputs=NULL;
 	this->parameters=NULL;
@@ -35,22 +38,28 @@
 }
 /*}}}*/
-/*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/
-Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels)
+/*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{1*/
+Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)
 	:TriaRef(nummodels)
 	,TriaHook(nummodels,index+1,iomodel->numberofelements+1){
-	/*id: */
-	this->id=tria_id;
-
-	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
-	this->parameters=NULL;
-
-	/*intialize inputs and results: */
-	this->inputs=new Inputs();
-	this->results=new Results();
-
-	/*initialize pointers:*/
-	this->nodes=NULL;
-	this->matice=NULL;
-	this->matpar=NULL;
+		
+		int i;
+		/*id: */
+		this->id=tria_id;
+		this->sid=tria_sid;
+
+		//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
+		this->parameters=NULL;
+
+		/*Build horizontalneighborsids list: */
+		for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
+
+		/*intialize inputs and results: */
+		this->inputs=new Inputs();
+		this->results=new Results();
+
+		/*initialize pointers:*/
+		this->nodes=NULL;
+		this->matice=NULL;
+		this->matpar=NULL;
 
 }
@@ -84,4 +93,5 @@
 	/*deal with Tria fields: */
 	tria->id=this->id;
+	tria->sid=this->sid;
 	if(this->inputs){
 		tria->inputs=(Inputs*)this->inputs->Copy();
@@ -104,4 +114,7 @@
 	tria->matice=(Matice*)tria->hmatice->delivers();
 	tria->matpar=(Matpar*)tria->hmatpar->delivers();
+
+	/*neighbors: */
+	for(i=0;i<3;i++)tria->horizontalneighborsids[i]=this->horizontalneighborsids[i];
 
 	return tria;
@@ -133,4 +146,5 @@
 	/*marshall Tria data: */
 	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
 	memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
 
@@ -174,4 +188,7 @@
 	xfree((void**)&marshalled_results);
 
+	/*marshall horizontal neighbors: */
+	memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
+
 	*pmarshalled_dataset=marshalled_dataset;
 	return;
@@ -190,4 +207,5 @@
 
 	return sizeof(id)
+	  +sizeof(sid)
 	  +hnodes_size
 	  +sizeof(numanalyses)
@@ -197,4 +215,5 @@
 	  +inputs->MarshallSize()
 	  +results->MarshallSize()
+	  +3*sizeof(int)
 	  +sizeof(int); //sizeof(int) for enum type
 }
@@ -213,4 +232,5 @@
 	 *object data (thanks to DataSet::Demarshall):*/
 	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
 	memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
 
@@ -244,4 +264,7 @@
 	/*parameters: may not exist even yet, so let Configure handle it: */
 	this->parameters=NULL;
+
+	/*neighbors: */
+	memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int);
 
 	/*return: */
@@ -2614,4 +2637,7 @@
 	if (results) results->DeepEcho();
 	else printf("results=NULL\n");
+
+	printf("neighboor sids: \n");
+	printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
 	
 	return;
@@ -2654,4 +2680,7 @@
 	if (results) results->Echo();
 	else printf("results=NULL\n");
+
+	printf("neighboor sids: \n");
+	printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
 }
 /*}}}*/
@@ -2715,4 +2744,12 @@
 	/*return TriaRef field*/
 	return this->element_type;
+
+}
+/*}}}*/
+/*FUNCTION Tria::GetHorizontalNeighboorSids {{{1*/
+int* Tria::GetHorizontalNeighboorSids(){
+
+	/*return TriaRef field*/
+	return &this->horizontalneighborsids[0];
 
 }
@@ -3340,4 +3377,11 @@
 }
 /*}}}*/
+/*FUNCTION Tria::Sid {{{1*/
+int    Tria::Sid(){
+	
+	return sid;
+
+}
+/*}}}*/
 /*FUNCTION Tria::InputArtificialNoise{{{1*/
 void  Tria::InputArtificialNoise(int enum_type,double min,double max){
@@ -3548,4 +3592,8 @@
 		for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1];
 		this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs));
+	}
+	if(iomodel->grounding_line_migration){
+		for(i=0;i<3;i++)nodeinputs[i]=iomodel->bathymetry[tria_vertex_ids[i]-1];
+		this->inputs->AddInput(new TriaVertexInput(BathymetryEnum,nodeinputs));
 	}
 	if (iomodel->drag_coefficient) {
@@ -4172,4 +4220,19 @@
 }
 /*}}}*/
+/*FUNCTION Tria::IsNodeOnShelf {{{1*/
+bool   Tria::IsNodeOnShelf(){
+
+	int  i;
+	bool shelf=false;
+
+	for(i=0;i<3;i++){
+		if (nodes[i]->IsOnShelf()){
+			shelf=true;
+			break;
+		}
+	}
+	return shelf;
+}
+/*}}}*/
 /*FUNCTION Tria::IsOnWater {{{1*/
 bool   Tria::IsOnWater(){
@@ -4335,4 +4398,74 @@
 }
 /*}}}*/
+/*FUNCTION Tria::MigrateGroundingLine{{{1*/
+void  Tria::MigrateGroundingLine(void){
+
+	int     i;
+
+	double  h[3];
+	double  s[3];
+	double  b[3];
+	double  ba[3];
+	double  bed_hydro;
+	int     isonshelf[3];
+	Input  *surface_input  = NULL;
+	Input  *bed_input      = NULL;
+	double  rho_water;
+	double  rho_ice;
+	double  density;
+	double *values         = NULL;
+	bool    elementonshelf = false;
+
+
+	/*Recover info at the vertices: */
+	surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
+	bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
+	if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
+
+	GetParameterListOnVertices(&h[0],ThicknessEnum);
+	GetParameterListOnVertices(&s[0],SurfaceEnum);
+	GetParameterListOnVertices(&b[0],BedEnum);
+	GetParameterListOnVertices(&ba[0],BathymetryEnum);
+	for(i=0;i<3;i++)isonshelf[i]=nodes[i]->IsOnShelf();
+
+	/*material parameters: */
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+
+	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
+	for(i=0;i<3;i++){
+		if (isonshelf[i]){
+			/*This node is on the shelf. See if its bed is going under the bathymetry: */
+			if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway.
+				/*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */
+				b[i]=ba[i];
+				s[i]=b[i]+h[i];
+			}
+			else{
+				/*do nothing, we are still floating.*/
+			}
+		}
+		else{
+			/*This node is on the sheet, near the grounding line. See if wants to unground. To 
+			 * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */
+			bed_hydro=-density*h[i];
+			if (bed_hydro>ba[i]){
+				/*We are now floating, bed and surface are determined from hydrostatic equilibrium: */
+				s[i]=(1-density)*h[i];
+				b[i]=-density*h[i];
+			}
+			else{
+				/*do nothing, we are still grounded.*/
+			}
+		}
+	}
+
+	/*Surface and bed are updated. Update inputs:*/
+	surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i];
+	bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i];
+
+}
+/*}}}*/
 /*FUNCTION Tria::MinVel{{{1*/
 void  Tria::MinVel(double* pminvel, bool process_units){
@@ -4426,4 +4559,5 @@
 	int     numrows     = 0;
 	int     numnodes    = 0;
+	int     temp_numnodes    = 0;
 
 	/*Go through all the results objects, and update the counters: */
@@ -4433,5 +4567,6 @@
 		numrows++;
 		/*now, how many vertices and how many nodal values for this result? :*/
-		numnodes=elementresult->NumberOfNodalValues(); //ask result object.
+		temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
+		if(temp_numnodes>numnodes)numnodes=temp_numnodes;
 	}
 
@@ -5196,2 +5331,57 @@
 }
 /*}}}*/
+/*FUNCTION Tria::UpdateShelfStatus{{{1*/
+void  Tria::UpdateShelfStatus(void){
+
+	int     i;
+
+	double  h[3];
+	double  s[3];
+	double  b[3];
+	double  ba[3];
+	Input  *surface_input  = NULL;
+	Input  *bed_input      = NULL;
+	double  rho_water;
+	double  rho_ice;
+	double  density;
+	bool    elementonshelf = false;
+
+
+	/*Recover info at the vertices: */
+	surface_input   =inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
+	bed_input   =inputs->GetInput(BedEnum);   _assert_(bed_input);
+	if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!");
+
+	GetParameterListOnVertices(&h[0],ThicknessEnum);
+	GetParameterListOnVertices(&s[0],SurfaceEnum);
+	GetParameterListOnVertices(&b[0],BedEnum);
+	GetParameterListOnVertices(&ba[0],BathymetryEnum);
+
+	/*material parameters: */
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+
+	/*go through vertices, and figure out if they are grounded or not, then update their status: */
+	for(i=0;i<3;i++){
+		if(b[i]<=ba[i]){
+			nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
+			nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
+		}
+		else{
+			nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
+			nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
+		}
+	}
+
+	/*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
+	elementonshelf=true;
+	for(i=0;i<3;i++){
+		if(nodes[i]->IsOnSheet()){
+			elementonshelf=false;
+			break;
+		}
+	}
+    this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 7088)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 7089)
@@ -30,8 +30,10 @@
 
 		int  id;
+		int  sid;
 
 		Node   **nodes;    // 3 nodes
 		Matice  *matice;   // 1 material ice
 		Matpar  *matpar;   // 1 material parameter
+		int      horizontalneighborsids[3];
 
 		Parameters *parameters;   //pointer to solution parameters
@@ -41,5 +43,5 @@
 		/*Tria constructors, destructors {{{1*/
 		Tria();
-		Tria(int tria_id,int i, IoModel* iomodel,int nummodels);
+		Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels);
 		~Tria();
 		/*}}}*/
@@ -78,6 +80,8 @@
 		void   CreatePVector(Vec pg, Vec pf);
 		int    GetNodeIndex(Node* node);
+		int    Sid();
 		bool   IsOnBed();
 		bool   IsOnShelf(); 
+		bool   IsNodeOnShelf(); 
 		bool   IsOnWater(); 
 		void   GetSolutionFromInputs(Vec solution);
@@ -109,4 +113,5 @@
 		void   MaxVy(double* pmaxvy, bool process_units);
 		void   MaxVz(double* pmaxvz, bool process_units);
+		void   MigrateGroundingLine();
 		void   MinVel(double* pminvel, bool process_units);
 		void   MinVx(double* pminvx, bool process_units);
@@ -124,5 +129,7 @@
 		double SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
+		void   UpdateShelfStatus(void);
 		double TimeAdapt();
+		int*   GetHorizontalNeighboorSids(void);
 		/*}}}*/
 		/*Tria specific routines:{{{1*/
Index: /issm/trunk/src/c/objects/Inputs/BoolInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 7089)
@@ -148,6 +148,6 @@
 /*FUNCTION BoolInput::SpawnResult{{{1*/
 ElementResult* BoolInput::SpawnResult(int step, double time){
-
-	_error_(" not supported yet!");
+	
+	return new BoolElementResult(this->enum_type,this->value,step,time);
 
 }
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 7089)
@@ -399,6 +399,6 @@
 
 	*pvalues=this->values;
-	*pnum_values=3;
-
-}
-/*}}}*/
+	if(pnum_values)*pnum_values=3;
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/objects/IoModel.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/IoModel.cpp	(revision 7089)
@@ -38,4 +38,5 @@
 	xfree((void**)&this->z);
 	xfree((void**)&this->elements);
+	xfree((void**)&this->elementconnectivity);
 	xfree((void**)&this->elements_type);
 	xfree((void**)&this->vertices_type);
@@ -59,4 +60,5 @@
 	xfree((void**)&this->surface);
 	xfree((void**)&this->bed);
+	xfree((void**)&this->bathymetry);
 	xfree((void**)&this->vx_obs);
 	xfree((void**)&this->vy_obs);
@@ -191,4 +193,5 @@
 	IoModelFetchData(&this->waitonlock,iomodel_handle,"waitonlock");
 	IoModelFetchData(&this->kff,iomodel_handle,"kff");
+	IoModelFetchData(&this->grounding_line_migration,iomodel_handle,"grounding_line_migration");
 
 	/*!Get thermal parameters: */
@@ -214,4 +217,6 @@
 		IoModelFetchData(&this->qmu_save_femmodel,iomodel_handle,"qmu_save_femmodel"); 
 	}
+
+
 	/*i/o: */
 	IoModelFetchData(&this->io_gather,iomodel_handle,"io_gather");
@@ -242,4 +247,5 @@
 	this->z=NULL;
 	this->elements=NULL;
+	this->elementconnectivity=NULL;
 	this->elements_type=NULL;
 	this->vertices_type=NULL;
@@ -275,4 +281,5 @@
 	this->surface=NULL;
 	this->bed=NULL;
+	this->bathymetry=NULL;
 	this->elementoniceshelf=NULL;
 	this->elementonwater=NULL;
Index: /issm/trunk/src/c/objects/IoModel.h
===================================================================
--- /issm/trunk/src/c/objects/IoModel.h	(revision 7088)
+++ /issm/trunk/src/c/objects/IoModel.h	(revision 7089)
@@ -32,4 +32,5 @@
 		double* elements_type;
 		double* vertices_type;
+		double* elementconnectivity;
 
 		/*3d mesh: */
@@ -83,4 +84,5 @@
 		double* surface;
 		double* bed;
+		double* bathymetry;
 		double* elementoniceshelf;
 		double* elementonwater;
@@ -189,4 +191,7 @@
 		int      qmu_save_femmodel;
 
+		/*grounding line migration: */
+		int      grounding_line_migration;
+
 		/*exterior partitioning data, to be carried around: */
 		bool*   my_elements;
Index: /issm/trunk/src/c/objects/Loads/Riftfront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Riftfront.cpp	(revision 7088)
+++ /issm/trunk/src/c/objects/Loads/Riftfront.cpp	(revision 7089)
@@ -87,5 +87,5 @@
 	this->length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
 	this->fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
-	this->state=*(iomodel->riftinfo+RIFTINFOSIZE*i+11);
+	this->state=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+11);
 
 	//intialize inputs, and add as many inputs per element as requested: 
Index: /issm/trunk/src/c/objects/Node.h
===================================================================
--- /issm/trunk/src/c/objects/Node.h	(revision 7088)
+++ /issm/trunk/src/c/objects/Node.h	(revision 7089)
@@ -21,5 +21,5 @@
 class Node: public Object ,public Update{
 
-	private: 
+	public: 
 
 		int	    id;  //unique arbitrary id.
@@ -31,5 +31,4 @@
 		int            analysis_type;
 		
-	public:
 
 		/*Node constructors, destructors {{{1*/
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 7088)
+++ /issm/trunk/src/c/objects/objects.h	(revision 7089)
@@ -57,5 +57,6 @@
 #include "./ElementResults/DoubleElementResult.h"
 #include "./ElementResults/TriaVertexElementResult.h"
-#include "./ElementResults/PentaVertexElementResult.h"
+#include "./ElementResults/PentaVertexElementResult.h" 
+#include "./ElementResults/BoolElementResult.h"
 
 /*ExternalResults: */
Index: /issm/trunk/src/c/solutions/transient2d_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/transient2d_core.cpp	(revision 7088)
+++ /issm/trunk/src/c/solutions/transient2d_core.cpp	(revision 7089)
@@ -22,4 +22,5 @@
 	bool   time_adapt=false;
 	int    output_frequency;
+	bool   grounding_line_migration;
 
 	/*intermediary: */
@@ -35,4 +36,5 @@
 	femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum);
 	femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum);
+	femmodel->parameters->FindParam(&grounding_line_migration,GroundingLineMigrationEnum);
 
 	/*initialize: */
@@ -58,5 +60,10 @@
 		_printf_(VerboseSolution(),"%s\n","   computing new thickness");
 		prognostic_core(femmodel);
-	
+
+		if (grounding_line_migration){
+			_printf_(VerboseSolution(),"%s\n","   computing new grounding line position");
+			GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		}
+
 		if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){
 			_printf_(VerboseSolution(),"%s\n","   saving results\n");
@@ -68,4 +75,5 @@
 			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
 			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
+			if(grounding_line_migration)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time);
 
 			/*unload results*/