Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 16434)
@@ -140,4 +140,174 @@
 /*}}}*/
 
+/*FUNCTION Seg::GetSize{{{*/
+IssmDouble Seg::GetSize(void){
+
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble x1,y1,x2,y2;
+
+	/*Get xyz list: */
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	x1=xyz_list[0][0]; y1=xyz_list[0][1];
+	x2=xyz_list[1][0]; y2=xyz_list[1][1];
+
+	return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
+}
+/*}}}*/
+/*FUNCTION Seg::CreateKMatrixFreeSurfaceTop {{{*/
+ElementMatrix* Seg::CreateKMatrixFreeSurfaceTop(void){
+
+	/*Intermediaries */
+	int        stabilization;
+	IssmDouble Jdet,D_scalar,dt,h;
+	IssmDouble vx,vel;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	h=this->GetSize();
+
+	/* Start  looping on the number of gaussian points: */
+	GaussSeg *gauss=new GaussSeg(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+
+		D_scalar=gauss->weight*Jdet;
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D_scalar,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+
+		GetNodalFunctions(B,gauss);
+		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
+
+		D_scalar=dt*gauss->weight*Jdet*vx;
+		TripleMultiply(B,1,numnodes,1,
+					&D_scalar,1,1,0,
+					Bprime,1,numnodes,0,
+					&Ke->values[0],1);
+
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=fabs(vx)+1.e-8;
+			D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
+		}
+		else if(stabilization==1){
+			/*SSA*/
+			vx_input->GetInputAverage(&vx);
+			D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
+		}
+		if(stabilization==1 || stabilization==2){
+			TripleMultiply(Bprime,1,numnodes,1,
+						&D_scalar,1,1,0,
+						Bprime,1,numnodes,0,
+						&Ke->values[0],1);
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Seg::CreateKMatrixFreeSurfaceBase {{{*/
+ElementMatrix* Seg::CreateKMatrixFreeSurfaceBase(void){
+
+	/*Intermediaries */
+	int        stabilization;
+	IssmDouble Jdet,D_scalar,dt,h;
+	IssmDouble vx,vel;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	h=this->GetSize();
+
+	/* Start  looping on the number of gaussian points: */
+	GaussSeg *gauss=new GaussSeg(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+
+		D_scalar=gauss->weight*Jdet;
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D_scalar,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+
+		GetNodalFunctions(B,gauss);
+		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
+
+		D_scalar=dt*gauss->weight*Jdet*vx;
+		TripleMultiply(B,1,numnodes,1,
+					&D_scalar,1,1,0,
+					Bprime,1,numnodes,0,
+					&Ke->values[0],1);
+
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=fabs(vx)+1.e-8;
+			D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
+		}
+		else if(stabilization==1){
+			/*SSA*/
+			vx_input->GetInputAverage(&vx);
+			D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
+		}
+		if(stabilization==1 || stabilization==2){
+			TripleMultiply(Bprime,1,numnodes,1,
+						&D_scalar,1,1,0,
+						Bprime,1,numnodes,0,
+						&Ke->values[0],1);
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
 /*FUNCTION Seg::CreateMassMatrix {{{*/
 ElementMatrix* Seg::CreateMassMatrix(void){
@@ -229,2 +399,97 @@
 }
 /*}}}*/
+/*FUNCTION Seg::CreatePVectorFreeSurfaceTop {{{*/
+ElementVector* Seg::CreatePVectorFreeSurfaceTop(void){
+
+	/*Intermediaries */
+	IssmDouble Jdet,dt;
+	IssmDouble ms,surface,vy;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vy_input     = inputs->GetInput(VyEnum);                         _assert_(vy_input);
+	Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
+	Input* surface_input= inputs->GetInput(SurfaceEnum);                    _assert_(surface_input);
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* Start  looping on the number of gaussian points: */
+	GaussSeg* gauss=new GaussSeg(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		vy_input->GetInputValue(&vy,gauss);
+		ms_input->GetInputValue(&ms,gauss);
+		surface_input->GetInputValue(&surface,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vy)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Seg::CreatePVectorFreeSurfaceBase {{{*/
+ElementVector* Seg::CreatePVectorFreeSurfaceBase(void){
+
+	/*Intermediaries */
+	IssmDouble Jdet,dt;
+	IssmDouble mb,mb_correction,bed,vy;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vy_input            = inputs->GetInput(VyEnum);                         _assert_(vy_input);
+	Input* mb_input            = inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
+	Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+	Input* bed_input           = inputs->GetInput(BedEnum);                        _assert_(bed_input);
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* Start  looping on the number of gaussian points: */
+	GaussSeg* gauss=new GaussSeg(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		vy_input->GetInputValue(&vy,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		bed_input->GetInputValue(&bed,gauss);
+		if(mb_correction_input)
+		 mb_correction_input->GetInputValue(&mb_correction,gauss);
+		else
+		 mb_correction=0.;
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vy)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16434)
@@ -192,4 +192,9 @@
 		/*Seg specific routines:*/
 		ElementMatrix* CreateMassMatrix(void);
+		ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
+		ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
+		ElementVector* CreatePVectorFreeSurfaceTop(void);
+		ElementVector* CreatePVectorFreeSurfaceBase(void);
+		IssmDouble     GetSize(void);
 };
 #endif  /* _SEG_H */
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 16434)
@@ -48,4 +48,31 @@
 
 /*Reference Element numerics*/
+/*FUNCTION SegRef::GetBprimeMasstransport{{{*/
+void SegRef::GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussSeg* gauss){
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ dN/dx ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		Bprime[i] = dbasis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}
+/*}}}*/
 /*FUNCTION SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){{{*/
 void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 16434)
@@ -22,4 +22,5 @@
 		/*Management*/
 		void SetElementType(int type,int type_counter);
+		void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussSeg* gauss);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16434)
@@ -208,4 +208,5 @@
 	/*retreive parameters: */
 	int analysis_type;
+	int meshtype;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
@@ -246,5 +247,4 @@
 		#endif
 		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
-			int meshtype;
 			parameters->FindParam(&meshtype,MeshTypeEnum);
 			if(meshtype==Mesh2DverticalEnum){
@@ -258,4 +258,26 @@
 		case MasstransportAnalysisEnum:
 			return CreateKMatrixMasstransport();
+			break;
+		case FreeSurfaceTopAnalysisEnum:
+			parameters->FindParam(&meshtype,MeshTypeEnum);
+			switch(meshtype){
+				case Mesh2DverticalEnum:
+					return CreateKMatrixFreeSurfaceTop1D();
+				case Mesh3DEnum:
+					return CreateKMatrixFreeSurfaceTop();
+				default:
+					_error_("Mesh not supported yet");
+			}
+			break;
+		case FreeSurfaceBaseAnalysisEnum:
+			parameters->FindParam(&meshtype,MeshTypeEnum);
+			switch(meshtype){
+				case Mesh2DverticalEnum:
+					return CreateKMatrixFreeSurfaceBase1D();
+				case Mesh3DEnum:
+					return CreateKMatrixFreeSurfaceBase();
+				default:
+					_error_("Mesh not supported yet");
+			}
 			break;
 		#endif
@@ -396,4 +418,5 @@
 
 	/*retrive parameters: */
+	int meshtype;
 	int analysis_type;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
@@ -436,5 +459,4 @@
 		#endif
 		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
-			int meshtype;
 			parameters->FindParam(&meshtype,MeshTypeEnum);
 			if(meshtype==Mesh2DverticalEnum){
@@ -448,4 +470,26 @@
 		case MasstransportAnalysisEnum:
 			return CreatePVectorMasstransport();
+			break;
+		case FreeSurfaceTopAnalysisEnum:
+			parameters->FindParam(&meshtype,MeshTypeEnum);
+			switch(meshtype){
+				case Mesh2DverticalEnum:
+					return CreatePVectorFreeSurfaceTop1D();
+				case Mesh3DEnum:
+					return CreatePVectorFreeSurfaceTop();
+				default:
+					_error_("Mesh not supported yet");
+			}
+			break;
+		case FreeSurfaceBaseAnalysisEnum:
+			parameters->FindParam(&meshtype,MeshTypeEnum);
+			switch(meshtype){
+				case Mesh2DverticalEnum:
+					return CreatePVectorFreeSurfaceBase1D();
+				case Mesh3DEnum:
+					return CreatePVectorFreeSurfaceBase();
+				default:
+					_error_("Mesh not supported yet");
+			}
 			break;
 		#endif
@@ -833,5 +877,4 @@
 IssmDouble Tria::GetArea(void){
 
-	IssmDouble area=0;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble x1,y1,x2,y2,x3,y3;
@@ -1822,4 +1865,10 @@
 			InputUpdateFromSolutionMasstransport(solution);
 			break;
+		case FreeSurfaceTopAnalysisEnum:
+			InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
+			break;
+		case FreeSurfaceBaseAnalysisEnum:
+			InputUpdateFromSolutionOneDof(solution,BedEnum);
+			break;
 		default:
 			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -2163,4 +2212,26 @@
 }
 /*}}}*/
+/*FUNCTION Tria::HasEdgeOnSurface {{{*/
+bool Tria::HasEdgeOnSurface(){
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum==3.)  _error_("Two edges on surface not supported yet...");
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}
+/*}}}*/
 /*FUNCTION Tria::EdgeOnBedIndices{{{*/
 void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
@@ -2184,4 +2255,25 @@
 }
 /*}}}*/
+/*FUNCTION Tria::EdgeOnSurfaceIndices{{{*/
+void Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			*pindex1 = indices[i][0];
+			*pindex2 = indices[i][1];
+			return;
+		}
+	}
+
+	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on surface");
+}
+/*}}}*/
 /*FUNCTION Tria::EdgeOnBedIndex{{{*/
 int Tria::EdgeOnBedIndex(void){
@@ -2201,4 +2293,23 @@
 	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
 	_error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
+/*FUNCTION Tria::EdgeOnSurfaceIndex{{{*/
+int Tria::EdgeOnSurfaceIndex(void){
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			return i;
+		}
+	}
+
+	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on surface");
 }
 /*}}}*/
@@ -7450,4 +7561,20 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop1D {{{*/
+ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop1D(void){
+
+	if(!HasEdgeOnSurface()) return NULL;
+
+	int index1,index2;
+	this->EdgeOnSurfaceIndices(&index1,&index2);
+
+	Seg* seg=(Seg*)SpawnSeg(index1,index2); 
+	ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceTop();
+	delete seg->material; delete seg;
+
+	/*clean up and return*/
+	return Ke;
+}
+/*}}}*/
 /*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/
 ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){
@@ -7547,4 +7674,20 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase1D {{{*/
+ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase1D(void){
+
+	if(!HasEdgeOnBed()) return NULL;
+
+	int index1,index2;
+	this->EdgeOnBedIndices(&index1,&index2);
+
+	Seg* seg=(Seg*)SpawnSeg(index1,index2); 
+	ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceBase();
+	delete seg->material; delete seg;
+
+	/*clean up and return*/
+	return Ke;
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorMasstransport{{{*/
 ElementVector* Tria::CreatePVectorMasstransport(void){
@@ -7699,4 +7842,20 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreatePVectorFreeSurfaceTop1D {{{*/
+ElementVector* Tria::CreatePVectorFreeSurfaceTop1D(void){
+
+	if(!HasEdgeOnSurface()) return NULL;
+
+	int index1,index2;
+	this->EdgeOnSurfaceIndices(&index1,&index2);
+
+	Seg* seg=(Seg*)SpawnSeg(index1,index2); 
+	ElementVector* pe=seg->CreatePVectorFreeSurfaceTop();
+	delete seg->material; delete seg;
+
+	/*clean up and return*/
+	return pe;
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/
 ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){
@@ -7746,4 +7905,20 @@
 	xDelete<IssmDouble>(basis);
 	delete gauss;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorFreeSurfaceBase1D {{{*/
+ElementVector* Tria::CreatePVectorFreeSurfaceBase1D(void){
+
+	if(!HasEdgeOnBed()) return NULL;
+
+	int index1,index2;
+	this->EdgeOnBedIndices(&index1,&index2);
+
+	Seg* seg=(Seg*)SpawnSeg(index1,index2); 
+	ElementVector* pe=seg->CreatePVectorFreeSurfaceBase();
+	delete seg->material; delete seg;
+
+	/*clean up and return*/
 	return pe;
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16434)
@@ -86,6 +86,9 @@
 		bool        IsOnBed();
 		bool        HasEdgeOnBed();
+		bool        HasEdgeOnSurface();
+		void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
 		void        EdgeOnBedIndices(int* pindex1,int* pindex);
 		int         EdgeOnBedIndex();
+		int         EdgeOnSurfaceIndex();
 		bool        IsFloating();
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
@@ -191,5 +194,7 @@
 		ElementMatrix* CreateKMatrixMasstransport_DG(void);
 		ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
+		ElementMatrix* CreateKMatrixFreeSurfaceTop1D(void);
 		ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
+		ElementMatrix* CreateKMatrixFreeSurfaceBase1D(void);
 		ElementMatrix* CreateMassMatrix(void);
 		ElementMatrix* CreateBasalMassMatrix(void);
@@ -205,5 +210,7 @@
 		ElementVector* CreatePVectorMasstransport_DG(void);
 		ElementVector* CreatePVectorFreeSurfaceTop(void);
+		ElementVector* CreatePVectorFreeSurfaceTop1D(void);
 		ElementVector* CreatePVectorFreeSurfaceBase(void);
+		ElementVector* CreatePVectorFreeSurfaceBase1D(void);
 		ElementVector* CreatePVectorSlope(void);
 		ElementVector* CreateBasalPVectorSlope(void);
Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16434)
@@ -47,4 +47,5 @@
 		void GetInputValue(int* pvalue);
 		void GetInputValue(IssmDouble* pvalue);
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16434)
@@ -52,4 +52,5 @@
 		void GetInputValue(int* pvalue);
 		void GetInputValue(IssmDouble* pvalue);
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16434)
@@ -49,4 +49,5 @@
 		void GetInputValue(int* pvalue){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16434)
@@ -46,4 +46,5 @@
 		void GetInputValue(int* pvalue);
 		void GetInputValue(IssmDouble* pvalue);
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/Input.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16434)
@@ -29,4 +29,5 @@
 		virtual void GetInputValue(int* pvalue)=0;
 		virtual void GetInputValue(IssmDouble* pvalue)=0;
+		virtual void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss)=0;
 		virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss)=0;
 		virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0;
Index: /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16434)
@@ -47,4 +47,5 @@
 		void GetInputValue(int* pvalue);
 		void GetInputValue(IssmDouble* pvalue);
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16434)
@@ -46,4 +46,5 @@
 		void GetInputValue(int* pvalue){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 16434)
@@ -85,4 +85,17 @@
 
 /*Object functions*/
+/*FUNCTION SegInput::GetInputAverage{{{*/
+void SegInput::GetInputAverage(IssmDouble* pvalue){
+
+	int        numnodes  = this->NumberofNodes();
+	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
+	IssmDouble value     = 0.;
+
+	for(int i=0;i<numnodes;i++) value+=values[i];
+	value = value/numnodesd;
+
+	*pvalue=value;
+}
+/*}}}*/
 /*FUNCTION SegInput::GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){{{*/
 void SegInput::GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 16434)
@@ -57,5 +57,5 @@
 		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
 		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
-		void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");};
+		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
 		void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16434)
@@ -50,4 +50,5 @@
 		void GetInputValue(int* pvalue){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16434)
@@ -97,5 +97,5 @@
 }
 /*}}}*/
-/*FUNCTION SegInput::SpawnSegInput{{{*/
+/*FUNCTION TriaInput::SpawnSegInput{{{*/
 Input* TriaInput::SpawnSegInput(int index1,int index2){
 
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16433)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16434)
@@ -47,4 +47,5 @@
 		void GetInputValue(int* pvalue){_error_("not implemented yet");}
 		void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}
+		void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp	(revision 16433)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp	(revision 16434)
@@ -30,5 +30,5 @@
 	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.);
 	iomodel->FetchDataToInput(elements,VxEnum);
-	iomodel->FetchDataToInput(elements,VyEnum);
+	iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
 	if(iomodel->meshtype==Mesh3DEnum){
 		iomodel->FetchDataToInput(elements,VzEnum);