Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 23800)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 23801)
@@ -1070,4 +1070,105 @@
 }
 /*}}}*/
+IssmDouble Penta::GetIcefrontArea(){/*{{{*/
+	
+	IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
+	IssmDouble	Haverage,frontarea;
+	IssmDouble  x1,y1,x2,y2,distance;
+	IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
+	int* indices=NULL;
+	IssmDouble* H=NULL;;
+	int nrfrontbed,numiceverts;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&bed[0],BedEnum);
+	GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+	GetInputListOnVertices(&bases[0],BaseEnum);
+	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+
+	if(!this->IsOnBase()) return 0;
+	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
+	nrfrontbed=0;
+	for(int i=0;i<NUMVERTICES2D;i++){
+		/*Find if bed<0*/
+		if(bed[i]<0.) nrfrontbed++;
+	}
+
+	if(nrfrontbed==3){
+		/*2. Find coordinates of where levelset crosses 0*/
+		int         numiceverts;
+		IssmDouble  s[2],x[2],y[2];
+		this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
+		_assert_(numiceverts); 
+
+		/*3 Write coordinates*/
+		IssmDouble  xyz_list[NUMVERTICES][3];
+		::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+		int counter = 0;
+		if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
+			for(int i=0;i<numiceverts;i++){
+				for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
+					x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+					y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+					counter++;
+				}
+			}
+		}
+		else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
+
+			for(int i=0;i<NUMVERTICES2D;i++){
+				if(lsf[indices[i]]==0.){
+					x[counter]=xyz_list[indices[i]][0];
+					y[counter]=xyz_list[indices[i]][1];
+					counter++;
+				}
+				if(counter==2) break;
+			}
+			if(counter==1){
+				/*We actually have only 1 vertex on levelset, write a single point as a segment*/
+				x[counter]=x[0];
+				y[counter]=y[0];
+				counter++;
+			}
+		}
+		else{
+			_error_("not sure what's going on here...");
+		}
+		x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
+		distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
+		
+		int numthk=numiceverts+2;
+		H=xNew<IssmDouble>(numthk);
+		for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
+
+		switch(numiceverts){
+			case 1: // average over triangle
+				H[0]=Haux[0];
+				H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+				H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+				Haverage=(H[1]+H[2])/2;
+				break;
+			case 2: // average over quadrangle
+				H[0]=Haux[0];
+				H[1]=Haux[1];
+				H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+				H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+				Haverage=(H[2]+H[3])/2;
+				break;
+			default:
+				_error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
+				break;
+		}
+		frontarea=distance*Haverage;
+	}
+	else return 0;	
+	
+	xDelete<int>(indices);
+	xDelete<IssmDouble>(H);
+	
+	_assert_(frontarea>0);
+	return frontarea;
+}
+/*}}}*/
 void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 
@@ -1131,4 +1232,80 @@
 
 	return lower_penta;
+}
+/*}}}*/
+void       Penta::GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
+
+	/* GetLevelsetIntersection computes: 
+	 * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
+	 * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
+
+	/*Intermediaries*/
+	int i, numiceverts, numnoiceverts;
+	int ind0, ind1, lastindex;
+	int indices_ice[NUMVERTICES2D],indices_noice[NUMVERTICES2D];
+	IssmDouble lsf[NUMVERTICES];
+	int* indices = xNew<int>(NUMVERTICES2D);
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&lsf[0],levelset_enum);
+
+	/* Determine distribution of ice over element.
+	 * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
+	lastindex=0;
+	for(i=0;i<NUMVERTICES2D;i++){ // go backwards along vertices, and check for sign change
+		ind0=(NUMVERTICES2D-i)%NUMVERTICES2D;
+		ind1=(ind0-1+NUMVERTICES2D)%NUMVERTICES2D;
+		if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
+			if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
+				lastindex=ind1;
+			else
+				lastindex=ind0;
+			break;
+		}
+	}
+
+	numiceverts=0;
+	numnoiceverts=0;
+	for(i=0;i<NUMVERTICES2D;i++){
+		ind0=(lastindex+i)%NUMVERTICES2D;
+		if(lsf[i]<=level){
+			indices_ice[numiceverts]=i;
+			numiceverts++;
+		}
+		else{
+			indices_noice[numnoiceverts]=i;
+			numnoiceverts++;
+		}
+	}
+	//merge indices 
+	for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
+	for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
+
+	switch (numiceverts){
+		case 0: // no vertex has ice: element is ice free, no intersection
+			for(i=0;i<2;i++)
+				fraction[i]=0.;
+			break;
+		case 1: // one vertex has ice:
+			for(i=0;i<2;i++){
+				fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
+			}
+			break;
+		case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
+			for(i=0;i<2;i++){
+				fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
+			}
+			break;
+		case NUMVERTICES2D: // all vertices have ice: return triangle area
+			for(i=0;i<2;i++)
+				fraction[i]=1.;
+			break;
+		default:
+			_error_("Wrong number of ice vertices in Penta::GetLevelsetIntersectionBase!");
+			break;
+	}
+
+	*pindices=indices;
+	*pnumiceverts=numiceverts;
 }
 /*}}}*/
@@ -2310,4 +2487,63 @@
 	if(this->hneighbors) this->hneighbors->reset();
 
+}
+/*}}}*/
+void       Penta::RignotMeltParameterization(){/*{{{*/
+
+   IssmDouble A, B, alpha, beta;
+	IssmDouble bed,qsg,qsg_basin,TF,yts;
+	int numbasins;
+	IssmDouble basinid[NUMVERTICES];
+	IssmDouble* basin_icefront_area=NULL;
+
+	/* Coefficients */
+	A    = 3e-4;        
+	B    = 0.15;        
+	alpha = 0.39;
+	beta = 1.18;
+	
+	/*Get inputs*/
+	Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
+	Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);		 _assert_(qsg_input);
+	Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
+	GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);	
+	
+	this->FindParam(&yts, ConstantsYtsEnum);
+	this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+	this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
+
+	IssmDouble meltrates[NUMVERTICES2D];  //frontal melt-rate
+	
+	/* Start looping on the number of vertices: */
+	GaussPenta* gauss=new GaussPenta();
+	for(int iv=0;iv<NUMVERTICES2D;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get variables */
+		bed_input->GetInputValue(&bed,gauss);
+		qsg_input->GetInputValue(&qsg,gauss);
+		TF_input->GetInputValue(&TF,gauss);
+
+		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 */
+			qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
+
+			/* calculate melt rates */
+			meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
+		}	
+
+		if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
+		if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum));
+	
+	this->InputExtrude(CalvingMeltingrateEnum,-1);
+   
+	/*Cleanup and return*/
+	xDelete<IssmDouble>(basin_icefront_area);
+	delete gauss;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 23800)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 23801)
@@ -75,8 +75,10 @@
 		void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
 		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
+		IssmDouble		GetIcefrontArea();
 		void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void           GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
 		void           GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
+		void				GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
 		Node*          GetNode(int node_number);
 		int            GetNodeIndex(Node* node);
@@ -150,4 +152,5 @@
 		void           ResetFSBasalBoundaryCondition(void);
 		void           ResetHooks();
+		void				RignotMeltParameterization();
 		void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M);
 		void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23800)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23801)
@@ -3420,6 +3420,6 @@
 
 	/* Coefficients */
-	A    = 3e-4;        //
-	B    = 0.15;        //
+	A    = 3e-4;        
+	B    = 0.15;      
 	alpha = 0.39;
 	beta = 1.18;
@@ -3458,5 +3458,4 @@
 		if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
 		if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
-		
 	}
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23800)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23801)
@@ -1498,7 +1498,5 @@
 void FemModel::IcefrontAreax(){/*{{{*/
 
-	int numvertices      = this->GetElementsWidth();
 	int numbasins;
-	IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
 	this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
 	IssmDouble* basin_icefront_area           = xNewZeroInit<IssmDouble>(numbasins);
@@ -1510,6 +1508,9 @@
 		for(int i=0;i<this->elements->Size();i++){
 			Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+			if(!element->IsOnBase()) continue;
+			int  numvertices = element->GetNumberOfVertices();
+			IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
 			element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
-			for(int j=0;j<numvertices;j++){
+			for(int j=0;j<3;j++){
 				if(BasinId[j]==basin){
 					local_icefront_area+=element->GetIcefrontArea();
@@ -1517,15 +1518,15 @@
 				}
 			}
+			xDelete<IssmDouble>(BasinId);
 		}
 		ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
 		ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	
+		
 		basin_icefront_area[basin-1]=total_icefront_area;
 	}
 	
 	this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
+	xDelete<IssmDouble>(basin_icefront_area);
 	
-	xDelete<IssmDouble>(basin_icefront_area);
-	xDelete<IssmDouble>(BasinId);
 }/*}}}*/
 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/