Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24729)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24730)
@@ -632,5 +632,5 @@
 //		_error_("S");
 //	}
-
+	
 	h=element->CharacteristicLength();
 	
@@ -936,6 +936,7 @@
 	/*Intermediaries */
 	IssmDouble Jdet,D_scalar;
-	IssmDouble vx,vy;
+	IssmDouble vx,vy,dvxdx,dvydy;
 	IssmDouble* xyz_list = NULL;
+	IssmDouble dvx[2],dvy[2];
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -964,8 +965,15 @@
 		vxaverage_input->GetInputValue(&vx,gauss);
 		vyaverage_input->GetInputValue(&vy,gauss);
+		vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+		dvxdx=dvx[0];
+		dvydy=dvy[1];
 
 		D_scalar = gauss->weight*Jdet;
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
+				/*\phi_i \phi_j \nabla\cdot v*/
+				Ke->values[i*numnodes+j] += -D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+				/*\phi_i v\cdot\nabla\phi_j*/
 				Ke->values[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
 			}
@@ -1021,4 +1029,101 @@
 	return Me;
 }/*}}}*/
+ElementVector* MasstransportAnalysis::CreateFctPVector(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries */
+	int			stabilization,dim,domaintype;
+	int         melt_style,point1;
+	bool        mainlyfloating;
+	IssmDouble  fraction1,fraction2;
+	IssmDouble  Jdet;
+	IssmDouble  ms,mb,gmb,fmb;
+	IssmDouble  vx,vy,dvxdx,dvydy;
+	IssmDouble  dvx[2],dvy[2];
+	IssmDouble  gllevelset,phi=1.;
+	IssmDouble* xyz_list = NULL;
+	Gauss*      gauss     = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DverticalEnum:   dim = 1; break;
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 2; break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+	  
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
+	element->FindParam(&stabilization,MasstransportStabilizationEnum);
+	Input2* gmb_input        = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
+	Input2* fmb_input        = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
+	Input2* gllevelset_input = element->GetInput2(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
+	Input2* ms_input         = element->GetInput2(SmbMassBalanceEnum);                       _assert_(ms_input);
+	Input2* vxaverage_input  = element->GetInput2(VxAverageEnum);										_assert_(vxaverage_input);
+	Input2* vyaverage_input  = element->GetInput2(VyAverageEnum);										_assert_(vyaverage_input);
+	
+	/*Recover portion of element that is grounded*/
+	phi=element->GetGroundedPortion(xyz_list);
+	if(melt_style==SubelementMelt2Enum){
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+	   gauss = element->NewGauss(point1,fraction1,fraction2,3);
+	}
+	else{
+		gauss = element->NewGauss(3);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		ms_input->GetInputValue(&ms,gauss);
+		gmb_input->GetInputValue(&gmb,gauss);
+		fmb_input->GetInputValue(&fmb,gauss);
+		gllevelset_input->GetInputValue(&gllevelset,gauss);
+
+		if(melt_style==SubelementMelt1Enum){
+			if (phi>0.999999999) mb=gmb;
+			else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+		}
+		else if(melt_style==SubelementMelt2Enum){
+			if(gllevelset>0.) mb=gmb;
+			else mb=fmb;
+		}
+		else if(melt_style==NoMeltOnPartiallyFloatingEnum){
+			if (phi<0.00000001) mb=fmb;
+			else mb=gmb;
+		}
+		else if(melt_style==FullMeltOnPartiallyFloatingEnum){
+			if (phi<0.99999999) mb=fmb;
+			else mb=gmb;
+		}
+		else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb)*basis[i];
+	
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	delete gauss;
+	return pe;
+}/*}}}*/
 void           MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
 
@@ -1049,4 +1154,24 @@
 	}
 }/*}}}*/
+void           MasstransportAnalysis::FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel){/*{{{*/
+
+	/*Output*/
+	Vector<IssmDouble>* pf = NULL;
+
+	/*Initialize P vector*/
+	AllocateSystemMatricesx(NULL,NULL,NULL,&pf,femmodel);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementVector* pe      = this->CreateFctPVector(element);
+		if(pe) pe->AddToGlobal(pf);
+		delete pe;
+	}
+	pf->Assemble();
+
+	/*Assign output pointer*/
+	*ppf=pf;
+}/*}}}*/
 void           MasstransportAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 24729)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 24730)
@@ -38,5 +38,7 @@
 		ElementMatrix* CreateFctKMatrix(Element* element);
 		ElementMatrix* CreateMassMatrix(Element* element);
+		ElementVector* CreateFctPVector(Element* element);
 		void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
+		void				FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel);
 		void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
 };