Index: /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 127)
+++ /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 128)
@@ -24,7 +24,4 @@
 		return DiagnosticVertAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"diagnostic_basevert")==0){
-		return DiagnosticBaseVertAnalysisEnum();
-	}
 	else if (strcmp(analysis_type,"prognostic")==0){
 		return PrognosticAnalysisEnum();
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 127)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 128)
@@ -20,5 +20,4 @@
 int DiagnosticHorizAnalysisEnum(void){  return          23; }
 int DiagnosticVertAnalysisEnum(void){   return          24; }
-int DiagnosticBaseVertAnalysisEnum(void){ return        25; }
 int PrognosticAnalysisEnum(void){       return          26; }
 int ThermalSteadyAnalysisEnum(void)	{   return          27; }
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 127)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 128)
@@ -33,5 +33,4 @@
 int DiagnosticHorizAnalysisEnum(void);
 int DiagnosticVertAnalysisEnum(void);
-int DiagnosticBaseVertAnalysisEnum(void);
 int PrognosticAnalysisEnum(void);
 int ThermalSteadyAnalysisEnum(void);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 127)
+++ /issm/trunk/src/c/Makefile.am	(revision 128)
@@ -163,4 +163,7 @@
 					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
 					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
+					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
 					./ModelProcessorx/Control/CreateParametersControl.cpp\
 					./Dofx/Dofx.h\
@@ -378,4 +381,7 @@
 					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
 					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
+					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
 					./ModelProcessorx/Control/CreateParametersControl.cpp\
 					./Dofx/Dofx.h\
Index: /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp	(revision 127)
+++ /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp	(revision 128)
@@ -21,10 +21,7 @@
 		CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle);
 	}
-	/*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
-		CreateConstraintsDiagnosticBaseVert(pconstraints,model,model_handle);
-	}
 	else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
 		CreateConstraintsDiagnosticVert(pconstraints,model,model_handle);
-	}
+	}/*
 	else if (strcmp(model->analysis_type,"melting")==0){
 		CreateConstraintsMelting(pconstraints,model,model_handle);
Index: /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp	(revision 127)
+++ /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp	(revision 128)
@@ -24,13 +24,10 @@
 	
 	}
-	/*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
+	else if ((strcmp(model->analysis_type,"diagnostic_vert")==0)){
 
-		CreateElementsNodesAndMaterialsDataSetsDiagnosticBaseVert(pelements,pnodes,pmaterials, model,model->analysis_type);
-
+		CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle);
+	
 	}
-	else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
-
-		CreateElementsNodesAndMaterialsDataSetsDiagnosticVert(pelements,pnodes,pmaterials, model,model->analysis_type);
-	}
+	/*
 	else if (strcmp(model->analysis_type,"melting")==0){
 
@@ -46,5 +43,5 @@
 	}*/
 	else{
-		throw ErrorException(__FUNCT__," analysis_type not supported yet!");
+		throw ErrorException(__FUNCT__,exprintf("%s%s%s"," analysis_type: ",model->analysis_type," not supported yet!"));
 	}
 }
Index: /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp	(revision 127)
+++ /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp	(revision 128)
@@ -22,12 +22,8 @@
 		CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
 	}
-	/*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
-
-		CreateLoadsDiagnosticBaseVert(ploads,model,model_handle);
-	}
 	else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
 
 		CreateLoadsDiagnosticVert(ploads,model,model_handle);
-	}
+	}/*
 	else if (strcmp(model->analysis_type,"melting")==0){
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 127)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 128)
@@ -59,4 +59,6 @@
 	double tria_b[3];
 	double tria_k[3];
+	double tria_melting[3];
+	double tria_accumulation[3];
 	int    tria_friction_type;
 	double tria_p;
@@ -278,5 +280,5 @@
 
 			/*Create tria element using its constructor:*/
-			tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
+			tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
 
 			/*Add tria element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 128)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 128)
@@ -0,0 +1,26 @@
+/*
+ * CreateConstraintsDiagnosticHoriz.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateConstraintsDiagnosticHoriz"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../Model.h"
+
+
+void	CreateConstraintsDiagnosticVert(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
+
+
+	DataSet* constraints = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum());
+
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 128)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 128)
@@ -0,0 +1,433 @@
+/*
+ * CreateElementsNodesAndMaterialsDiagnosticVert.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateElementsNodesAndMaterialsDiagnosticVert"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../Model.h"
+
+
+void	CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
+
+
+	/*output: int* epart, int* my_grids, double* my_bordergrids*/
+
+
+	int i,j,k,n;
+	extern int my_rank;
+	extern int num_procs;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    materials = NULL;
+	
+	/*Objects: */
+	Node*       node   = NULL;
+	Penta*      penta = NULL;
+	Matice*     matice  = NULL;
+	Matpar*     matpar  = NULL;
+
+	int         analysis_type;
+	
+	/*output: */
+	int* epart=NULL; //element partitioning.
+	int* npart=NULL; //node partitioning.
+	int* my_grids=NULL;
+	double* my_bordergrids=NULL;
+
+
+	/*intermediary: */
+	const int elements_width=6;
+	double B_avg;
+			
+	/*matice constructor input: */
+	int    matice_mid;
+	double matice_B;
+	double matice_n;
+	
+	/*penta constructor input: */
+
+	int penta_id;
+	int penta_mid;
+	int penta_mparid;
+	int penta_g[6];
+	double penta_h[6];
+	double penta_s[6];
+	double penta_b[6];
+	double penta_k[6];
+	int penta_friction_type;
+	double penta_p;
+	double penta_q;
+	int penta_shelf;
+	int penta_onbed;
+	int penta_onsurface;
+	double penta_meanvel;/*!scaling ratio for velocities*/
+	double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
+	int penta_collapse;
+	double penta_melting[6];
+	double penta_accumulation[6];
+	double penta_geothermalflux[6];
+	int penta_artdiff;
+	int penta_thermal_steadystate;
+	double penta_viscosity_overshoot;
+
+	/*matpar constructor input: */
+	int	matpar_mid;
+	double matpar_rho_ice; 
+	double matpar_rho_water;
+	double matpar_heatcapacity;
+	double matpar_thermalconductivity;
+	double matpar_latentheat;
+	double matpar_beta;
+	double matpar_meltingpoint;
+	double matpar_mixed_layer_capacity;
+	double matpar_thermal_exchange_velocity;
+	double matpar_g;
+
+	/* node constructor input: */
+	int node_id;
+	int node_partitionborder=0;
+	double node_x[3];
+	int node_onbed;
+	int node_onsurface;
+	int node_upper_node_id;
+	int node_numdofs;
+
+		
+	#ifdef _PARALLEL_
+	/*Metis partitioning: */
+	int  range;
+	Vec  gridborder;
+	int  my_numgrids;
+	int* all_numgrids=NULL;
+	int  gridcount;
+	int  count;
+	#endif
+	int  first_grid_index;
+
+	/*Rifts:*/
+	int* riftsnumpenaltypairs;
+	double** riftspenaltypairs;
+	int* riftsfill;
+	double* riftsfriction;
+	double* riftpenaltypairs=NULL;
+	int el1,el2;
+	 
+	/*Penalty partitioning: */
+	int num_grids3d_collapsed;
+	double* double_penalties_grids3d_collapsed=NULL;
+	double* double_penalties_grids3d_noncollapsed=NULL;
+	int     grid_ids[6];
+	int     num_grid_ids;
+	int     grid_id;
+	
+
+	/*Get analysis_type: */
+	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum());
+	nodes     = new DataSet(NodesEnum());
+	materials = new DataSet(MaterialsEnum());
+
+	/*Width of elements: */
+	if(strcmp(model->meshtype,"2d")==0){
+		throw ErrorException(__FUNCT__," 2d mesh not supported for vertical velocity");
+	}
+
+
+	#ifdef _PARALLEL_
+	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
+	if(strcmp(model->meshtype,"2d")==0){
+		/*load elements: */
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+	}
+	else{
+		/*load elements2d: */
+		ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
+	}
+
+	MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
+
+	/*Free elements and elements2d: */
+	xfree((void**)&model->elements);
+	xfree((void**)&model->elements2d);
+		
+
+	/*Deal with rifts, they have to be included into one partition only, not several: */
+	FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
+	
+	for(i=0;i<model->numrifts;i++){
+		riftpenaltypairs=model->riftspenaltypairs[i];
+		for(j=0;j<model->riftsnumpenaltypairs[i];j++){
+			el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
+			el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
+			epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
+		}
+	}
+	/*Free rifts: */
+	xfree((void**)&riftsnumpenaltypairs); 
+	for(i=0;i<model->numrifts;i++){
+		double* temp=riftspenaltypairs[i];
+		xfree((void**)&temp);
+	}
+	xfree((void**)&riftspenaltypairs);
+	xfree((void**)&riftsfill);
+	xfree((void**)&riftsfriction);
+
+	/*Used later on: */
+	my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
+	#endif
+
+
+	/*Create 3d elements: */
+	/*Fetch data needed: */
+	ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
+	ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
+	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
+	ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
+	ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
+	ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
+	ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
+	ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
+	
+	for (i=0;i<model->numberofelements;i++){
+	#ifdef _PARALLEL_
+	/*We are using our element partition to decide which elements will be created on this node: */
+	if(my_rank==epart[i]){
+	#endif
+
+		
+		/*name and id: */
+		penta_id=i+1; //matlab indexing.
+		penta_mid=i+1; //refers to the corresponding material property card
+		penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
+
+		/*vertices,thickness,surface,bed and drag: */
+		for(j=0;j<6;j++){
+			penta_g[j]=(int)*(model->elements+elements_width*i+j);
+			penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1)); 
+			penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1)); 
+			penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1)); 
+			penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1))/model->yts; 
+			penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1))/model->yts; 
+		}
+
+		/*diverse: */
+		penta_shelf=(int)*(model->elementoniceshelf+i);
+		penta_onbed=(int)*(model->elementonbed+i);
+		penta_onsurface=(int)*(model->elementonsurface+i);
+		penta_collapse=1;
+
+		/*Create Penta using its constructor:*/
+		penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
+				penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
+				penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
+				penta_thermal_steadystate,penta_viscosity_overshoot); 
+
+		/*Add penta element to elements dataset: */
+		elements->AddObject(penta);
+
+
+		/*Deal with material:*/
+		matice_mid=i+1; //same as the material id from the geom2 elements.
+
+		/*Create matice using its constructor:*/
+		matice= new Matice(matice_mid,matice_B,matice_n);
+
+		/*Add matice element to materials dataset: */
+		materials->AddObject(matice);
+		
+		#ifdef _PARALLEL_
+		/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
+		 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+		 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
+		 will hold which grids belong to this partition*/
+		my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
+		my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
+		my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
+		my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
+		my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
+		my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
+		#endif
+
+	#ifdef _PARALLEL_
+	}//if(my_rank==epart[i])
+	#endif
+
+	}//for (i=0;i<numberofelements;i++)
+
+	/*Free data: */
+	xfree((void**)&model->elements);
+	xfree((void**)&model->thickness);
+	xfree((void**)&model->surface);
+	xfree((void**)&model->bed);
+	xfree((void**)&model->elementoniceshelf);
+	xfree((void**)&model->elementonbed);
+	xfree((void**)&model->elementonsurface);
+	xfree((void**)&model->elements_type);
+	xfree((void**)&model->melting);
+	xfree((void**)&model->accumulation);
+
+
+	/*Add one constant material property to materials: */
+	matpar_mid=model->numberofelements+1; //put it at the end of the materials
+	matpar_g=model->g; 
+	matpar_rho_ice=model->rho_ice; 
+	matpar_rho_water=model->rho_water; 
+	matpar_thermalconductivity=model->thermalconductivity; 
+	matpar_heatcapacity=model->heatcapacity; 
+	matpar_latentheat=model->latentheat; 
+	matpar_beta=model->beta; 
+	matpar_meltingpoint=model->meltingpoint; 
+	matpar_mixed_layer_capacity=model->mixed_layer_capacity; 
+	matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; 
+
+	/*Create matpar object using its constructor: */
+	matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
+			matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
+			matpar_thermal_exchange_velocity,matpar_g);
+		
+	/*Add to materials datset: */
+	materials->AddObject(matpar);
+	
+	#ifdef _PARALLEL_
+		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
+		 *element partition, and which are on its border with other nodes:*/
+		gridborder=NewVec(model->numberofnodes);
+
+		for (i=0;i<model->numberofnodes;i++){
+			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
+		}
+		VecAssemblyBegin(gridborder);
+		VecAssemblyEnd(gridborder);
+
+		#ifdef _DEBUG_
+		VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
+		#endif
+		
+		VecToMPISerial(&my_bordergrids,gridborder);
+
+		#ifdef _DEBUG_
+		if(my_rank==0){
+			for (i=0;i<model->numberofnodes;i++){
+				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
+			}
+		}
+		#endif
+	#endif
+
+
+	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
+	 We can therefore determine  which grids are internal to this node's partition 
+	 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
+	 that, go and create the grids*/
+
+	/*Create nodes from x,y,z, as well as the spc values on those grids: */
+		
+	/*First fetch data: */
+	if (strcmp(model->meshtype,"3d")==0){
+		ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
+		ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
+	}
+	ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
+	ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
+	ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+
+	
+	/*Get number of dofs per node: */
+	DistributeNumDofs(&node_numdofs,analysis_type);
+
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((my_grids[i]==1)){
+	#endif
+
+		node_id=i+1; //matlab indexing
+			
+		
+		
+		#ifdef _PARALLEL_
+		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+			node_partitionborder=1;
+		}
+		else{
+			node_partitionborder=0;
+		}
+		#else
+			node_partitionborder=0;
+		#endif
+
+
+		node_x[0]=model->x[i];
+		node_x[1]=model->y[i];
+		node_x[2]=model->z[i];
+
+		
+		node_onbed=(int)model->gridonbed[i];
+		node_onsurface=(int)model->gridonsurface[i];	
+		if (isnan(model->uppernodes[i])){
+			node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+		}
+		else{
+			node_upper_node_id=(int)model->uppernodes[i];
+		}
+
+		/*Create node using its constructor: */
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+
+		/*Add node to nodes dataset: */
+		nodes->AddObject(node);
+
+	#ifdef _PARALLEL_
+	} //if((my_grids[i]==1))
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	elements->Presort();
+	nodes->Presort();
+	materials->Presort();
+
+	/*Clean fetched data: */
+	xfree((void**)&model->deadgrids);
+	xfree((void**)&model->x);
+	xfree((void**)&model->y);
+	xfree((void**)&model->z);
+	xfree((void**)&model->gridonbed);
+	xfree((void**)&model->gridonsurface);
+	xfree((void**)&model->uppernodes);
+		
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pmaterials=materials;
+
+	/*Keep partitioning information into model*/
+	model->epart=epart;
+	model->my_grids=my_grids;
+	model->my_bordergrids=my_bordergrids;
+
+	/*Free ressources:*/
+	#ifdef _PARALLEL_
+	xfree((void**)&all_numgrids);
+	VecFree(&gridborder);
+	#endif
+
+}
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp	(revision 128)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp	(revision 128)
@@ -0,0 +1,28 @@
+/*! \file CreateLoadsDiagnosticVert.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateLoadsDiagnosticVert"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../Model.h"
+
+
+void	CreateLoadsDiagnosticVert(DataSet** ploads, Model* model,ConstDataHandle model_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum());
+	
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
+
+
Index: /issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 127)
+++ /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 128)
@@ -170,8 +170,13 @@
 	/*Create of fem datasets: specialised drivers: */
 	
-	/*diagnostic:*/
+	/*diagnostic horizontal*/
 	void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
 	void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle);
+
+	/*diagnostic vertical*/
+	void	CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
+	void	CreateConstraintsDiagnosticVert(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
+	void    CreateLoadsDiagnosticVert(DataSet** ploads, Model* model, ConstDataHandle model_handle);
 	
 	/*control:*/
Index: /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp	(revision 127)
+++ /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp	(revision 128)
@@ -32,5 +32,4 @@
 	loads->Configure(elements, loads, nodes, materials);
 
-
 	/*Get size of matrix: */
 	gsize=nodes->NumberOfDofs();
@@ -39,5 +38,5 @@
 	if(kflag)Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode);
 	if(pflag)pg=NewVec(gsize);
-	
+
 	/*Fill stiffness matrix and right hand side vector, from elements: */
 	if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type);
@@ -59,5 +58,4 @@
 	}
 	
-
 	/*Assign output pointers: */
 	*pKgg=Kgg;
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 127)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 128)
@@ -287,4 +287,9 @@
 
 	}
+	else if ((analysis_type==DiagnosticVertAnalysisEnum())){
+
+		CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type);
+
+	}
 	else{
 		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
@@ -546,4 +551,144 @@
 
 #undef __FUNCT__ 
+#define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
+void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type){
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=6;
+	const int    NDOF1=1;
+	const int    numdofs=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+
+	/* 3d gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* fourth_gauss_vert_coord  =  NULL;
+	double* area_gauss_weights           =  NULL;
+	double* vert_gauss_weights           =  NULL;
+	int     ig1,ig2;
+	double  gauss_weight1,gauss_weight2;
+	double  gauss_l1l2l3l4[4];
+	int     order_area_gauss;
+	int     num_vert_gauss;
+	int     num_area_gauss;
+	double  gauss_weight;
+
+	/* matrices: */
+	double  Ke_gg[numdofs][numdofs];
+	double  Ke_gg_gaussian[numdofs][numdofs];
+	double  Jdet;
+	double  B[NDOF1][numgrids];
+	double  Bprime[NDOF1][numgrids];
+	double  DL_scalar;
+	
+	/*Collapsed formulation: */
+	Tria*  tria=NULL;
+
+	/*If this element is on the bed rock, spawn a Tria element out of the bed triangle, and use it 
+	 * to create the stifness matrix for the base velocity: */
+	if(onbed){
+		tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
+		tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
+		delete tria;
+	}
+
+	/*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
+	 * matrix: */
+	if(onsurface){
+		tria=SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
+		tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type);
+		delete tria;
+	}
+
+	/*Now, onto the formulation for the vertical velocity: */
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set Ke_gg to 0: */
+	for(i=0;i<numdofs;i++){
+		for(j=0;j<numdofs;j++){
+			Ke_gg[i][j]=0.0;
+		}
+	}
+
+	/*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
+	  get tria gaussian points as well as segment gaussian points. For tria gaussian 
+	  points, order of integration is 2, because we need to integrate the product tB*D*B' 
+	  which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
+	  points, same deal, which yields 3 gaussian points.*/
+
+	order_area_gauss=2;
+	num_vert_gauss=2;
+
+	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
+	#ifdef _DEBUG_ 
+	for (i=0;i<num_area_gauss;i++){
+		_printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
+	}
+	for (i=0;i<num_vert_gauss;i++){
+		_printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
+	}
+	#endif
+	
+	/* Start  looping on the number of gaussian points: */
+	for (ig1=0; ig1<num_area_gauss; ig1++){
+		for (ig2=0; ig2<num_vert_gauss; ig2++){
+		
+			/*Pick up the gaussian point: */
+			gauss_weight1=*(area_gauss_weights+ig1);
+			gauss_weight2=*(vert_gauss_weights+ig2);
+			gauss_weight=gauss_weight1*gauss_weight2;
+			
+			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
+			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
+			gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
+			gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+
+			/*Get B and Bprime matrices: */
+			GetB_vert(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
+			GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
+
+			/* Get Jacobian determinant: */
+			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+
+			DL_scalar=gauss_weight*Jdet;
+
+			/*  Do the triple product tB*D*Bprime: */
+			TripleMultiply( &B[0][0],1,numgrids,1,
+					&DL_scalar,1,1,0,
+					&Bprime[0][0],1,numgrids,0,
+					&Ke_gg_gaussian[0][0],0);
+
+			/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+			for( i=0; i<numdofs; i++){
+				for(j=0;j<numdofs;j++){
+					Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+				}
+			}	
+		} //for (ig2=0; ig2<num_vert_gauss; ig2++)
+	} //for (ig1=0; ig1<num_area_gauss; ig1++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+	
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&fourth_gauss_vert_coord);
+	xfree((void**)&area_gauss_weights);
+	xfree((void**)&vert_gauss_weights);
+}
+
+#undef __FUNCT__ 
 #define __FUNCT__ "Penta::CreatePVector"
 void  Penta::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
@@ -553,5 +698,8 @@
 
 		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
-
+	}
+	else if ((analysis_type==DiagnosticVertAnalysisEnum())){
+
+		CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
 	}
 	else{
@@ -741,4 +889,6 @@
 	double   tria_c[3];
 	double   tria_k[3];
+	double   tria_melting[3];
+	double   tria_accumulation[3];
 	
 	/*configuration: */
@@ -763,4 +913,12 @@
 	tria_k[2]=k[g2];
 
+	tria_melting[0]=melting[g0];
+	tria_melting[1]=melting[g1];
+	tria_melting[2]=melting[g2];
+	
+	tria_accumulation[0]=accumulation[g0];
+	tria_accumulation[1]=accumulation[g1];
+	tria_accumulation[2]=accumulation[g2];
+
 	tria_node_ids[0]=node_ids[g0];
 	tria_node_ids[1]=node_ids[g1];
@@ -775,5 +933,5 @@
 	tria_node_offsets[2]=node_offsets[g2];
 
-	tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
+	tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
 
 	tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets);
@@ -1418,2 +1576,190 @@
 
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta:GetB_vert"
+void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+
+
+	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
+	where hi is the interpolation function for grid i.*/
+
+	int i;
+	const int NDOF3=3;
+	const int numgrids=6;
+	double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
+
+	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
+	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
+
+	#ifdef _DEBUG_ 
+	for (i=0;i<numgrids;i++){
+		_printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
+	}
+	#endif
+
+	/*Build B: */
+	for (i=0;i<numgrids;i++){
+		B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i];  
+	}
+	
+}
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta:GetBPrime_vert"
+void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+
+	// Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i
+
+	int i;
+				
+	GetNodalFunctions(B, gauss_l1l2l3l4);
+
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
+void  Penta::CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type){
+
+	int i;
+
+
+	/* node data: */
+	const int    numgrids=6;
+	const int    NDOF1=1;
+	const int    numdofs=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+	
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* fourth_gauss_vert_coord  =  NULL;
+	double* area_gauss_weights      =  NULL;
+	double* vert_gauss_weights      =  NULL;
+	double  gauss_l1l2l3l4[4];
+	int     order_area_gauss;
+	int     num_vert_gauss;
+	int     num_area_gauss;
+	int     ig1,ig2;
+	double  gauss_weight1,gauss_weight2;
+	double  gauss_weight;
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*element vector at the gaussian points: */
+	double  pe_g[numdofs];
+	double  pe_g_gaussian[numdofs];
+	double l1l2l3l4l5l6[6];
+
+	/*Spawning: */
+	Tria* tria=NULL;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double* velocity_param=NULL;
+	double vx_list[numgrids];
+	double vy_list[numgrids];
+	double du[3];
+	double dv[3];
+	double dudx,dvdy;
+
+	/*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 
+	 *diagnostic base vertical stifness: */
+	tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
+	tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
+	delete tria;
+
+	/*Now, handle the standard penta element: */
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set pe_g to 0: */
+	for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+
+	/* recover input parameters: */
+	velocity_param=ParameterInputsRecover(inputs,"velocity");
+
+	if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
+
+	/*Initialize vx_list and vy_list: */
+	for(i=0;i<numgrids;i++){
+		/*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity 
+		 *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
+		vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0]; 
+		vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
+	}
+
+	/*Get gaussian points and weights :*/
+	order_area_gauss=2;
+	num_vert_gauss=2;
+
+	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
+	#ifdef _DEBUG_ 
+	for (i=0;i<num_area_gauss;i++){
+		_printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
+	}
+	for (i=0;i<num_vert_gauss;i++){
+		_printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
+	}
+	#endif
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig1=0; ig1<num_area_gauss; ig1++){
+		for (ig2=0; ig2<num_vert_gauss; ig2++){
+		
+			/*Pick up the gaussian point: */
+			gauss_weight1=*(area_gauss_weights+ig1);
+			gauss_weight2=*(vert_gauss_weights+ig2);
+			gauss_weight=gauss_weight1*gauss_weight2;
+			
+			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
+			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
+			gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
+			gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+	
+			/*Get velocity derivative, with respect to x and y: */
+			GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
+			GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
+			dudx=du[0];
+			dvdy=dv[1];
+			
+
+			/* Get Jacobian determinant: */
+			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+			#ifdef _DEBUG_ 
+			_printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet);
+			#endif
+		
+			 /*Get nodal functions: */
+			GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
+
+			/*Build pe_g_gaussian vector: */
+			for (i=0;i<numgrids;i++){
+				pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l2l3l4l5l6[i];
+			}
+
+			/*Add pe_g_gaussian vector to pe_g: */
+			for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
+
+		} //for (ig2=0; ig2<num_vert_gauss; ig2++)
+	} //for (ig1=0; ig1<num_area_gauss; ig1++)
+
+	/*Add pe_g to global vector pg: */
+	VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&fourth_gauss_vert_coord);
+	xfree((void**)&area_gauss_weights);
+	xfree((void**)&vert_gauss_weights);
+
+
+}
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 127)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 128)
@@ -37,4 +37,6 @@
 		double b[6];
 		double k[6];
+		double melting[6];
+		double accumulation[6];
 		int    friction_type;
 		double p;
@@ -46,6 +48,4 @@
 		double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
 		int    collapse;
-		double melting[6];
-		double accumulation[6];
 		double geothermalflux[6];
 		int    artdiff;
@@ -73,4 +73,5 @@
 		void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
 		void  CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type);
+		void  CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type);
 		void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
 		void  UpdateFromInputs(ParameterInputs* inputs);
@@ -94,4 +95,6 @@
 		void  GetB(double* pB, double* xyz_list, double* gauss_l1l2l3l4);
 		void  GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4);
+		void  GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
+		void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
 		void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4);
 		void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4);
@@ -100,4 +103,5 @@
 		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
 		void  CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type);
+		void  CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type);
 		void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
 		void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 127)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 128)
@@ -32,6 +32,7 @@
 }
 
-Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],
-				int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,double tria_viscosity_overshoot){
+Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3],
+				double tria_accumulation[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
+				double tria_viscosity_overshoot){
 	
 	int i;
@@ -48,4 +49,6 @@
 		b[i]=tria_b[i];
 		k[i]=tria_k[i];
+		melting[i]=tria_melting[i];
+		accumulation[i]=tria_accumulation[i];
 	}
 	matice=NULL;
@@ -82,4 +85,6 @@
 	printf("   b=[%g,%g,%g]\n",b[0],b[1],b[2]);
 	printf("   k=[%g,%g,%g]\n",k[0],k[1],k[2]);
+	printf("   melting=[%g,%g,%g]\n",melting[0],melting[1],melting[2]);
+	printf("   accumulation=[%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2]);
 	printf("   friction_type: %i\n",friction_type);
 	printf("   p: %g\n",p);
@@ -129,4 +134,6 @@
 	memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b);
 	memcpy(marshalled_dataset,&k,sizeof(k));marshalled_dataset+=sizeof(k);
+	memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting);
+	memcpy(marshalled_dataset,&accumulation,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
 	memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
 	memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
@@ -157,4 +164,6 @@
 		+sizeof(b)
 		+sizeof(k)
+		+sizeof(melting)
+		+sizeof(accumulation)
 		+sizeof(friction_type)
 		+sizeof(onbed)
@@ -197,4 +206,6 @@
 	memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
 	memcpy(&k,marshalled_dataset,sizeof(k));marshalled_dataset+=sizeof(k);
+	memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting);
+	memcpy(&accumulation,marshalled_dataset,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
 	memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
 	memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
@@ -2042,2 +2053,327 @@
 }
 
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
+void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
+
+	int i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdofs=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+	
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+
+	/* matrices: */
+	double L[1][3];
+	double DL_scalar;
+
+	double Ke_gg[3][3]; //stiffness matrix 
+	double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
+	double Jdet;
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set Ke_gg to 0: */
+	for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+	
+		//Get L matrix if viscous basal drag present:
+		GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
+
+		DL_scalar=gauss_weight*Jdet;
+		
+		/*  Do the triple producte tL*D*L: */
+		TripleMultiply( &L[0][0],1,3,1,
+					&DL_scalar,1,1,0,
+					&L[0][0],1,3,0,
+					&Ke_gg_gaussian[0][0],0);
+
+		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+	} //for (ig=0; ig<num_gauss; ig++
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+		
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
+void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
+
+	int i,j;
+	
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdofs=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	
+	/* surface normal: */
+	double x4,y4,z4;
+	double x5,y5,z5;
+	double x6,y6,z6;
+	double v46[3];
+	double v56[3];
+	double normal[3];
+	double norm_normal;
+	double nz;
+
+	/*Matrices: */
+	double DL_scalar;
+	double L[3];
+	double Jdet;
+
+	/* local element matrices: */
+	double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
+	double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
+	
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set Ke_gg to 0: */
+	for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/*Build normal vector to the surface:*/
+	
+	x4=xyz_list[0][0];
+	y4=xyz_list[0][1];
+	z4=xyz_list[0][2];
+
+	x5=xyz_list[1][0];
+	y5=xyz_list[1][1];
+	z5=xyz_list[1][2];
+
+	x6=xyz_list[2][0];
+	y6=xyz_list[2][1];
+	z6=xyz_list[2][2];
+
+	v46[0]=x4-x6;
+	v46[1]=y4-y6;
+	v46[2]=z4-z6;
+
+	v56[0]=x5-x6;
+	v56[1]=y5-y6;
+	v56[2]=z5-z6;
+
+	normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6);
+	normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6);
+	normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6);
+
+	norm_normal=sqrt(pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2));
+	nz=1.0/norm_normal*normal[2];
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+	
+		//Get L matrix if viscous basal drag present:
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
+
+		/**********************Do not forget the sign**********************************/
+		DL_scalar=- gauss_weight*Jdet*nz; 
+		/******************************************************************************/
+	
+		/*  Do the triple producte tL*D*L: */
+		TripleMultiply( L,1,3,1,
+					&DL_scalar,1,1,0,
+					L,1,3,0,
+					&Ke_gg_gaussian[0][0],0);
+
+		/* Add the Ke_gg_gaussian, onto Ke_gg: */
+		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+
+		
+	} //for (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+	
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+}
+		
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
+void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type){
+
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdofs=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+	
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*nodal functions: */
+	double l1l2l3[3];
+
+	/*element vector at the gaussian points: */
+	double  pe_g[numdofs];
+	double  pe_g_gaussian[numdofs];
+
+	/* matrices: */
+	double L[numgrids];
+
+	/*input parameters for structural analysis (diagnostic): */
+	double* velocity_param=NULL;
+	double  vx_list[numgrids];
+	double  vy_list[numgrids];
+	double  vx,vy;
+	double  meltingvalue;
+	double  slope[2];
+	double  dbdx,dbdy;
+
+	/* recover input parameters: */
+	velocity_param=ParameterInputsRecover(inputs,"velocity");
+	if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set pe_g to 0: */
+	for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+
+	/*Initialize vx_list and vy_list: */
+	for(i=0;i<numgrids;i++){
+		if(velocity_param){
+			/*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity 
+			 *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
+			vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0]; 
+			vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
+		}
+		else{
+			vx_list[i]=0;
+			vy_list[i]=0;
+		}
+	}
+		
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/*For icesheets: */
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/*Get melting at gaussian point: */
+		GetParameterValue(&meltingvalue, &melting[0],gauss_l1l2l3);
+
+		/*Get velocity at gaussian point: */
+		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
+		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+
+		/*Get bed slope: */
+		GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
+		dbdx=slope[0];
+		dbdy=slope[1];
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+		
+		//Get L matrix if viscous basal drag present:
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
+
+		
+		/*Build gaussian vector: */
+		for(i=0;i<numgrids;i++){
+			pe_g_gaussian[i]=Jdet*gauss_weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
+		}
+	
+		/*Add pe_g_gaussian vector to pe_g: */
+		for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
+	
+	}
+
+	/*Add pe_g to global vector pg: */
+	VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
+
+	cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 127)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 128)
@@ -36,4 +36,6 @@
 		double b[3];
 		double k[3];
+		double melting[3];
+		double accumulation[3];
 		int    friction_type;
 		double p;
@@ -48,5 +50,5 @@
 
 		Tria();
-		Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],
+		Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],double melting[3],double accumulation[3],
 				int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot);
 		~Tria();
@@ -68,4 +70,6 @@
 		void  CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type);
 		void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type);
+		void  CreateKMatrixDiagnosticBaseVert(Mat Kgg,ParameterInputs* inputs,int analysis_type);
+		void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type);
 		void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
 		void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
@@ -87,4 +91,5 @@
 
 		void  CreatePVectorDiagnosticHoriz(Vec pg,ParameterInputs* inputs,int analysis_type);
+		void  CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type);
 		Matpar* GetMatPar();
 		int   GetShelf();
Index: /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp	(revision 127)
+++ /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp	(revision 128)
@@ -19,4 +19,7 @@
 		numdofs=2;
 	}
+	else if (analysis_type==DiagnosticVertAnalysisEnum()){
+		numdofs=1;
+	}
 	else if (analysis_type==ControlAnalysisEnum()){
 		numdofs=2;
Index: sm/trunk/src/m/solutions/cielo/cielo/cielodiagnostic_core_linear.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/cielo/cielodiagnostic_core_linear.m	(revision 127)
+++ 	(revision )
@@ -1,48 +1,0 @@
-function [ru_g m]=cielodiagnostic_core_linear(m,params,inputs,analysis)
-%INPUT ru_g =cielodiagnostic_core_linear(m,params,inputs,analysis)
-
-	%check inputs
-	if isempty(inputs),
-		clear inputs;
-	end
-
-	%recover flag_y_s from m.y_s
-	if isempty(m.y_s),
-		flag_y_s= 0;
-	else
-		flag_y_s= 1;
-	end
-
-	%stiffness and load generation only:
-	params.kflag=1; params.pflag=1;
-	params.ktflag=0;
-
-	%system matrices
-	[rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis);
-	
-	%Reduce tangent matrix from g size to f size
-	[rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); 
-
-	%Cleanup memory
-	rK_gg=IMdb('drop rK_gg');
-	if ~isempty(rdK_gg), rdK_gg=IMdb('drop rdK_gg');end;
-	
-	%Reduce load from g size to f size
-	[rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s );
-	
-	%Cleanup memory
-	[rp_g rK_fs]=IMdb('drop rp_g rK_fs');
-	
-	%Solve	
-	[ru_f]=Solver(rK_ff,rp_f,{},params.solverstring);
-
-	%Cleanup memory
-	[rK_ff rp_f]=IMdb('drop rK_ff rp_f');
-	
-	%Merge back to g set
-	ru_g= Mergesolvec( ru_f, m.G_mn, m.y_s ); 
-	
-	%Cleanup memory
-	ru_f=IMdb('drop ru_f ');
-
-end %end function
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 127)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 128)
@@ -4,18 +4,10 @@
 	t1=clock;
 
-	analysis_dh='diagnostic_horiz';
-	analysis_dbv='diagnostic_basevert';
-	analysis_dv='diagnostic_vert';
-
 	%Build all models requested for diagnostic simulation
 	%remark, partitions are all identical.
-	md.analysis_type=analysis_dh; 
-	m_dh=CreateFemModel(md);
+	md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
 
 	if strcmpi(md.type,'3d'),
-		md.analysis_type=analysis_dbv; 
-		m_dbv=CreateFemModel(md);
-		md.analysis_type=analysis_dv; 
-		m_dv=CreateFemModel(md);
+		md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
 	end
 
@@ -31,32 +23,14 @@
 	
 		%extrude velocities for collapsed penta elements
-		u_g=VelocityExtrude(md,u_g);
+		disp(sprintf('\n%s',['extruding horizontal velocities...']));
+		u_g=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
 
-		%Compute depth averaged velocity and add it to inputs
-		velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);
+		disp(sprintf('\n%s',['computing vertical velocities...']));
+		u_g_vert=diagnostic_core_linear(m_dv,struct('velocity',u_g));
 
-		%swith u_g and velocity_average to cluster partitioning
-		u_g=SwitchPartitioning(u_g,'cluster',part,tpart,[1 2]); 
-		velocity_average=SwitchPartitioning(velocity_average,'cluster',part,tpart,[1 2]); 
-		
-		disp(sprintf('\n%s',['computing basal vertical velocities...']));
-		SetUset(m_dbv);
-		u_g_basevert=diagnostic_core_linear(m_dbv,m_dbv.parameters,struct('velocity',u_g,'velocity_average',velocity_average),analysis_dbv); 
-		
-		%use u_g_basevert to constrain vertical velocity
-		SetUset(m_dv);
-		m_dv.y_g=u_g_basevert;
-
-		%reduce constraining vector to s-set (the one we solve on)
-		m_dv.y_s = Reducevectorg( m_dv.y_g);
-
-		%run core linear to solve for vertical velocity
-		disp(sprintf('\n%s',['computing vertical velocities...']));
-		u_g_vert=diagnostic_core_linear(m_dv,m_dv.parameters,struct('velocity',u_g),analysis_dv); 
-		
-		%load results onto model: carefull, u_g and u_g_vert are in cluster partitioning
-		md.vx=u_g(indx)*md.yts;
-		md.vy=u_g(indy)*md.yts;
-		md.vz=u_g_vert(indz)*md.yts;
+		%load results onto model: 
+		md.vx=u_g(1:2:end)*md.yts;
+		md.vy=u_g(2:2:end)*md.yts;
+		md.vz=u_g_vert*md.yts;
 		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
 		
@@ -65,12 +39,8 @@
 	else
 		
-		%Build partitioning vectors to recover solution
-		indx=m_dh.part(1:2:end);
-		indy=m_dh.part(2:2:end);
-
 		%load results onto model
-		md.vx=u_g(indx)*md.yts;
-		md.vy=u_g(indy)*md.yts;
-		md.vz=zeros(md.numberofgrids,1);
+		md.vx=u_g(1:2:end)*md.yts;
+		md.vy=u_g(2:2:end)*md.yts;
+		md.vz=zeros(m_dv.nodesets.gsize,1);
 		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
 	
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_linear.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core_linear.m	(revision 128)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_linear.m	(revision 128)
@@ -0,0 +1,25 @@
+function u_g=diagnostic_core_linear(m,inputs)
+%DIAGNOSTIC_CORE_LINEAR - linear solution sequence
+%
+%   Usage:
+%      u_g=diagnostic_core_linear(m,inputs)
+
+	%stiffness and load generation only:
+	m.parameters.kflag=1; m.parameters.pflag=1;
+
+	%system matrices
+	[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs);
+
+	%Reduce tangent matrix from g size to f size
+	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
+	
+	%Reduce load from g size to f size
+	[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
+
+	%Solve	
+	[u_f]=Solver(K_ff,p_f,[],m.parameters);
+
+	%Merge back to g set
+	[u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); 
+	
+end %end function
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 127)
+++ /issm/trunk/src/mex/Makefile.am	(revision 128)
@@ -38,5 +38,6 @@
 				TriMeshProcessRifts\
 				TriMeshRefine\
-				UpdateFromInputs
+				UpdateFromInputs\
+				VelocityExtrude
 
 endif
@@ -166,2 +167,5 @@
 			  UpdateFromInputs/UpdateFromInputs.h
 
+VelocityExtrude_SOURCES = VelocityExtrude/VelocityExtrude.cpp\
+			  VelocityExtrude/VelocityExtrude.h
+
Index: /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.cpp
===================================================================
--- /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.cpp	(revision 127)
+++ /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.cpp	(revision 128)
@@ -35,5 +35,4 @@
 	/*write output : */
 	WriteData(UGOUT,ug,0,0,"Vector",NULL);
-
 	/*Free ressources: */
 	delete elements;
@@ -42,7 +41,8 @@
 	delete materials;
 	VecFree(&ug);
-
+	
 	/*end module: */
 	MODULEEND();
+
 }
 
Index: /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.h
===================================================================
--- /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.h	(revision 127)
+++ /issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.h	(revision 128)
@@ -18,6 +18,6 @@
 /* serial input macros: */
 #define ELEMENTS (mxArray*)prhs[0]
-#define LOADS (mxArray*)prhs[1]
-#define NODES (mxArray*)prhs[2]
+#define NODES (mxArray*)prhs[1]
+#define LOADS (mxArray*)prhs[2]
 #define MATERIALS (mxArray*)prhs[3]
 #define UG (mxArray*)prhs[4]
Index: /issm/trunk/todo
===================================================================
--- /issm/trunk/todo	(revision 127)
+++ /issm/trunk/todo	(revision 128)
@@ -1,4 +1,5 @@
 Create sorting routine for ids.
-Update all inputs. Look at that in more details, ie with respect to control methods. When do we do update of inputs?
-Handle input update for Nodes and Materials. 
+Update all inputs. Look at that in more details, ie with respect to control methods. When do we do update of inputs?  Handle input update for Nodes and Materials. 
 Define default cluster in /etc/cielo.rc
+VelocityExtrude: debug 
+Add vx and vy as input parameters to diagnostic.