Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3595)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3596)
@@ -110,4 +110,6 @@
 /*FUNCTION Tria::Configure {{{1*/
 void  Tria::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){
+
+	int i;
 
 	DataSet* loadsin=NULL;
@@ -511,4 +513,10 @@
 	  /* 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);
+
+#ifdef _ISSM_DEBUG_ 
+	for (i=0;i<num_gauss;i++){
+		printf("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),*(gauss_weights+i));
+	}
+#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -3798,4 +3806,7 @@
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_ 
+		printf("Element id %i Jacobian determinant: %g\n",GetId(),Jdet);
+#endif
 
 		/* Get nodal functions value at gaussian point:*/
@@ -3923,4 +3934,10 @@
 	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
 
+	#ifdef _ISSM_DEBUG_ 
+	for (i=0;i<3;i++){
+		printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh3[0][i],dh1dh3[1][i]);
+	}
+	#endif
+
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
@@ -3957,4 +3974,10 @@
 	/*Get dh1dh2dh3 in actual coordinate system: */
 	GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
+
+#ifdef _ISSM_DEBUG_ 
+	for (i=0;i<3;i++){
+		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
+	}
+#endif
 
 	/*Build B_prog: */
@@ -4535,4 +4558,9 @@
 	/* 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, 4);
+#ifdef _ISSM_DEBUG_ 
+	for (i=0;i<num_gauss;i++){
+		printf("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),*(gauss_weights+i));
+	}
+#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -4564,4 +4592,7 @@
 		/* Get nodal functions value at gaussian point:*/
 		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_
+		printf("viscositycomp %g thickness %g dvx [%g %g] dvy [%g %g]  dadjx [%g %g] dadjy[%g %g]\n",viscosity_complement,thickness,dvx[0],dvx[1],dvy[0],dvy[1],dadjx[0],dadjx[1],dadjy[0],dadjy[1]);
+#endif
 
 		/*Get nodal functions derivatives*/
@@ -4768,16 +4799,28 @@
 		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
 		GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
+		#ifdef _ISSM_DEBUG_ 
+			printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
+		#endif
 
 		/*recover lambda and mu: */
 		GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
 		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
+		#ifdef _ISSM_DEBUG_ 
+			printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
+		#endif
 			
 		/*recover vx and vy: */
 		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
 		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+		#ifdef _ISSM_DEBUG_ 
+			printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
+		#endif
 
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-	
+		#ifdef _ISSM_DEBUG_ 
+		printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
+		#endif
+		
 		/* Get nodal functions value at gaussian point:*/
 		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
@@ -4990,4 +5033,7 @@
 		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
 		GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_ 
+		printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
+#endif
 
 		/*recover lambda mu and xi: */
@@ -4995,4 +5041,7 @@
 		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
 		GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_ 
+		printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
+#endif
 
 		/*recover vx vy and vz: */
@@ -5000,7 +5049,20 @@
 		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
 		GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_ 
+		printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
+
+		/*Get normal vecyor to the bed */
+		SurfaceNormal(&surface_normal[0],xyz_list);
+
+		bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
+		bed_normal[1]=-surface_normal[1];
+		bed_normal[2]=-surface_normal[2];
+#endif
 
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+#ifdef _ISSM_DEBUG_ 
+		printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
+#endif
 
 		/* Get nodal functions value at gaussian point:*/