Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17324)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17325)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 1
+//#define FSANALYTICAL 3
 
 /*Model processing*/
@@ -413,16 +413,4 @@
 				_error_("not supported yet");
 			}
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					if(nodeonbed[i]>0. && groundedice_ls[i]>0. && nodeonFS[i]>0.){
-						if(vertices_type[i] == FSApproximationEnum){
-							for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
-						}
-						else{
-							_error_("not supported");
-						}
-					}
-				}
-			}
 			if(iomodel->meshtype==Mesh3DEnum){
 				IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17324)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17325)
@@ -2164,5 +2164,5 @@
 	int        numindices;
 	int       *indices = NULL;
-	IssmDouble slopex,slopey;
+	IssmDouble slopex,slopey,groundedice;
 	IssmDouble xz_plane[6];
 
@@ -2177,4 +2177,5 @@
 	Input* slopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(slopex_input);
 	Input* slopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(slopey_input);
+	Input* groundedicelevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedicelevelset_input);
 
 	/*Loop over basal nodes and update their CS*/
@@ -2185,11 +2186,23 @@
 		slopex_input->GetInputValue(&slopex,gauss);
 		slopey_input->GetInputValue(&slopey,gauss);
-
-		/*New X axis          New Z axis*/
-		xz_plane[0]=1.;       xz_plane[3]=-slopex;  
-		xz_plane[1]=0.;       xz_plane[4]=-slopey;  
-		xz_plane[2]=slopex;   xz_plane[5]=1.;          
+		groundedicelevelset_input->GetInputValue(&groundedice,gauss);
+
+		if(groundedice>0){
+			/*New X axis          New Z axis*/
+			xz_plane[0]=1.;       xz_plane[3]=-slopex;  
+			xz_plane[1]=0.;       xz_plane[4]=-slopey;  
+			xz_plane[2]=slopex;   xz_plane[5]=1.;          
+			this->nodes[indices[i]]->DofInSSet(2); //vz 
+		}
+		else{
+			/*New X axis          New Z axis*/
+			xz_plane[0]=1.;       xz_plane[3]=0;  
+			xz_plane[1]=0.;       xz_plane[4]=0;  
+			xz_plane[2]=0;        xz_plane[5]=1.;          
+			this->nodes[indices[i]]->DofInFSet(2); //vz
+		}
 
 		XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
+		
 	}