Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16858)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16859)
@@ -37,6 +37,9 @@
 		void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 		/*L1L2*/
+		ElementVector* CreatePVectorL1L2(Element* element);
+		ElementVector* CreatePVectorL1L2DrivingStress(Element* element);
+		ElementVector* CreatePVectorL1L2Front(Element* element);
 		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
 		/*HO*/
 		ElementMatrix* CreateKMatrixHO(Element* element);
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16858)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16859)
@@ -838,6 +838,8 @@
 			return NULL;
 		case SSAApproximationEnum: 
 			return CreatePVectorSSA(element);
+		case L1L2ApproximationEnum: 
+			return CreatePVectorL1L2(element);
 		case HOApproximationEnum: 
 			return CreatePVectorHO(element);
 		case FSApproximationEnum: 
@@ -1338,6 +1340,144 @@
 }/*}}}*/
 
 /*L1L2*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorL1L2DrivingStress(basalelement);
+	ElementVector* pe2=CreatePVectorL1L2Front(basalelement);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2DrivingStress(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  thickness,Jdet,slope[2];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe    = element->NewElementVector(L1L2ApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 
+	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2Front(Element* element){/*{{{*/
+
+	/*If no front, return NULL*/
+	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
+	IssmDouble  surface_under_water,base_under_water,pressure;
+	IssmDouble *xyz_list = NULL;
+	IssmDouble *xyz_list_front = NULL;
+	IssmDouble  normal[2];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector(L1L2ApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* bed_input       = element->GetInput(BedEnum);       _assert_(bed_input);
+	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+	element->NormalSection(&normal[0],xyz_list_front);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		bed_input->GetInputValue(&bed,gauss);
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
+		base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
+		water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
+		ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
+		pressure = ice_pressure + water_pressure;
+
+		for (int i=0;i<numnodes;i++){
+			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
+			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_front);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+	return NULL;
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,meshtype;