Index: /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp	(revision 16012)
+++ /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp	(revision 16013)
@@ -20,8 +20,9 @@
 
 	if(VerboseSolution()) _printf0_("call computational core:\n");
-	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
-	solutionsequence_linear(femmodel);
-	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
-	solutionsequence_linear(femmodel);
+	//femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
+	//solutionsequence_linear(femmodel);
+	//femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
+	//solutionsequence_linear(femmodel);
+	surfaceslope_core(femmodel);
 
 	if(VerboseSolution()) _printf0_("call computational core:\n");
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16012)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16013)
@@ -992,5 +992,5 @@
 	vector[1]=xyz_list[1][1] - xyz_list[0][1];
 
-	norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
+	norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
 
 	normal[0]= + vector[1]/norm;
@@ -7255,5 +7255,5 @@
 	GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
 	for(int i=0;i<numnodes;i++){
-		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]);
+		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
 		Nx[i] = Nx[i]/norm;
 		Ny[i] = Ny[i]/norm;
@@ -7261,5 +7261,5 @@
 
 	/*Build HNx and HNy*/
-	GetInputListOnNodes(H, ThicknessEnum);
+	GetInputListOnNodes(H,ThicknessEnum);
 	for(int i=0;i<numnodes;i++){
 		HNx[i]=H[i]*Nx[i];
@@ -7292,5 +7292,5 @@
 							basis[i]*basis[j]*(dhnx[0]+dhny[1])
 							+ basis[i]*(dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
-							+ gamma*(basis[j]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)*
+							+ gamma*(basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)*
 							(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny)
 							);
@@ -7490,5 +7490,5 @@
 	GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
 	for(int i=0;i<numnodes;i++){
-		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]);
+		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
 		Nx[i] = Nx[i]/norm;
 		Ny[i] = Ny[i]/norm;
@@ -7496,5 +7496,5 @@
 
 	/*Build HNx and HNy*/
-	GetInputListOnNodes(H, ThicknessEnum);
+	GetInputListOnNodes(H,ThicknessEnum);
 	for(int i=0;i<numnodes;i++){
 		HNx[i]=H[i]*Nx[i];