Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 24129)
+++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 24130)
@@ -351,4 +351,18 @@
 				_error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
 		}
+	}
+	else{
+		_error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
+	}
+
+}
+/*}}}*/
+int   Numericalflux::GetNumberOfNodesOneSide(void){/*{{{*/
+
+	if(this->flux_degree==P0DGEnum){
+		return 1;
+	}
+	else if(this->flux_degree==P1DGEnum){
+		return 2;
 	}
 	else{
@@ -695,15 +709,18 @@
 
 	/* Intermediaries*/
-	IssmDouble DL1,DL2,Jdet,vx,vy,UdotN;
+	IssmDouble A1,A2,Jdet,vx,vy,UdotN;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble normal[2];
 
 	/*Fetch number of nodes for this flux*/
-	int numnodes = this->GetNumberOfNodes();
+	int numnodes       = this->GetNumberOfNodes();
+	int numnodes_plus  = this->GetNumberOfNodesOneSide();
+	int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
+	_assert_(numnodes==numnodes_plus+numnodes_minus);
 
 	/*Initialize variables*/
-	ElementMatrix *Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
-	IssmDouble    *B      = xNew<IssmDouble>(numnodes);
-	IssmDouble    *Bprime = xNew<IssmDouble>(numnodes);
+	ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
+	IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
+	IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
 
 	/*Retrieve all inputs and parameters*/
@@ -722,6 +739,6 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
-		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
+		tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
+		tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -729,18 +746,52 @@
 		UdotN=vx*normal[0]+vy*normal[1];
 		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL1=gauss->weight*Jdet*dt*UdotN/2;
-		DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
-
-		for(int i=0;i<numnodes;i++){
-			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j]+=DL1*B[i]*Bprime[j];
-				Ke->values[i*numnodes+j]+=DL2*B[i]*B[j];
+		A1=gauss->weight*Jdet*dt*UdotN/2;
+		A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
+
+		/*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
+		 *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
+		 *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
+		 *
+		 *Term 2 (stabilization)  |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
+		 *                                      = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
+		 *     | A++ | A+- |
+		 * K = |-----------|
+		 *     | A-+ | A-- |
+		 *
+		 *These 4 terms for each expressions are added independently*/
+
+		/*First term A++*/
+		for(int i=0;i<numnodes_plus;i++){
+			for(int j=0;j<numnodes_plus;j++){
+				Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
+				Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
 			}
 		}
+		/*Second term A+-*/
+		for(int i=0;i<numnodes_plus;i++){
+			for(int j=0;j<numnodes_minus;j++){
+				Ke->values[i*numnodes+numnodes_plus+j] +=  A1*(basis_minus[j]*basis_plus[i]);
+				Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
+			}
+		}
+		/*Third term A-+*/
+		for(int i=0;i<numnodes_minus;i++){
+			for(int j=0;j<numnodes_plus;j++){
+				Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
+				Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
+			}
+		}
+		/*Fourth term A-+*/
+		for(int i=0;i<numnodes_minus;i++){
+			for(int j=0;j<numnodes_minus;j++){
+				Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
+				Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] +=  A2*(basis_minus[j]*basis_minus[i]);
+			}
+		}
 	}
 
 	/*Clean up and return*/
-   xDelete<IssmDouble>(B);
-   xDelete<IssmDouble>(Bprime);
+   xDelete<IssmDouble>(basis_plus);
+   xDelete<IssmDouble>(basis_minus);
 	delete gauss;
 	return Ke;
Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 24129)
+++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 24130)
@@ -64,4 +64,5 @@
 		void GetNodesSidList(int* sidlist);
 		int  GetNumberOfNodes(void);
+		int  GetNumberOfNodesOneSide(void);
 		bool IsPenalty(void);
 		void PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};