Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5853)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5854)
@@ -687,6 +687,7 @@
 void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
 
-	/*retrive parameters: */
+	/*retrieve parameters: */
 	int analysis_type;
+	ElementMatrix* Ke=NULL;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
@@ -722,5 +723,7 @@
 			break;
 		case MeltingAnalysisEnum:
-			CreateKMatrixMelting( Kgg);
+			Ke=CreateKMatrixMelting();
+			if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+			delete Ke;
 			break;
 		default:
@@ -2321,5 +2324,5 @@
 	  bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
 	  the stiffness matrix. */
-	if (!IsOnBed()) return NULL;
+	if (!IsOnBed() || IsOnWater()) return NULL;
 
 	/*Depth Averaging B*/
@@ -2821,21 +2824,13 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixMelting {{{1*/
-void  Penta::CreateKMatrixMelting(Mat Kgg){
-
-	Tria* tria=NULL;
-
-	/*If on water, skip: */
-	if(IsOnWater())return;
-
-	if (!IsOnBed()){
-		return;
-	}
-	else{
-
-		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->CreateKMatrixMelting(Kgg);
-		delete tria->matice; delete tria;
-		return;
-	}
+ElementMatrix* Penta::CreateKMatrixMelting(void){
+
+	if (!IsOnBed() || IsOnWater()) return NULL;
+
+	Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+	ElementMatrix* Ke=tria->CreateKMatrixMelting();
+
+	delete tria->matice; delete tria;
+	return Ke;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5853)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5854)
@@ -140,5 +140,5 @@
 		ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
 		void	  CreateKMatrixDiagnosticVert( Mat Kgg);
-		void	  CreateKMatrixMelting(Mat Kggg);
+		ElementMatrix* CreateKMatrixMelting(void);
 		void	  CreateKMatrixPrognostic(Mat Kggg);
 		void	  CreateKMatrixSlope(Mat Kggg);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5853)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5854)
@@ -3154,23 +3154,11 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrixMelting {{{1*/
-void  Tria::CreateKMatrixMelting(Mat Kgg){
-
-	/*indexing: */
-	int i,j;
+ElementMatrix* Tria::CreateKMatrixMelting(void){
 
 	const int  numdof=NUMVERTICES*NDOF1;
-	int*       doflist=NULL;
-
-	/*Grid data: */
+	int i,j,ig;
 	double     xyz_list[NUMVERTICES][3];
-
-	/*Material constants */
 	double     heatcapacity,latentheat;
-
-	/* gaussian points: */
-	int     ig;
 	GaussTria *gauss=NULL;
-
-	/*matrices: */
 	double     Jdet;
 	double     D_scalar;
@@ -3180,12 +3168,12 @@
 	double     Ke_gaussian[numdof][numdof]={0.0};
 
-	/*Recover constants of ice */
+	/*Initialize Element matrix and return if necessary*/
+	if(IsOnWater()) return NULL;
+	ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	latentheat=matpar->GetLatentHeat();
 	heatcapacity=matpar->GetHeatCapacity();
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
 
 	/* Start looping on the number of gauss  (nodes on the bedrock) */
@@ -3195,29 +3183,17 @@
 		gauss->GaussPoint(ig);
 
+		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
-
-		/*Get L matrix : */
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-		/*Calculate DL on gauss point */
 		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
 
-		/*  Do the triple product tL*D*L: */
 		MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
 		MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
 
-		for(i=0;i<NUMVERTICES;i++){
-			for(j=0;j<NUMVERTICES;j++){
-				K_terms[i][j]+=Ke_gaussian[i][j];
-			}
-		}
-	}
-
-	/*Add Ke_gg to global matrix Kgg: */
-	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
+	}
 
 	/*Clean up and return*/
 	delete gauss;
-	xfree((void**)&doflist);
+	return Ke;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5853)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5854)
@@ -133,5 +133,5 @@
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		void	  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);
-		void	  CreateKMatrixMelting(Mat Kgg);
+		ElementMatrix* CreateKMatrixMelting(void);
 		ElementMatrix* CreateKMatrixPrognostic(void);
 		ElementMatrix* CreateKMatrixPrognostic_CG(void);