Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24118)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24119)
@@ -329,4 +329,260 @@
 	delete gauss;
 
+}
+/*}}}*/
+void       Penta::CalvingFluxLevelset(){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
+		IssmDouble flux_per_area=0;
+		this->inputs->AddInput(new PentaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+	}
+	else{
+		int               domaintype,index1,index2;
+		const IssmPDouble epsilon = 1.e-15;
+		IssmDouble        s1,s2;
+		IssmDouble        gl[NUMVERTICES];
+		IssmDouble        xyz_front[2][3];
+
+		IssmDouble *xyz_list = NULL;
+		this->GetVerticesCoordinates(&xyz_list);
+
+		/*Recover parameters and values*/
+		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+		/*Be sure that values are not zero*/
+		if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+		if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+		if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+		int pt1 = 0;
+		int pt2 = 1;
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[2]/(gl[2]-gl[1]);
+			s2=gl[2]/(gl[2]-gl[0]);
+			if(gl[2]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+			/*Portion of the segments*/
+			s1=gl[0]/(gl[0]-gl[1]);
+			s2=gl[0]/(gl[0]-gl[2]);
+			if(gl[0]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[1]/(gl[1]-gl[0]);
+			s2=gl[1]/(gl[1]-gl[2]);
+			if(gl[1]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		}
+		else{
+			_error_("case not possible");
+		}
+
+		/*Some checks in debugging mode*/
+		_assert_(s1>=0 && s1<=1.); 
+		_assert_(s2>=0 && s2<=1.); 
+
+		/*Get normal vector*/
+		IssmDouble normal[3];
+		this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+		normal[0] = -normal[0];
+		normal[1] = -normal[1];
+
+		/*Get inputs*/
+		IssmDouble flux = 0.;
+		IssmDouble area = 0.;
+		IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
+		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+		Input* calvingratex_input=NULL;
+		Input* calvingratey_input=NULL;
+		calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+		calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+
+		/*Start looping on Gaussian points*/
+		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			thickness_input->GetInputValue(&thickness,gauss);
+			calvingratex_input->GetInputValue(&calvingratex,gauss);
+			calvingratey_input->GetInputValue(&calvingratey,gauss);
+			this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+			flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
+			area += Jdet*gauss->weight*thickness; 
+
+			flux_per_area=flux/area;
+		}
+
+		this->inputs->AddInput(new PentaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+		
+		/*Clean up and return*/
+		delete gauss;
+	}
+}
+/*}}}*/
+void       Penta::CalvingMeltingFluxLevelset(){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
+		IssmDouble flux_per_area=0;
+		this->inputs->AddInput(new PentaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
+	}
+	else{
+		int               domaintype,index1,index2;
+		const IssmPDouble epsilon = 1.e-15;
+		IssmDouble        s1,s2;
+		IssmDouble        gl[NUMVERTICES];
+		IssmDouble        xyz_front[2][3];
+
+		IssmDouble *xyz_list = NULL;
+		this->GetVerticesCoordinates(&xyz_list);
+
+		/*Recover parameters and values*/
+		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+		/*Be sure that values are not zero*/
+		if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+		if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+		if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+		int pt1 = 0;
+		int pt2 = 1;
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[2]/(gl[2]-gl[1]);
+			s2=gl[2]/(gl[2]-gl[0]);
+			if(gl[2]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+			/*Portion of the segments*/
+			s1=gl[0]/(gl[0]-gl[1]);
+			s2=gl[0]/(gl[0]-gl[2]);
+			if(gl[0]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[1]/(gl[1]-gl[0]);
+			s2=gl[1]/(gl[1]-gl[2]);
+			if(gl[1]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		}
+		else{
+			_error_("case not possible");
+		}
+
+		/*Some checks in debugging mode*/
+		_assert_(s1>=0 && s1<=1.); 
+		_assert_(s2>=0 && s2<=1.); 
+
+		/*Get normal vector*/
+		IssmDouble normal[3];
+		this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+		normal[0] = -normal[0];
+		normal[1] = -normal[1];
+
+		/*Get inputs*/
+		IssmDouble flux = 0.;
+		IssmDouble area = 0.;
+		IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
+		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+		Input* calvingratex_input=NULL;
+		Input* calvingratey_input=NULL;
+		Input* vx_input=NULL;
+		Input* vy_input=NULL;
+		Input* meltingrate_input=NULL;
+		calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+		calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+		meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
+
+		/*Start looping on Gaussian points*/
+		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			thickness_input->GetInputValue(&thickness,gauss);
+			calvingratex_input->GetInputValue(&calvingratex,gauss);
+			calvingratey_input->GetInputValue(&calvingratey,gauss);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			vel=vx*vx+vy*vy;
+			meltingrate_input->GetInputValue(&meltingrate,gauss);	
+			meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
+			meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
+			this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+			flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
+			area += Jdet*gauss->weight*thickness; 
+
+			flux_per_area=flux/area;
+		}
+
+		this->inputs->AddInput(new PentaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
+		
+		/*Clean up and return*/
+		delete gauss;
+	}
 }
 /*}}}*/
@@ -1504,4 +1760,122 @@
 	/*Clean up and return*/
 	return groundedarea;
+}
+/*}}}*/
+IssmDouble Penta::IcefrontMassFluxLevelset(bool scaled){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
+
+	/*Scaled not implemented yet...*/
+	_assert_(!scaled); 
+
+	int               domaintype,index1,index2;
+	const IssmPDouble epsilon = 1.e-15;
+	IssmDouble        s1,s2;
+	IssmDouble        gl[NUMVERTICES];
+	IssmDouble        xyz_front[2][3];
+
+	IssmDouble *xyz_list = NULL;
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+	int pt1 = 0;
+	int pt2 = 1;
+	if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[2]/(gl[2]-gl[1]);
+		s2=gl[2]/(gl[2]-gl[0]);
+		if(gl[2]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+	}
+	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+		/*Portion of the segments*/
+		s1=gl[0]/(gl[0]-gl[1]);
+		s2=gl[0]/(gl[0]-gl[2]);
+		if(gl[0]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+	}
+	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[1]/(gl[1]-gl[0]);
+		s2=gl[1]/(gl[1]-gl[2]);
+		if(gl[1]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+	}
+	else{
+		_error_("case not possible");
+	}
+
+
+	/*Some checks in debugging mode*/
+	_assert_(s1>=0 && s1<=1.); 
+	_assert_(s2>=0 && s2<=1.); 
+
+	/*Get normal vector*/
+	IssmDouble normal[3];
+	this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+	normal[0] = -normal[0];
+	normal[1] = -normal[1];
+	
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble vx,vy,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+	vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
+	}
+
+	return flux;
 }
 /*}}}*/
@@ -2148,4 +2522,13 @@
 }
 /*}}}*/
+Gauss*     Penta::NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz){/*{{{*/
+
+	IssmDouble  area_coordinates[2][3];
+
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+
+	return new GaussPenta(area_coordinates,order_horiz);
+}
+/*}}}*/
 Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
 	return new GaussPenta(vertex1,vertex2,order);
@@ -2297,4 +2680,19 @@
 
 	for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
+}
+/*}}}*/
+void       Penta::NormalSectionBase(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
+
+	/*Build unit outward pointing vector*/
+	IssmDouble vector[2];
+	IssmDouble norm;
+
+	vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+	vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+
+	norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
+
+	normal[0]= + vector[1]/norm;
+	normal[1]= - vector[0]/norm;
 }
 /*}}}*/
@@ -2503,5 +2901,7 @@
 void       Penta::RignotMeltParameterization(){/*{{{*/
 
-   IssmDouble A, B, alpha, beta;
+	if(!this->IsOnBase()) return;
+   
+	IssmDouble A, B, alpha, beta;
 	IssmDouble bed,qsg,qsg_basin,TF,yts;
 	int numbasins;
@@ -2543,5 +2943,5 @@
 
 			/* calculate melt rates */
-			meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
+			meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
 		}	
 
@@ -2685,4 +3085,5 @@
 	if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
 	if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
+	
 	Tria* tria=(Tria*)SpawnTria(0,1,2);
 	switch(this->material->ObjectEnum()){
@@ -3006,4 +3407,255 @@
 	return dt;
 }/*}}}*/
+IssmDouble Penta::TotalCalvingFluxLevelset(bool scaled){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
+
+	/*Scaled not implemented yet...*/
+	_assert_(!scaled); 
+
+	int               domaintype,index1,index2;
+	const IssmPDouble epsilon = 1.e-15;
+	IssmDouble        s1,s2;
+	IssmDouble        gl[NUMVERTICES];
+	IssmDouble        xyz_front[2][3];
+
+	IssmDouble *xyz_list = NULL;
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+	int pt1 = 0;
+	int pt2 = 1;
+	if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[2]/(gl[2]-gl[1]);
+		s2=gl[2]/(gl[2]-gl[0]);
+		if(gl[2]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+	}
+	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+		/*Portion of the segments*/
+		s1=gl[0]/(gl[0]-gl[1]);
+		s2=gl[0]/(gl[0]-gl[2]);
+		if(gl[0]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+	}
+	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[1]/(gl[1]-gl[0]);
+		s2=gl[1]/(gl[1]-gl[2]);
+		if(gl[1]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+	}
+	else{
+		_error_("case not possible");
+	}
+
+
+	/*Some checks in debugging mode*/
+	_assert_(s1>=0 && s1<=1.); 
+	_assert_(s2>=0 && s2<=1.); 
+
+	/*Get normal vector*/
+	IssmDouble normal[3];
+	this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+	normal[0] = -normal[0];
+	normal[1] = -normal[1];
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble calvingratex,calvingratey,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* calvingratex_input=NULL;
+	Input* calvingratey_input=NULL;
+	calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+	calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		calvingratex_input->GetInputValue(&calvingratex,gauss);
+		calvingratey_input->GetInputValue(&calvingratey,gauss);
+		this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
+	}
+
+	return flux;
+		
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
+IssmDouble Penta::TotalCalvingMeltingFluxLevelset(bool scaled){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
+
+	/*Scaled not implemented yet...*/
+	_assert_(!scaled); 
+
+	int               domaintype,index1,index2;
+	const IssmPDouble epsilon = 1.e-15;
+	IssmDouble        s1,s2;
+	IssmDouble        gl[NUMVERTICES];
+	IssmDouble        xyz_front[2][3];
+
+	IssmDouble *xyz_list = NULL;
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+	int pt1 = 0;
+	int pt2 = 1;
+	if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[2]/(gl[2]-gl[1]);
+		s2=gl[2]/(gl[2]-gl[0]);
+		if(gl[2]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+	}
+	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+		/*Portion of the segments*/
+		s1=gl[0]/(gl[0]-gl[1]);
+		s2=gl[0]/(gl[0]-gl[2]);
+		if(gl[0]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+	}
+	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+		/*Portion of the segments*/
+		s1=gl[1]/(gl[1]-gl[0]);
+		s2=gl[1]/(gl[1]-gl[2]);
+		if(gl[1]<0.){
+			pt1 = 1; pt2 = 0;
+		}
+
+		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+	}
+	else{
+		_error_("case not possible");
+	}
+
+
+	/*Some checks in debugging mode*/
+	_assert_(s1>=0 && s1<=1.); 
+	_assert_(s2>=0 && s2<=1.); 
+
+	/*Get normal vector*/
+	IssmDouble normal[3];
+	this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+	normal[0] = -normal[0];
+	normal[1] = -normal[1];
+	
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* calvingratex_input=NULL;
+	Input* calvingratey_input=NULL;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* meltingrate_input=NULL;
+	calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+	calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+	vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+	vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
+	meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		calvingratex_input->GetInputValue(&calvingratex,gauss);
+		calvingratey_input->GetInputValue(&calvingratey,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		meltingrate_input->GetInputValue(&meltingrate,gauss);	
+		meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
+		meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
+		this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
+	}
+
+	return flux;
+	
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
 IssmDouble Penta::TotalFloatingBmb(bool scaled){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24118)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24119)
@@ -50,4 +50,6 @@
 		void           CalvingRateVonmises();
 		void           CalvingRateLevermann();
+		void           CalvingFluxLevelset();
+		void           CalvingMeltingFluxLevelset();
 		IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
 		void           ComputeBasalStress(void);
@@ -93,4 +95,5 @@
 		void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
 		IssmDouble     GroundedArea(bool scaled);
+		IssmDouble		IcefrontMassFluxLevelset(bool scaled);
 		IssmDouble     IceVolume(bool scaled);
 		IssmDouble     IceVolumeAboveFloatation(bool scaled);
@@ -122,4 +125,5 @@
 		Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
 		Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
+		Gauss*         NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz);
 		Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
 		Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
@@ -139,4 +143,5 @@
 		void	         NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
 		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
+		void           NormalSectionBase(IssmDouble* normal,IssmDouble* xyz_list);
 		void	         NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
 		int            NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
@@ -163,4 +168,6 @@
 		int            TensorInterpolation(){_error_("not implemented yet");};
 		IssmDouble     TimeAdapt();
+		IssmDouble		TotalCalvingFluxLevelset(bool scaled);
+		IssmDouble		TotalCalvingMeltingFluxLevelset(bool scaled);
 		IssmDouble     TotalFloatingBmb(bool scaled);
 		IssmDouble     TotalGroundedBmb(bool scaled);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24118)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24119)
@@ -492,127 +492,129 @@
 	}
 	else{
-	int               domaintype,index1,index2;
-	const IssmPDouble epsilon = 1.e-15;
-	IssmDouble        s1,s2;
-	IssmDouble        gl[NUMVERTICES];
-	IssmDouble        xyz_front[2][3];
-
-
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
-
-	/*Recover parameters and values*/
-	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
-
-	/*Be sure that values are not zero*/
-	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
-	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
-	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
-
-	if(domaintype==Domain2DverticalEnum){
-		_error_("not implemented");
-	}
-	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
-		int pt1 = 0;
-		int pt2 = 1;
-		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-
-			/*Portion of the segments*/
-			s1=gl[2]/(gl[2]-gl[1]);
-			s2=gl[2]/(gl[2]-gl[0]);
-			if(gl[2]<0.){
-				pt1 = 1; pt2 = 0;
+		int               domaintype,index1,index2;
+		const IssmPDouble epsilon = 1.e-15;
+		IssmDouble        s1,s2;
+		IssmDouble        gl[NUMVERTICES];
+		IssmDouble        xyz_front[2][3];
+
+		IssmDouble *xyz_list = NULL;
+		this->GetVerticesCoordinates(&xyz_list);
+
+		/*Recover parameters and values*/
+		parameters->FindParam(&domaintype,DomainTypeEnum);
+		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+		/*Be sure that values are not zero*/
+		if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+		if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+		if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+		if(domaintype==Domain2DverticalEnum){
+			_error_("not implemented");
+		}
+		else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
+			int pt1 = 0;
+			int pt2 = 1;
+			if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+				/*Portion of the segments*/
+				s1=gl[2]/(gl[2]-gl[1]);
+				s2=gl[2]/(gl[2]-gl[0]);
+				if(gl[2]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+				xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+				xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+				xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+				xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+				xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+				xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
 			}
-			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
-		}
-		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-
-			/*Portion of the segments*/
-			s1=gl[0]/(gl[0]-gl[1]);
-			s2=gl[0]/(gl[0]-gl[2]);
-			if(gl[0]<0.){
-				pt1 = 1; pt2 = 0;
+			else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+				/*Portion of the segments*/
+				s1=gl[0]/(gl[0]-gl[1]);
+				s2=gl[0]/(gl[0]-gl[2]);
+				if(gl[0]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+
+				xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+				xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+				xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+				xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+				xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+				xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
 			}
-
-			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
-		}
-		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-
-			/*Portion of the segments*/
-			s1=gl[1]/(gl[1]-gl[0]);
-			s2=gl[1]/(gl[1]-gl[2]);
-			if(gl[1]<0.){
-				pt1 = 1; pt2 = 0;
+			else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+				/*Portion of the segments*/
+				s1=gl[1]/(gl[1]-gl[0]);
+				s2=gl[1]/(gl[1]-gl[2]);
+				if(gl[1]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+
+				xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+				xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+				xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+				xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+				xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+				xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
 			}
-
-			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+			else{
+				_error_("case not possible");
+			}
+
+		}
+		else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
+
+		/*Some checks in debugging mode*/
+		_assert_(s1>=0 && s1<=1.); 
+		_assert_(s2>=0 && s2<=1.); 
+
+		/*Get normal vector*/
+		IssmDouble normal[3];
+		this->NormalSection(&normal[0],&xyz_front[0][0]);
+		normal[0] = -normal[0];
+		normal[1] = -normal[1];
+
+		/*Get inputs*/
+		IssmDouble flux = 0.;
+		IssmDouble area = 0.;
+		IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
+		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+		Input* calvingratex_input=NULL;
+		Input* calvingratey_input=NULL;
+		if(domaintype==Domain2DhorizontalEnum){
+			calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+			calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
 		}
 		else{
-			_error_("case not possible");
-		}
-
-	}
-	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
-
-	/*Some checks in debugging mode*/
-	_assert_(s1>=0 && s1<=1.); 
-	_assert_(s2>=0 && s2<=1.); 
-
-	/*Get normal vector*/
-	IssmDouble normal[3];
-	this->NormalSection(&normal[0],&xyz_front[0][0]);
-	normal[0] = -normal[0];
-	normal[1] = -normal[1];
-
-	/*Get inputs*/
-	IssmDouble flux = 0.;
-	IssmDouble area = 0.;
-	IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
-	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
-	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* calvingratex_input=NULL;
-	Input* calvingratey_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-		calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-	}
-	else{
-		calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-		calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-	}
-
-	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-		thickness_input->GetInputValue(&thickness,gauss);
-		calvingratex_input->GetInputValue(&calvingratex,gauss);
-		calvingratey_input->GetInputValue(&calvingratey,gauss);
-		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
-
-		flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
-		area += Jdet*gauss->weight*thickness; 
-
-		flux_per_area=flux/area;
-	}
-	
-	this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+			calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+			calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+		}
+
+		/*Start looping on Gaussian points*/
+		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			thickness_input->GetInputValue(&thickness,gauss);
+			calvingratex_input->GetInputValue(&calvingratex,gauss);
+			calvingratey_input->GetInputValue(&calvingratey,gauss);
+			this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
+
+			flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
+			area += Jdet*gauss->weight*thickness; 
+
+			flux_per_area=flux/area;
+		}
+
+		this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+		
+		/*Clean up and return*/
+		delete gauss;
 	}
 }
@@ -626,139 +628,142 @@
 	}
 	else{
-	int               domaintype,index1,index2;
-	const IssmPDouble epsilon = 1.e-15;
-	IssmDouble        s1,s2;
-	IssmDouble        gl[NUMVERTICES];
-	IssmDouble        xyz_front[2][3];
-
-
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
-
-	/*Recover parameters and values*/
-	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
-
-	/*Be sure that values are not zero*/
-	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
-	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
-	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
-
-	if(domaintype==Domain2DverticalEnum){
-		_error_("not implemented");
-	}
-	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
-		int pt1 = 0;
-		int pt2 = 1;
-		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-
-			/*Portion of the segments*/
-			s1=gl[2]/(gl[2]-gl[1]);
-			s2=gl[2]/(gl[2]-gl[0]);
-			if(gl[2]<0.){
-				pt1 = 1; pt2 = 0;
+		int               domaintype,index1,index2;
+		const IssmPDouble epsilon = 1.e-15;
+		IssmDouble        s1,s2;
+		IssmDouble        gl[NUMVERTICES];
+		IssmDouble        xyz_front[2][3];
+
+
+		IssmDouble *xyz_list = NULL;
+		this->GetVerticesCoordinates(&xyz_list);
+
+		/*Recover parameters and values*/
+		parameters->FindParam(&domaintype,DomainTypeEnum);
+		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+		/*Be sure that values are not zero*/
+		if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+		if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+		if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+		if(domaintype==Domain2DverticalEnum){
+			_error_("not implemented");
+		}
+		else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
+			int pt1 = 0;
+			int pt2 = 1;
+			if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+				/*Portion of the segments*/
+				s1=gl[2]/(gl[2]-gl[1]);
+				s2=gl[2]/(gl[2]-gl[0]);
+				if(gl[2]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+				xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+				xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+				xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+				xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+				xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+				xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
 			}
-			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
-		}
-		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-
-			/*Portion of the segments*/
-			s1=gl[0]/(gl[0]-gl[1]);
-			s2=gl[0]/(gl[0]-gl[2]);
-			if(gl[0]<0.){
-				pt1 = 1; pt2 = 0;
+			else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+				/*Portion of the segments*/
+				s1=gl[0]/(gl[0]-gl[1]);
+				s2=gl[0]/(gl[0]-gl[2]);
+				if(gl[0]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+
+				xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+				xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+				xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+				xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+				xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+				xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
 			}
-
-			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
-		}
-		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-
-			/*Portion of the segments*/
-			s1=gl[1]/(gl[1]-gl[0]);
-			s2=gl[1]/(gl[1]-gl[2]);
-			if(gl[1]<0.){
-				pt1 = 1; pt2 = 0;
+			else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+				/*Portion of the segments*/
+				s1=gl[1]/(gl[1]-gl[0]);
+				s2=gl[1]/(gl[1]-gl[2]);
+				if(gl[1]<0.){
+					pt1 = 1; pt2 = 0;
+				}
+
+				xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+				xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+				xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+				xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+				xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+				xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
 			}
-
-			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+			else{
+				_error_("case not possible");
+			}
+
+		}
+		else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
+
+		/*Some checks in debugging mode*/
+		_assert_(s1>=0 && s1<=1.); 
+		_assert_(s2>=0 && s2<=1.); 
+
+		/*Get normal vector*/
+		IssmDouble normal[3];
+		this->NormalSection(&normal[0],&xyz_front[0][0]);
+		normal[0] = -normal[0];
+		normal[1] = -normal[1];
+
+		/*Get inputs*/
+		IssmDouble flux = 0.;
+		IssmDouble area = 0.;
+		IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
+		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+		Input* calvingratex_input=NULL;
+		Input* calvingratey_input=NULL;
+		Input* vx_input=NULL;
+		Input* vy_input=NULL;
+		Input* meltingrate_input=NULL;
+		if(domaintype==Domain2DhorizontalEnum){
+			calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+			calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+			vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+			vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+			meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
 		}
 		else{
-			_error_("case not possible");
-		}
-
-	}
-	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
-
-	/*Some checks in debugging mode*/
-	_assert_(s1>=0 && s1<=1.); 
-	_assert_(s2>=0 && s2<=1.); 
-
-	/*Get normal vector*/
-	IssmDouble normal[3];
-	this->NormalSection(&normal[0],&xyz_front[0][0]);
-	normal[0] = -normal[0];
-	normal[1] = -normal[1];
-
-	/*Get inputs*/
-	IssmDouble flux = 0.;
-	IssmDouble area = 0.;
-	IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
-	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
-	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* calvingratex_input=NULL;
-	Input* calvingratey_input=NULL;
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	Input* meltingrate_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-		calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-		meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
-	}
-	else{
-		calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-		calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-	}
-
-	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-		thickness_input->GetInputValue(&thickness,gauss);
-		calvingratex_input->GetInputValue(&calvingratex,gauss);
-		calvingratey_input->GetInputValue(&calvingratey,gauss);
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vel=vx*vx+vy*vy;
-		meltingrate_input->GetInputValue(&meltingrate,gauss);	
-		meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
-		meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
-		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
-
-		flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
-		area += Jdet*gauss->weight*thickness; 
-
-		flux_per_area=flux/area;
-	}
-	
-	this->inputs->AddInput(new TriaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
+			calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+			calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+		}
+
+		/*Start looping on Gaussian points*/
+		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			thickness_input->GetInputValue(&thickness,gauss);
+			calvingratex_input->GetInputValue(&calvingratex,gauss);
+			calvingratey_input->GetInputValue(&calvingratey,gauss);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			vel=vx*vx+vy*vy;
+			meltingrate_input->GetInputValue(&meltingrate,gauss);	
+			meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
+			meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
+			this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
+
+			flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
+			area += Jdet*gauss->weight*thickness; 
+
+			flux_per_area=flux/area;
+		}
+
+		this->inputs->AddInput(new TriaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
+		
+		/*Clean up and return*/
+		delete gauss;
 	}
 }
@@ -2299,5 +2304,5 @@
 		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
 	}
-
+	delete gauss;
 	return flux;
 }
@@ -2428,5 +2433,5 @@
 		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
 	}
-
+	delete gauss;
 	return flux;
 }
@@ -3673,5 +3678,5 @@
 		TF_input->GetInputValue(&TF,gauss);
 
-		if(basinid[iv]==0 || basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
+		if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
 		else{
 			/* change the unit of qsg (m^3/d -> m/d) with ice front area */
@@ -3679,5 +3684,5 @@
 
 			/* calculate melt rates */
-			meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
+			meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
 		}	
 
@@ -4257,5 +4262,5 @@
 	/*Get inputs*/
 	IssmDouble flux = 0.;
-	IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
+	IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
 	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
 	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 24118)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 24119)
@@ -433,4 +433,35 @@
 }
 /*}}}*/
+GaussPenta::GaussPenta(IssmDouble area_coordinates[2][3],int order_horiz){/*{{{*/
+
+	/*Intermediaties*/
+	IssmPDouble *seg_horiz_coords  = NULL;
+	IssmPDouble *seg_horiz_weights = NULL;
+
+	/*get the gauss points using the product of two line rules*/
+	GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
+
+	/*Allocate GaussPenta fields*/
+	numgauss=order_horiz;
+	coords1=xNew<IssmDouble>(numgauss);
+	coords2=xNew<IssmDouble>(numgauss);
+	coords3=xNew<IssmDouble>(numgauss);
+	coords4=xNew<IssmDouble>(numgauss);
+	weights=xNew<IssmDouble>(numgauss);
+
+	/*Quads: get the gauss points using the product of two line rules  */
+	for(int i=0;i<order_horiz;i++){
+		coords1[i]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
+		coords2[i]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
+		coords3[i]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
+		coords4[i]=0.;
+		weights[i]=seg_horiz_weights[i];
+	}
+
+	/*clean-up*/
+	xDelete<IssmPDouble>(seg_horiz_coords);
+	xDelete<IssmPDouble>(seg_horiz_weights);
+}
+/*}}}*/
 GaussPenta::~GaussPenta(){/*{{{*/
 	xDelete<IssmDouble>(weights);
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h	(revision 24118)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h	(revision 24119)
@@ -37,4 +37,5 @@
 		GaussPenta(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
 		GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert);
+		GaussPenta(IssmDouble area_coordinates[2][3],int order_horiz);
 		~GaussPenta();