Index: issm/trunk/src/c/analyses/FreeSurfaceTopAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 26744)
+++ issm/trunk/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 27035)
@@ -80,4 +80,5 @@
 	iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
 	if(iomodel->domaintype!=Domain2DhorizontalEnum){
 		iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
@@ -111,5 +112,5 @@
 }/*}}}*/
 ElementMatrix* FreeSurfaceTopAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
-_error_("Not implemented");
+	_error_("Not implemented");
 }/*}}}*/
 ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
@@ -120,5 +121,5 @@
 	IssmDouble *xyz_list  = NULL;
 	IssmDouble  Jdet,D_scalar,dt,h;
-	IssmDouble  vel,vx,vy;
+	IssmDouble  vel,vx,vy,tau;
 
 	/*Get top element*/
@@ -195,5 +196,19 @@
 		}
 
-		if(stabilization==2){
+		if(stabilization==1){
+			/*artifical diffusion*/
+			if(dim==1){
+				vx_input->GetInputAverage(&vx);
+				D[0]=h/2.*fabs(vx);
+			}
+			else{
+				vx_input->GetInputAverage(&vx);
+				vy_input->GetInputAverage(&vy);
+
+				D[0*dim+0]=h/2.0*fabs(vx);
+				D[1*dim+1]=h/2.0*fabs(vy);
+			}
+		}
+		else if(stabilization==2){
 			/*Streamline upwinding*/
 			if(dim==1){
@@ -209,16 +224,14 @@
 			}
 		}
-		else if(stabilization==1){
-			/*SSA*/
+		else if(stabilization==5){
+			/*SUPG*/
 			if(dim==1){
 				vx_input->GetInputAverage(&vx);
-				D[0]=h/2.*fabs(vx);
+				tau=h/(2.*fabs(vx)+1e-10);
 			}
 			else{
 				vx_input->GetInputAverage(&vx);
 				vy_input->GetInputAverage(&vy);
-
-				D[0*dim+0]=h/2.0*fabs(vx);
-				D[1*dim+1]=h/2.0*fabs(vy);
+				tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
 			}
 		}
@@ -229,7 +242,7 @@
 					for(int j=0;j<numnodes;j++){
 						Ke->values[i*numnodes+j] += (
-									dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
-									dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
-									);
+								dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
+								dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
+								);
 					}
 				}
@@ -237,4 +250,19 @@
 			else{
 				for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
+			}
+		}
+		else if(stabilization==5){
+			D_scalar=gauss->weight*Jdet*dt;
+			if(dim==2){
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=tau*D_scalar*
+							(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
+							(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
+					}
+				}
+			}
+			else{
+				for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
 			}
 		}
@@ -251,8 +279,9 @@
 }/*}}}*/
 ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
+
 	/*Intermediaries*/
-	int         domaintype,dim;
+	int         domaintype,dim,stabilization;
 	IssmDouble  Jdet,dt;
-	IssmDouble  ms,surface,vz;
+	IssmDouble  ms,surface,vx,vy,vz,tau;
 	Element*    topelement = NULL;
 	IssmDouble *xyz_list  = NULL;
@@ -284,16 +313,27 @@
 	ElementVector* pe    = topelement->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
 	topelement->GetVerticesCoordinates(&xyz_list);
 	topelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* ms_input      = topelement->GetInput(SmbMassBalanceEnum);  _assert_(ms_input);
-	Input* surface_input = topelement->GetInput(SurfaceEnum);                     _assert_(surface_input);
-	Input* vz_input      = NULL;
+	Input *ms_input      = topelement->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
+	Input *surface_input = topelement->GetInput(SurfaceEnum);        _assert_(surface_input);
+	Input *vz_input      = NULL;
+	Input *vx_input      = NULL;
+	Input *vy_input      = NULL;
 	switch(dim){
-		case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break;
-		case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break;
+		case 1: 
+			vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
+			vz_input = topelement->GetInput(VyEnum)	; _assert_(vz_input); 
+			break;
+		case 2: 
+			vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
+			vy_input = topelement->GetInput(VyEnum); _assert_(vy_input);
+			vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); 
+			break;
 		default: _error_("not implemented");
 	}
+	IssmDouble h = topelement->CharacteristicLength();
 
 	/*Initialize mb_correction to 0, do not forget!:*/
@@ -312,7 +352,24 @@
 	}
 
+	if(stabilization==5){
+		/*SUPG*/
+		topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		if(dim==1){
+			vx_input->GetInputAverage(&vx);
+			tau=h/(2.*fabs(vx)+1e-10);
+			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
+		}
+		else{
+			vx_input->GetInputAverage(&vx);
+			vy_input->GetInputAverage(&vy);
+			tau=h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
+			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+		}
+	}
+
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
@@ -321,5 +378,5 @@
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-	   _error_("not implemented yet");
+	_error_("not implemented yet");
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
@@ -329,4 +386,71 @@
 
 	element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
+
+	/*Now, we need to do some "processing"*/
+	int numvertices = element->GetNumberOfVertices();
+	int        migration_style;
+
+	IssmDouble* surface = xNew<IssmDouble>(numvertices);
+	IssmDouble* newsurface = xNew<IssmDouble>(numvertices);
+	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+	IssmDouble* base = xNew<IssmDouble>(numvertices);
+	IssmDouble* bed = xNew<IssmDouble>(numvertices);
+	IssmDouble* phi = xNew<IssmDouble>(numvertices);
+	IssmDouble* sealevel = xNew<IssmDouble>(numvertices);
+
+	IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
+	IssmDouble rho_ice      = element->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_water    = element->FindParam(MaterialsRhoSeawaterEnum);
+
+	bool isgroundingline;
+	element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
+	element->FindParam(&migration_style,GroundinglineMigrationEnum);
+	if(isgroundingline) element->GetInputListOnVertices(&bed[0],BedEnum);
+
+	element->GetInputListOnVertices(&base[0],BaseEnum);
+	element->GetInputListOnVertices(&surface[0],SurfaceEnum);
+	element->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
+	element->GetInputListOnVertices(&sealevel[0],SealevelEnum);
+
+	for(int i=0;i<numvertices;i++){
+		newsurface[i]=surface[i];
+		thickness[i]=surface[i]-base[i];
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(thickness[i])) _error_("NaN found in solution vector");
+		if(xIsInf<IssmDouble>(thickness[i])) _error_("Inf found in solution vector");
+
+		/* check for thickness<minthickness */
+		if(thickness[i]<minthickness){ 
+			thickness[i]=minthickness;
+			if(phi[i]>0.){
+				if(base[i]<=bed[i]) base[i] = bed[i];
+				newsurface[i] = base[i]+minthickness;
+			}else{
+				// assume floatation condition
+				newsurface[i] = (1.-rho_ice/rho_water)*minthickness;
+				base[i] = -rho_ice/rho_water*minthickness;
+			}
+		}
+
+		/* update thickness */
+		thickness[i]=newsurface[i]-base[i];
+		/* some checks */
+		if(thickness[i]<0.) _error_("thickness<0");
+		if(newsurface[i]<base[i]) _error_("surface<base");
+	}
+
+	/* update inputs */
+	element->AddInput(BaseEnum,base,element->GetElementType());
+	element->AddInput(SurfaceEnum,newsurface,element->GetElementType());
+	element->AddInput(ThicknessEnum,thickness,element->GetElementType());
+
+	/* Free resources */
+	xDelete<IssmDouble>(newsurface);
+	xDelete<IssmDouble>(surface);
+	xDelete<IssmDouble>(thickness);
+	xDelete<IssmDouble>(base);
+	xDelete<IssmDouble>(bed);
+	xDelete<IssmDouble>(phi);
+	xDelete<IssmDouble>(sealevel);
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
