Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8415)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8416)
@@ -1628,10 +1628,48 @@
 	if (!IsOnSurface()) return NULL;
 
-	/*Call Tria function*/
-	Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
-	ElementMatrix* Ke=tria->CreateKMatrixDiagnosticVertSurface();
-	delete tria->matice; delete tria;
-
-	/*clean up and return*/
+	/*Constants*/
+	const int numdof=NDOF1*NUMVERTICES;
+
+	/*Intermediaries */
+	int       i,j,ig;
+	double    xyz_list[NUMVERTICES][3];
+	double    xyz_list_tria[NUMVERTICES2D][3];
+	double    surface_normal[3];
+	double    Jdet2d,DL_scalar;
+	double    l1l6[NUMVERTICES];
+	GaussPenta *gauss=NULL;
+	double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j];
+	SurfaceNormal(&surface_normal[0],xyz_list_tria);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(3,4,5,2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
+		GetNodalFunctionsP1(&l1l6[0], gauss);
+
+		/**********************Do not forget the sign**********************************/
+		DL_scalar=- gauss->weight*Jdet2d*surface_normal[2]; 
+		/******************************************************************************/
+
+		TripleMultiply( l1l6,1,numdof,1,
+					&DL_scalar,1,1,0,
+					l1l6,1,numdof,0,
+					&Ke_g[0][0],0);
+
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
 	return Ke;
 }
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8415)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8416)
@@ -1087,52 +1087,4 @@
 
 	/*Clean up and return*/
-	return Ke;
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticVertSurface {{{1*/
-ElementMatrix* Tria::CreateKMatrixDiagnosticVertSurface(void){
-
-	/*Constants*/
-	const int numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
-	int       i,j,ig;
-	double    xyz_list[NUMVERTICES][3];
-	double    surface_normal[3];
-	double    Jdet,DL_scalar;
-	double    L[3];
-	GaussTria *gauss=NULL;
-	double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	SurfaceNormal(&surface_normal[0],xyz_list);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for (ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-		/**********************Do not forget the sign**********************************/
-		DL_scalar=- gauss->weight*Jdet*surface_normal[2]; 
-		/******************************************************************************/
-
-		TripleMultiply( L,1,3,1,
-					&DL_scalar,1,1,0,
-					L,1,3,0,
-					&Ke_g[0][0],0);
-
-		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
-	}
-
-	/*Clean up and return*/
-	delete gauss;
 	return Ke;
 }
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8415)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8416)
@@ -153,5 +153,4 @@
 		ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
-		ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
 		ElementMatrix* CreateKMatrixMelting(void);
 		ElementMatrix* CreateKMatrixHydrology(void);
Index: /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp	(revision 8415)
+++ /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp	(revision 8416)
@@ -149,4 +149,14 @@
 		coords4=(double*)xmalloc(numgauss*sizeof(double));
 		for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
+	}
+	/*Upper surface Tria*/
+	else if(index1==3 && index2==4 && index3==5){
+
+		/*Get GaussTria*/
+		GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
+
+		/*compute z coordinate*/
+		coords4=(double*)xmalloc(numgauss*sizeof(double));
+		for(int i=0;i<numgauss;i++) coords4[i]=1.0;
 	}
 	else{