Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25450)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25451)
@@ -20,4 +20,5 @@
 //#define FSANALYTICAL 10
 //#define LATERALFRICTION 1
+//#define DISCSLOPE 1 //testing for SSA
 
 /*Model processing*/
@@ -1708,5 +1709,27 @@
 
 	/* Start  looping on the number of gaussian points: */
+	#ifndef DISCSLOPE
 	Gauss* gauss=element->NewGauss(2);
+	#else
+	Gauss* gauss=NULL;
+	Gauss* gauss_subelem=NULL;
+	bool partly_floating=false;
+	bool mainlyfloating=false;
+	int point1;
+	IssmDouble fraction1,fraction2;
+	IssmDouble phi=element->GetGroundedPortion(xyz_list);
+	if(phi>0.00000001 && phi<0.99999999) partly_floating=true; 
+
+	int ig=-1;// needed for driving stress parameterization
+	if(partly_floating){
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+	   gauss=element->NewGauss(point1,fraction1,fraction2,3); //considering the entire element
+		gauss_subelem=element->NewGauss(fraction1,fraction2,3);//gauss on each subelement
+	}
+	else{
+		gauss=element->NewGauss(2);//original
+	}
+	#endif
+
 	while(gauss->next()){
 
@@ -1716,4 +1739,13 @@
 		thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
 		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		#ifdef DISCSLOPE
+		if(gauss_subelem && partly_floating){
+			ig++;
+			gauss_subelem->next();
+			_assert_(std::abs(gauss_subelem->weight-gauss->weight)<0.0000001);
+			/*Compute the discontinuous surface slope for this gauss point/subelement*/
+			this->ComputeSurfaceSlope(&slope[0],gauss_subelem,gauss,point1,fraction1,fraction2,ig,dim,element);
+		} 		
+		#endif
 
 		for(int i=0;i<numnodes;i++){
@@ -1730,4 +1762,7 @@
 	xDelete<IssmDouble>(basis);
 	delete gauss;
+	#ifdef DISCSLOPE
+	if(gauss_subelem) delete gauss_subelem;
+	#endif
 	return pe;
 }/*}}}*/
@@ -2033,4 +2068,284 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+#ifdef DISCSLOPE
+void StressbalanceAnalysis::ComputeSurfaceSlope(IssmDouble* slope,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element){/*{{{*/
+
+	/*Compute the surface slope for each subelement, for a given integration point (gauss)*/
+	int numnodes=element->GetNumberOfNodes();
+	IssmDouble rho_ice=element->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_water=element->FindParam(MaterialsRhoSeawaterEnum); //ocean
+	IssmDouble* H=xNew<IssmDouble>(numnodes);
+	IssmDouble* S=xNew<IssmDouble>(numnodes);
+	IssmDouble* S_subelem=xNew<IssmDouble>(numnodes);
+	IssmDouble H_f1,H_f2,S_f1,S_f2;
+
+	/*Get nodal deriviatives of the subelements related to the Gauss point*/
+	IssmDouble* dbasis_subelem=xNew<IssmDouble>(dim*numnodes);//CG basis for each subelement
+	this->NodalFunctionsDerivativesRGB(dbasis_subelem,gauss_DG,gauss_CG,point1,fraction1,fraction2,ig,dim,element);
+
+	/*Define thickness at the grounding line (on element edges)*/
+	element->GetInputListOnNodes(H,ThicknessEnum); // thickness on vertices
+	switch(point1){//{{{
+		case 0:
+			H_f1=H[0]+(H[1]-H[0])*fraction1;
+			H_f2=H[0]+(H[2]-H[0])*fraction2;
+			break;
+		case 1:
+			H_f1=H[1]+(H[2]-H[1])*fraction1;
+			H_f2=H[1]+(H[0]-H[1])*fraction2;
+			break;
+		case 2:
+			H_f1=H[2]+(H[0]-H[2])*fraction1;
+			H_f2=H[2]+(H[1]-H[2])*fraction2;
+			break;
+		default:
+			_error_("index "<<point1<<" not supported yet");
+	}//}}}
+
+	/*Define surface at the grounding (on element edges)*/
+	S_f1=H_f1*(1-rho_ice/rho_water);
+	S_f2=H_f2*(1-rho_ice/rho_water);
+	element->GetInputListOnNodes(S,SurfaceEnum); // surface on vertices
+
+	/*Define surface on the subelement vertices*/
+   if(ig<4){ // BLUE element itapopo only if order is = 3
+		switch(point1){ //{{{
+			case 0:
+				S_subelem[0]=S[0];
+				S_subelem[1]=S_f1;
+				S_subelem[2]=S_f2;
+				break;
+			case 1:
+				S_subelem[0]=S[1];
+				S_subelem[1]=S_f1;
+				S_subelem[2]=S_f2;
+				break;
+			case 2:
+				S_subelem[0]=S[2];
+				S_subelem[1]=S_f1;
+				S_subelem[2]=S_f2;
+				break;
+			default:
+				_error_("index "<<point1<<" not supported yet");
+		}//}}}	
+	}
+	if(ig>3 && ig<8){ // GREEN element
+		switch(point1){ //{{{
+			case 0:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[2];
+				S_subelem[2]=S_f2;
+				break;
+			case 1:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[0];
+				S_subelem[2]=S_f2;
+				break;
+			case 2:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[1];
+				S_subelem[2]=S_f2;
+				break;
+			default:
+				_error_("index "<<point1<<" not supported yet");
+		}//}}}	
+	}
+	if(ig>7){ // RED element
+		switch(point1){ //{{{
+			case 0:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[1];
+				S_subelem[2]=S[2];
+				break;
+			case 1:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[2];
+				S_subelem[2]=S[0];
+				break;
+			case 2:
+				S_subelem[0]=S_f1;
+				S_subelem[1]=S[0];
+				S_subelem[2]=S[1];
+				break;
+			default:
+				_error_("index "<<point1<<" not supported yet");
+		}//}}}	
+	}
+
+	//_printf_("\t"<<"H[1]\t"<<"H[2]\t"<<"H[3]\t"<<"S[1]\t"<<"S[2]\t"<<"S[3]\n");
+	//_printf_("\t"<<H[0]<<"\t"<<H[1]<<"\t"<<H[2]<<"\t"<<S[0]<<"\t"<<S[1]<<"\t"<<S[2]<<"\n");
+	//_printf_("\t H_f1 \t H_f2 \t S_f1 \t S_f2 \n");
+	//_printf_("\t"<< H_f1 <<"\t"<< H_f2 << "\t" << S_f1 << "\t" << S_f2<< "\n");
+	//_printf_("\t"<<"S_subelem[1]\t"<<"S_subelem[2]\t"<<"S_subelem[3]\n");
+	//_printf_("\t"<<S_subelem[0]<<"\t"<<S_subelem[1]<<"\t"<<S_subelem[2]<<"\n");
+
+	/*Compute slope*/
+	slope[0]=0;
+	slope[1]=0;
+	for(int i=0;i<numnodes;i++){
+		slope[0]+=S_subelem[i]*dbasis_subelem[numnodes*0+i]; //dSdx
+		slope[1]+=S_subelem[i]*dbasis_subelem[numnodes*1+i]; //dSdy
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(dbasis_subelem);
+	xDelete<IssmDouble>(H);
+	xDelete<IssmDouble>(S);
+	xDelete<IssmDouble>(S_subelem);
+
+}/*}}}*/
+void StressbalanceAnalysis::NodalFunctionsDerivativesRGB(IssmDouble* dbasis_subelem,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element){/*{{{*/
+	
+	/*Fetch number of nodes for this finite element*/
+   int numnodes = element->GetNumberOfNodes();
+   IssmDouble dbasis_ref[dim*numnodes];
+   IssmDouble Jinv[2][2];
+   //IssmDouble* dbasis_subelem = xNew<IssmDouble>(dim*numnodes);//CG basis for each subelement
+   IssmDouble* xyz_list_RGB = xNew<IssmDouble>(3*numnodes); // x,y,z per node
+   IssmDouble* xyz_list = NULL;
+
+   element->GetVerticesCoordinates(&xyz_list);
+   Tria* tria = xDynamicCast<Tria*>(element);
+
+   /*Get nodal functions derivatives*/
+   tria->GetNodalFunctionsDerivativesReference(dbasis_ref,gauss_DG,P1Enum); //gauss is not used here
+
+   if(ig<4){ // BLUE element itapopo only if order is = 3
+      switch(point1){//{{{
+         case 0:
+            /*Vertex 0 - local in the subelement*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]+(xyz_list[3*2+0]-xyz_list[3*0+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]+(xyz_list[3*2+1]-xyz_list[3*0+1])*fraction2; //y;
+            break;
+			 case 1:
+            /*Vertex 0 - local in the subelement*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]+(xyz_list[3*0+0]-xyz_list[3*1+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]+(xyz_list[3*0+1]-xyz_list[3*1+1])*fraction2; //y;
+            break;
+			case 2:
+            /*Vertex 0 - local in the subelement*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]+(xyz_list[3*1+0]-xyz_list[3*2+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]+(xyz_list[3*1+1]-xyz_list[3*2+1])*fraction2; //y;
+            break;
+         default:
+            _error_("index "<<point1<<" not supported yet");
+      }//}}}
+   }
+	if(ig>3 && ig<8){ // GREEN element
+      switch(point1){//{{{
+         case 0:
+            /*Vertex 0 - local in the subelement*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]+(xyz_list[3*2+0]-xyz_list[3*0+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]+(xyz_list[3*2+1]-xyz_list[3*0+1])*fraction2; //y;
+            break;
+			 case 1:
+            /*Vertex 0*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]+(xyz_list[3*0+0]-xyz_list[3*1+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]+(xyz_list[3*0+1]-xyz_list[3*1+1])*fraction2; //y;
+            break;
+			case 2:
+            /*Vertex 0*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]+(xyz_list[3*1+0]-xyz_list[3*2+0])*fraction2; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]+(xyz_list[3*1+1]-xyz_list[3*2+1])*fraction2; //y;
+            break;
+         default:
+            _error_("index "<<point1<<" not supported yet");
+      }//}}}
+   }
+	if(ig>7){ // RED element
+      switch(point1){//{{{
+         case 0:
+            /*Vertex 0*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]; //y;
+            break;
+         case 1:
+				/*Vertex 0*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]; //y;
+            break;
+			case 2:
+            /*Vertex 0*/
+            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
+            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
+            /*Vertex 1*/
+            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]; //x
+            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]; //y
+            /*Vertex 2*/
+            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]; //x;
+            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]; //y;
+            break;
+         default:
+            _error_("index "<<point1<<" not supported yet");
+      }//}}}
+   }
+
+	//std::cout<<"	ig="<<ig<<std::endl;
+	//std::cout<<"	"<<xyz_list_RGB[0*numnodes+0]<<"\t"<<xyz_list_RGB[1*numnodes+0]<<"\t"<<xyz_list_RGB[2*numnodes+0]<<std::endl; //x
+	//std::cout<<"	"<<xyz_list_RGB[0*numnodes+1]<<"\t"<<xyz_list_RGB[1*numnodes+1]<<"\t"<<xyz_list_RGB[2*numnodes+1]<<std::endl; //y
+
+   /*Get Jacobian*/
+	tria->GetJacobianInvert(&Jinv[0][0],xyz_list_RGB,gauss_CG); // gauss is not used here
+	
+	/*Build dbasis CG in the subelement*/
+	for(int i=0;i<numnodes;i++){
+		dbasis_subelem[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];//dPhi_i/dx
+		dbasis_subelem[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];//dPhi_i/dy
+	}
+
+   xDelete<IssmDouble>(xyz_list);
+   xDelete<IssmDouble>(xyz_list_RGB);
+   //xDelete<IssmDouble>(dbasis_subelem);
+
+}/*}}}*/
+#endif
 
 /*L1L2*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 25450)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 25451)
@@ -45,4 +45,7 @@
 		void           GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void           InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+		void				ComputeSurfaceSlope(IssmDouble* slope,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element);
+		void				NodalFunctionsDerivativesRGB(IssmDouble* dbasis_subelem,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element);
+
 		/*L1L2*/
 		ElementMatrix* CreateKMatrixL1L2(Element* element);
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp	(revision 25450)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp	(revision 25451)
@@ -366,5 +366,5 @@
 	 *
 	 */
-	int         ig;
+	int         ii;
 	IssmDouble x,y;
 	IssmDouble xy_list[3][2];
@@ -384,21 +384,21 @@
 	this->weights=xNew<IssmDouble>(this->numgauss);
 
-	for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points (BLUE)
-		this->coords1[ig]=gauss1->coords1[ig];
-		this->coords2[ig]=gauss1->coords2[ig];
-		this->coords3[ig]=gauss1->coords3[ig];
-		this->weights[ig]=gauss1->weights[ig]*r1*r2;
-	}
-	for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points (GREEN)
-		this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
-		this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
-		this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
-		this->weights[gauss1->numgauss+ig]=gauss2->weights[ig]*r1*(1-r2);
-	}
-	for(ig=0;ig<gauss3->numgauss;ig++){ // Add the second triangle gauss points (RED)
-		this->coords1[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords1[ig];
-		this->coords2[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords2[ig];
-		this->coords3[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords3[ig];
-		this->weights[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->weights[ig]*(1-r1);
+	for(ii=0;ii<gauss1->numgauss;ii++){ // Add the first triangle gauss points (BLUE)
+		this->coords1[ii]=gauss1->coords1[ii];
+		this->coords2[ii]=gauss1->coords2[ii];
+		this->coords3[ii]=gauss1->coords3[ii];
+		this->weights[ii]=gauss1->weights[ii]*r1*r2;
+	}
+	for(ii=0;ii<gauss2->numgauss;ii++){ // Add the second triangle gauss points (GREEN)
+		this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii];
+		this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii];
+		this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii];
+		this->weights[gauss1->numgauss+ii]=gauss2->weights[ii]*r1*(1-r2);
+	}
+	for(ii=0;ii<gauss3->numgauss;ii++){ // Add the second triangle gauss points (RED)
+		this->coords1[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords1[ii];
+		this->coords2[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords2[ii];
+		this->coords3[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords3[ii];
+		this->weights[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->weights[ii]*(1-r1);
 	}