Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6149)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6150)
@@ -830,5 +830,5 @@
 			GetSolutionFromInputsDiagnosticHoriz(solution);
 		}
-		else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){
+		else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
 			return; //the elements around will create the solution
 		}
@@ -4438,4 +4438,7 @@
 		InputUpdateFromSolutionDiagnosticPattynStokes(solution);
 	}
+	else if (approximation==MacAyealStokesApproximationEnum){
+		InputUpdateFromSolutionDiagnosticMacAyealStokes(solution);
+	}
 	else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
 		InputUpdateFromSolutionDiagnosticStokes(solution);
@@ -4627,4 +4630,98 @@
 	xfree((void**)&doflistp);
 	xfree((void**)&doflistm);
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(double* solution){
+
+	int i;
+
+	const int    numdofm=NDOF2*NUMVERTICES;
+	const int    numdofs=NDOF4*NUMVERTICES;
+	const int    numdof2d=NDOF2*NUMVERTICES2D;
+	int*         doflistm=NULL;
+	int*         doflists=NULL;
+	double       macayeal_values[numdofm];
+	double       stokes_values[numdofs];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vzmacayeal[NUMVERTICES];
+	double       vzstokes[NUMVERTICES];
+	double       vel[NUMVERTICES];
+	int          dummy;
+	double       pressure[NUMVERTICES];
+	double       xyz_list[NUMVERTICES][3];
+	double       stokesreconditioning;
+
+	double *vzmacayeal_ptr         = NULL;
+	Penta  *penta          = NULL;
+
+	/*OK, we have to add results of this element for macayeal 
+	 * and results from the penta at base for macayeal. Now recover results*/
+	penta=GetBasalElement();
+
+	/*Get dof listof this element (macayeal dofs) and of the penta at base (macayeal dofs): */
+	penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
+	GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
+	this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+
+	/*Get node data: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numdof2d;i++){
+		macayeal_values[i]=solution[doflistm[i]];
+		macayeal_values[i+numdof2d]=solution[doflistm[i]];
+	}
+	for(i=0;i<numdofs;i++){
+		stokes_values[i]=solution[doflists[i]];
+	}
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<NUMVERTICES;i++){
+		vx[i]=stokes_values[i*NDOF4+0]+macayeal_values[i*NDOF2+0];
+		vy[i]=stokes_values[i*NDOF4+1]+macayeal_values[i*NDOF2+1];
+		vzstokes[i]=stokes_values[i*NDOF4+2];
+		pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
+	}
+
+	/*Get Vz*/
+	Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);
+	if (vzmacayeal_input){
+		if (vzmacayeal_input->Enum()!=PentaVertexInputEnum){
+			ISSMERROR("Cannot compute Vel as VzMacAyeal is of type %s",EnumToString(vzmacayeal_input->Enum()));
+		}
+		vzmacayeal_input->GetValuesPtr(&vzmacayeal_ptr,&dummy);
+		for(i=0;i<NUMVERTICES;i++) vzmacayeal[i]=vzmacayeal_ptr[i];
+	}
+	else{
+		ISSMERROR("Cannot update solution as VzMacAyeal is not present");
+	}
+
+	/*Now Compute vel*/
+	for(i=0;i<NUMVERTICES;i++) {
+		vz[i]=vzmacayeal[i]+vzstokes[i];
+		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	}
+
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	this->inputs->ChangeEnum(VxEnum,VxOldEnum);
+	this->inputs->ChangeEnum(VyEnum,VyOldEnum);
+	this->inputs->ChangeEnum(VzEnum,VzOldEnum);
+	this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+	this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+	this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
+	this->inputs->AddInput(new PentaVertexInput(VzStokesEnum,vzstokes));
+	this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
+	this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+
+	/*Free ressources:*/
+	xfree((void**)&doflistm);
+	xfree((void**)&doflists);
 }
 /*}}}*/
@@ -4889,4 +4986,5 @@
 	double       vy[NUMVERTICES];
 	double       vz[NUMVERTICES];
+	double       vzmacayeal[NUMVERTICES];
 	double       vzpattyn[NUMVERTICES];
 	double       vzstokes[NUMVERTICES];
@@ -4943,5 +5041,5 @@
 	else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
 
-	/*Do some modifications if we actually have a PattynStokes element*/
+	/*Do some modifications if we actually have a PattynStokes or MacAyealStokes element*/
 	if(approximation==PattynStokesApproximationEnum){
 		Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
@@ -4957,4 +5055,17 @@
 		}
 	}
+	else if(approximation==MacAyealStokesApproximationEnum){
+		Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
+		if (vzstokes_input){
+			if (vzstokes_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as VzStokes is of type %s",EnumToString(vy_input->Enum()));
+			vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy);
+			for(i=0;i<NUMVERTICES;i++) vzstokes[i]=vzstokes_ptr[i];
+		}
+		else ISSMERROR("Cannot compute Vz as VzStokes in not present in MacAyealStokes element");
+		for(i=0;i<NUMVERTICES;i++){
+			vzmacayeal[i]=vz[i];
+			vz[i]=vzmacayeal[i]+vzstokes[i];
+		}
+	}
 
 	/*Now Compute vel*/
@@ -4964,5 +5075,5 @@
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
 	 *so the pressure is just the pressure at the z elevation: except it this is a PattynStokes element */
-	if(approximation!=PattynStokesApproximationEnum){
+	if(approximation!=PattynStokesApproximationEnum &&  approximation!=MacAyealStokesApproximationEnum){
 		rho_ice=matpar->GetRhoIce();
 		g=matpar->GetG();
@@ -4975,10 +5086,13 @@
 	this->inputs->ChangeEnum(VzEnum,VzOldEnum);
 
-	if(approximation!=PattynStokesApproximationEnum){
+	if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){
 		this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
 		this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
 	}
-	else{
+	else if(approximation==PattynStokesApproximationEnum){
 		this->inputs->AddInput(new PentaVertexInput(VzPattynEnum,vzpattyn));
+	}
+	else if(approximation==MacAyealStokesApproximationEnum){
+		this->inputs->AddInput(new PentaVertexInput(VzMacAyealEnum,vzpattyn));
 	}
 
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6149)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6150)
@@ -202,4 +202,5 @@
 		void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong);
+		void    InputUpdateFromSolutionDiagnosticMacAyealStokes( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticPattynStokes( double* solutiong);