Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 17510)
+++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 17511)
@@ -50,5 +50,5 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,D;
+	IssmDouble  Jdet,D,dt=1.e+9;
 	IssmDouble *xyz_list = NULL;
 
@@ -59,4 +59,5 @@
 		case Mesh2DverticalEnum: dim = 2; break;
 		case Mesh3DEnum:         dim = 3; break;
+		case Mesh3DtetrasEnum:    dim = 3; break;
 		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
 	}
@@ -68,5 +69,4 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -79,11 +79,19 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(Bprime,gauss);
 		GetB(B,element,dim,xyz_list,gauss);
-		D=gauss->weight*Jdet;
+		D=gauss->weight*Jdet*dt;
 
+		/*vertical diffusion*/
 		TripleMultiply(B,1,numnodes,1,
 					&D,1,1,0,
-					Bprime,1,numnodes,0,
+					B,1,numnodes,0,
+					&Ke->values[0],1);
+
+		/*Next value*/
+		D=gauss->weight*Jdet;
+		element->NodalFunctions(B,gauss);
+		TripleMultiply(B,numnodes,1,0,
+					&D,1,1,0,
+					B,1,numnodes,0,
 					&Ke->values[0],1);
 	} 
@@ -92,10 +100,46 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	return Ke;
 }/*}}}*/
 ElementVector* DepthAverageAnalysis::CreatePVector(Element* element){/*{{{*/
-	return NULL;
+
+	/*Intermediaries*/
+	int         input_enum;
+	IssmDouble  Jdet,scalar,value;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe     = element->NewElementVector();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&input_enum,InputToDepthaverageEnum);
+	Input* input = element->GetInput(input_enum); _assert_(input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		/* Build transient now */
+		input->GetInputValue(&value, gauss);
+		scalar=value*Jdet*gauss->weight;
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return pe;
 }/*}}}*/
 void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17510)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17511)
@@ -445,4 +445,31 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::InputUpdateFromSolutionOneDof{{{*/
+void  Tetra::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
+
+	/*Intermediary*/
+	int* doflist = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(int i=0;i<numnodes;i++){
+		values[i]=solution[doflist[i]];
+		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Add input to the element: */
+	this->inputs->AddInput(new TetraInput(enum_type,values,P1Enum));
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
+	xDelete<int>(doflist);
+}
+/*}}}*/
 /*FUNCTION Tetra::IsIceInElement    THIS ONE{{{*/
 bool   Tetra::IsIceInElement(){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17510)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17511)
@@ -47,5 +47,5 @@
 		/*Update virtual functions resolution: {{{*/
 		void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
-		void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
+		void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum);
 		void  InputUpdateFromVector(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
 #ifdef _HAVE_DAKOTA_
Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 17510)
+++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 17511)
@@ -26,5 +26,5 @@
 
 		case StressbalanceSolutionEnum:
-			numanalyses=5;
+			numanalyses=6;
 			analyses=xNew<int>(numanalyses);
 			analyses[0]=StressbalanceAnalysisEnum;
@@ -33,4 +33,5 @@
 			analyses[3]=L2ProjectionBaseAnalysisEnum;
 			analyses[4]=ExtrudeFromBaseAnalysisEnum;
+			analyses[5]=DepthAverageAnalysisEnum;
 			break;
 
Index: /issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- /issm/trunk-jpl/src/c/cores/cores.h	(revision 17510)
+++ /issm/trunk-jpl/src/c/cores/cores.h	(revision 17511)
@@ -30,4 +30,5 @@
 void controltao_core(FemModel* femmodel);
 void masstransport_core(FemModel* femmodel);
+void depthaverage_core(FemModel* femmodel);
 void extrudefrombase_core(FemModel* femmodel);
 void extrudefromtop_core(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp	(revision 17510)
+++ /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp	(revision 17511)
@@ -12,5 +12,5 @@
 void depthaverage_core(FemModel* femmodel){
 
-	if(VerboseSolution()) _printf0_("extruding solution from base...\n");
+	if(VerboseSolution()) _printf0_("depth averaging...\n");
 
 	/*Call on core computations: */
Index: /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp	(revision 17510)
+++ /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp	(revision 17511)
@@ -33,5 +33,5 @@
 	}
 	if(meshtype==Mesh2DverticalEnum){
-	      femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
+		femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
 		extrudefrombase_core(femmodel);
 	}