Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5802)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5803)
@@ -252,4 +252,6 @@
 	VzObsEnum,
 	VzOldEnum,
+	VzPattynEnum,
+	VzStokesEnum,
 	QmuVzEnum,
 	WeightsEnum,
Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5802)
+++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5803)
@@ -226,4 +226,6 @@
 		case VzObsEnum : return "VzObs";
 		case VzOldEnum : return "VzOld";
+		case VzPattynEnum : return "VzPattyn";
+		case VzStokesEnum : return "VzStokes";
 		case QmuVzEnum : return "QmuVz";
 		case WeightsEnum : return "Weights";
Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5802)
+++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5803)
@@ -224,4 +224,6 @@
 	else if (strcmp(name,"VzObs")==0) return VzObsEnum;
 	else if (strcmp(name,"VzOld")==0) return VzOldEnum;
+	else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum;
+	else if (strcmp(name,"VzStokes")==0) return VzStokesEnum;
 	else if (strcmp(name,"QmuVz")==0) return QmuVzEnum;
 	else if (strcmp(name,"Weights")==0) return WeightsEnum;
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 5802)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 5803)
@@ -31,4 +31,5 @@
 	/*Fetch data: */
 	IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
+	IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes");
 
 	/*Initialize counter*/
@@ -41,5 +42,9 @@
 		if(iomodel->my_vertices[i]){
 
-			if ((int)iomodel->spcvelocity[6*i+2]){
+			if ((int)iomodel->gridonstokes[i]){
+				constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,DiagnosticVertAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for Stokes
+				count++;
+			}
+			else if ((int)iomodel->spcvelocity[6*i+2]){
 				constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,
 								*(iomodel->spcvelocity+6*i+5)/iomodel->yts,DiagnosticVertAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
@@ -52,4 +57,5 @@
 	/*Free data: */
 	xfree((void**)&iomodel->spcvelocity);
+	xfree((void**)&iomodel->gridonstokes);
 
 	cleanup_and_return:
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5802)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5803)
@@ -367,12 +367,5 @@
 	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
 	if (analysis_type==DiagnosticHorizAnalysisEnum){
-		int approximation;
-		inputs->GetParameterValue(&approximation,ApproximationEnum);
-		if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
-			InputUpdateFromSolutionDiagnosticStokes( solution);
-		}
-		else{
-			InputUpdateFromSolutionDiagnosticHoriz( solution);
-		}
+		InputUpdateFromSolutionDiagnosticHoriz( solution);
 	}
 	else if (analysis_type==DiagnosticHutterAnalysisEnum){
@@ -901,6 +894,9 @@
 			GetSolutionFromInputsDiagnosticStokes(solution);
 		}
-		else{
+		else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
 			GetSolutionFromInputsDiagnosticHoriz(solution);
+		}
+		else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){
+			return; //the elements around will create the solution
 		}
 	}
@@ -3603,7 +3599,7 @@
 	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
-	Input* vz_input=NULL;
+	Input* vzstokes_input=NULL;
 	if(approximation==PattynStokesApproximationEnum){
-		vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
+		vzstokes_input=inputs->GetInput(VzStokesEnum); ISSMASSERT(vzstokes_input);
 	}
 
@@ -3618,5 +3614,5 @@
 		vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
 		if(approximation==PattynStokesApproximationEnum){
-			vz_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+			vzstokes_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 			dwdz=dw[2];
 		}
@@ -4009,10 +4005,5 @@
 	/*If the element is a coupling, do nothing: every grid is also on an other elements 
 	 * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
-	if(approximation==MacAyealPattynApproximationEnum){
-	 return;
-	}
-	else{
-		GetDofList(&doflist,approximation,GsetEnum);
-	}
+	GetDofList(&doflist,approximation,GsetEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -4460,4 +4451,10 @@
 	else if (approximation==PattynApproximationEnum){
 		InputUpdateFromSolutionDiagnosticPattyn(solution);
+	}
+	else if (approximation==PattynStokesApproximationEnum){
+		InputUpdateFromSolutionDiagnosticPattynStokes(solution);
+	}
+	else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
+		InputUpdateFromSolutionDiagnosticStokes(solution);
 	}
 	else if (approximation==MacAyealPattynApproximationEnum){
@@ -4726,4 +4723,98 @@
 
 /*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){
+
+	int i;
+
+	const int    numdofpervertexp=2;
+	const int    numdofpervertexs=4;
+	const int    numdofp=numdofpervertexp*NUMVERTICES;
+	const int    numdofs=numdofpervertexs*NUMVERTICES;
+	int*         doflistp=NULL;
+	int*         doflists=NULL;
+	double       pattyn_values[numdofp];
+	double       stokes_values[numdofs];
+	double       vx[NUMVERTICES];
+	double       vy[NUMVERTICES];
+	double       vz[NUMVERTICES];
+	double       vzpattyn[NUMVERTICES];
+	double       vzstokes[NUMVERTICES];
+	double       vel[NUMVERTICES];
+	int          dummy;
+	double       pressure[NUMVERTICES];
+	double       surface[NUMVERTICES];
+	double       rho_ice,g;
+	double       xyz_list[NUMVERTICES][3];
+
+	double *vzpattyn_ptr         = NULL;
+	Penta  *penta          = NULL;
+
+	/*OK, we have to add results of this element for pattyn 
+	 * and results from the penta at base for macayeal. Now recover results*/
+	penta=GetBasalElement();
+
+	/*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
+	GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
+	GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
+
+	/*Get node data: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numdofp;i++){
+		pattyn_values[i]=solution[doflistp[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*numdofpervertexs+0]+pattyn_values[i*numdofpervertexp+0];
+		vy[i]=stokes_values[i*numdofpervertexs+1]+pattyn_values[i*numdofpervertexp+1];
+		vzstokes[i]=stokes_values[i*numdofpervertexs+2];
+		pressure[i]=stokes_values[i*numdofpervertexs+3];
+	}
+
+	/*Get Vz*/
+	Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);
+	if (vzpattyn_input){
+		if (vzpattyn_input->Enum()!=PentaVertexInputEnum){
+			ISSMERROR("Cannot compute Vel as VzPattyn is of type %s",EnumToString(vzpattyn_input->Enum()));
+		}
+		vzpattyn_input->GetValuesPtr(&vzpattyn_ptr,&dummy);
+		for(i=0;i<NUMVERTICES;i++) vzpattyn[i]=vzpattyn_ptr[i];
+	}
+	else{
+		ISSMERROR("Cannot update solution as VzPattyn is not present");
+	}
+
+	/*Now Compute vel*/
+	for(i=0;i<NUMVERTICES;i++) {
+		vz[i]=vzpattyn[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**)&doflistp);
+	xfree((void**)&doflists);
+}
+/*}}}*/
 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{1*/
 void  Penta::InputUpdateFromSolutionDiagnosticHutter(double* solution){
@@ -4815,4 +4906,6 @@
 	double       vy[NUMVERTICES];
 	double       vz[NUMVERTICES];
+	double       vzpattyn[NUMVERTICES];
+	double       vzstokes[NUMVERTICES];
 	double       vel[NUMVERTICES];
 	int          dummy;
@@ -4824,4 +4917,7 @@
 	double*      vx_ptr=NULL;
 	double*      vy_ptr=NULL;
+	double*      vzstokes_ptr=NULL;
+
+	int         approximation;
 
 	/*Get dof list: */
@@ -4830,4 +4926,5 @@
 	/*Get node data: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
 
 	/*Use the dof list to index into the solution vector: */
@@ -4858,23 +4955,46 @@
 	else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
 
+	/*Do some modifications if we actually have a PattynStokes element*/
+	if(approximation==PattynStokesApproximationEnum){
+		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 PattynStokes element");
+		for(i=0;i<NUMVERTICES;i++){
+			vzpattyn[i]=vz[i];
+			vz[i]=vzpattyn[i]+vzstokes[i];
+		}
+	}
+
 	/*Now Compute vel*/
+
 	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.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: */
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-	GetParameterListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	 *so the pressure is just the pressure at the z elevation: except it this is a PattynStokes element */
+	if(!approximation==PattynStokesApproximationEnum){
+		rho_ice=matpar->GetRhoIce();
+		g=matpar->GetG();
+		GetParameterListOnVertices(&surface[0],SurfaceEnum);
+		for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	}
 
 	/*Now, we have to move the previous Vz inputs to old 
-	 * status, otherwise, we'll wipe them off: */
+	 * status, otherwise, we'll wipe them off and add the new inputs: */
 	this->inputs->ChangeEnum(VzEnum,VzOldEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
-
-	/*Add vz and vel as inputs to the penta element: */
+
+	if(!approximation==PattynStokesApproximationEnum){
+		this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+		this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+	}
+	else{
+		this->inputs->AddInput(new PentaVertexInput(VzPattynEnum,vzpattyn));
+	}
+
 	this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
 	this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
-	this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
 
 	/*Free ressources:*/
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5802)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5803)
@@ -4152,7 +4152,7 @@
 	Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
-	Input* vz_input=NULL;
+	Input* vzstokes_input=NULL;
 	if(approximation==PattynStokesApproximationEnum){
-		vz_input=inputs->GetInput(VzEnum);            ISSMASSERT(vz_input);
+		vzstokes_input=inputs->GetInput(VzStokesEnum);       ISSMASSERT(vzstokes_input);
 	}
 	Input* bed_input=inputs->GetInput(BedEnum);             ISSMASSERT(bed_input);
@@ -4172,5 +4172,5 @@
 		vy_input->GetParameterValue(&vy, gauss);
 		if(approximation==PattynStokesApproximationEnum){
-			vz_input->GetParameterValue(&vz, gauss);
+			vzstokes_input->GetParameterValue(&vz, gauss);
 		}
 		else vz=0;