Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15713)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15714)
@@ -4803,9 +4803,8 @@
 	const int    numdof=NDOF1*NUMVERTICES;
 
-	int          i;
-	int*         doflist=NULL;
-	IssmDouble       values[numdof];
-	IssmDouble       enthalpy;
-	GaussPenta   *gauss=NULL;
+	int*        doflist=NULL;
+	IssmDouble  values[numdof];
+	IssmDouble  enthalpy;
+	GaussPenta *gauss=NULL;
 
 	/*Get dof list: */
@@ -4814,5 +4813,5 @@
 
 	gauss=new GaussPenta();
-	for(i=0;i<NUMVERTICES;i++){
+	for(int i=0;i<NUMVERTICES;i++){
 		/*Recover temperature*/
 		gauss->GaussVertex(i);
@@ -5320,5 +5319,4 @@
 	int vnumnodes = this->NumberofNodesVelocity();
 	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
 
 	/*Prepare coordinate system list*/
@@ -5329,5 +5327,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -7795,7 +7793,4 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticVertVolume(void){
 
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
 	/*Intermediaries */
 	IssmDouble  Jdet;
@@ -7898,7 +7893,4 @@
 /*FUNCTION Penta::CreatePVectorCouplingSSAFSViscous {{{*/
 ElementVector* Penta::CreatePVectorCouplingSSAFSViscous(void){
-
-	/*Constants*/
-	const int   numdof=NUMVERTICES*NDOF4;
 
 	/*Intermediaries */
@@ -7960,7 +7952,4 @@
 /*FUNCTION Penta::CreatePVectorCouplingSSAFSFriction{{{*/
 ElementVector* Penta::CreatePVectorCouplingSSAFSFriction(void){
-
-	/*Constants*/
-	const int numdof=NUMVERTICES*NDOF4;
 
 	/*Intermediaries*/
@@ -8050,7 +8039,4 @@
 ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){
 
-	/*Constants*/
-	const int   numdof=NUMVERTICES*NDOF4;
-
 	/*Intermediaries */
 	int         i;
@@ -8111,7 +8097,4 @@
 /*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/
 ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){
-
-	/*Constants*/
-	const int numdof=NUMVERTICES*NDOF4;
 
 	/*Intermediaries*/
@@ -8279,5 +8262,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes(); _assert_(numnodes==6); 
-	int numdof   = numnodes*NDOF2;
 
 	/*Initialize Element vector*/
@@ -8406,5 +8388,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
 
 	/*Initialize vector*/
@@ -8479,5 +8460,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
 
 	/*Initialize Element vector and other vectors*/
@@ -8586,5 +8566,4 @@
 	int vnumnodes = this->NumberofNodesVelocity();
 	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
 
 	/*Prepare coordinate system list*/
@@ -8595,5 +8574,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -8642,7 +8621,4 @@
 /*FUNCTION Penta::PVectorGLSstabilization{{{*/
 void Penta::PVectorGLSstabilization(ElementVector* pe){
-
-	/*Constants*/
-	const int numdof=NDOF4*NUMVERTICES;
 
 	/*Intermediaries*/
@@ -8753,5 +8729,4 @@
 	int vnumnodes = this->NumberofNodesVelocity();
 	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
 
 	/*Prepare coordinate system list*/
@@ -8762,5 +8737,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -8826,5 +8801,4 @@
 	int vnumnodes = this->NumberofNodesVelocity();
 	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
 
 	/*Prepare coordinate system list*/
@@ -8835,5 +8809,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15713)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15714)
@@ -5004,11 +5004,8 @@
 ElementVector* Tria::CreatePVectorAdjointHoriz(void){
 
-	/*Constants*/
-	const int    numdof=NDOF2*NUMVERTICES;
-
 	/*Intermediaries */
-	int         i,resp;
-	int        *responses=NULL;
-	int         num_responses;
+	int        i,resp;
+	int       *responses=NULL;
+	int        num_responses;
 	IssmDouble Jdet;
 	IssmDouble obs_velocity_mag,velocity_mag;
@@ -5016,12 +5013,15 @@
 	IssmDouble epsvel=2.220446049250313e-16;
 	IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
-	IssmDouble scalex=0,scaley=0,scale=0,S=0;
+	IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.;
 	IssmDouble vx,vy,vxobs,vyobs,weight;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[3];
-	GaussTria* gauss=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
 
 	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -5043,5 +5043,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(4);
+	GaussTria* gauss=new GaussTria(4);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -5074,5 +5074,5 @@
 					 *        du     obs
 					 */
-					for (i=0;i<NUMVERTICES;i++){
+					for(i=0;i<numnodes;i++){
 						dux=vxobs-vx;
 						duy=vyobs-vy;
@@ -5093,5 +5093,5 @@
 					 *               obs
 					 */
-					for (i=0;i<NUMVERTICES;i++){
+					for(i=0;i<numnodes;i++){
 						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
 						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
@@ -5114,5 +5114,5 @@
 					 *            
 					 */
-					for (i=0;i<NUMVERTICES;i++){
+					for(i=0;i<numnodes;i++){
 						velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
 						obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
@@ -5134,5 +5134,5 @@
 					 *        du      S  2 sqrt(...)           obs
 					 */
-					for (i=0;i<NUMVERTICES;i++){
+					for(i=0;i<numnodes;i++){
 						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
 						dux=scale*(vxobs-vx);
@@ -5152,5 +5152,5 @@
 					 *        du                         |u| + eps  |u|                           u + eps
 					 */
-					for (i=0;i<NUMVERTICES;i++){
+					for(i=0;i<numnodes;i++){
 						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
 						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
@@ -5181,6 +5181,7 @@
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<int>(responses);
 	delete gauss;
-	xDelete<int>(responses);
 	return pe;
 }
@@ -6494,13 +6495,12 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 
-	/*Fetch number of nodes and dof for this finite element*/
+	/*Fetch number of nodes for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
 	IssmDouble     D[2][2];
 
@@ -6538,7 +6538,7 @@
 		D_scalar=gauss->weight*Jdettria;
 
-		TripleMultiply(basis,1,numdof,1,
+		TripleMultiply(basis,1,numnodes,1,
 					&D_scalar,1,1,0,
-					basis,1,numdof,0,
+					basis,1,numnodes,0,
 					&Ke->values[0],1);
 
@@ -6554,14 +6554,14 @@
 		D[1][1]=D_scalar*dvydy;
 		D[1][0]=0.;
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					B,2,numdof,0,
+					B,2,numnodes,0,
 					&Ke->values[0],1);
 
 		D[0][0]=D_scalar*vx;
 		D[1][1]=D_scalar*vy;
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					Bprime,2,numdof,0,
+					Bprime,2,numnodes,0,
 					&Ke->values[0],1);
 
@@ -6588,7 +6588,7 @@
 			D[0][1]=D_scalar*D[0][1];
 			D[1][1]=D_scalar*D[1][1];
-			TripleMultiply(Bprime,2,numdof,1,
+			TripleMultiply(Bprime,2,numnodes,1,
 						&D[0][0],2,2,0,
-						Bprime,2,numdof,0,
+						Bprime,2,numnodes,0,
 						&Ke->values[0],1);
 		}
@@ -6607,5 +6607,5 @@
 
 	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
+	const int    numnodes=NUMVERTICES;
 
 	/*Intermediaries */
@@ -6653,7 +6653,7 @@
 		DL_scalar=gauss->weight*Jdettria;
 
-		TripleMultiply(&basis[0],1,numdof,1,
+		TripleMultiply(&basis[0],1,numnodes,1,
 					&DL_scalar,1,1,0,
-					&basis[0],1,numdof,0,
+					&basis[0],1,numnodes,0,
 					&Ke->values[0],1);
 
@@ -6667,7 +6667,7 @@
 		DLprime[1][1]=DL_scalar*vy;
 
-		TripleMultiply( &B[0][0],2,numdof,1,
+		TripleMultiply( &B[0][0],2,numnodes,1,
 					&DLprime[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
+					&Bprime[0][0],2,numnodes,0,
 					&Ke->values[0],1);
 	}
@@ -6701,5 +6701,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element vector and other vectors*/
@@ -6733,5 +6732,5 @@
 		 basal_melting_correction_g=0.;
 
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
 	}
 
@@ -6746,5 +6745,5 @@
 
 	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
+	const int    numnodes=NUMVERTICES;
 
 	/*Intermediaries */
@@ -6778,5 +6777,5 @@
 		thickness_input->GetInputValue(&thickness_g,gauss);
 
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
 	}
 
@@ -6957,11 +6956,10 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
 	IssmDouble     D[2][2];
 
@@ -7005,14 +7003,14 @@
 		D[1][1]=D_scalar*dvydy;
 		D[1][1]=0.;
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					B,2,numdof,0,
+					B,2,numnodes,0,
 					&Ke->values[0],1);
 
 		D[0][0]=D_scalar*vx;
 		D[1][1]=D_scalar*vy;
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					Bprime,2,numdof,0,
+					Bprime,2,numnodes,0,
 					&Ke->values[0],1);
 
@@ -7039,7 +7037,7 @@
 			D[0][1]=D_scalar*D[0][1];
 			D[1][1]=D_scalar*D[1][1];
-			TripleMultiply(Bprime,2,numdof,1,
+			TripleMultiply(Bprime,2,numnodes,1,
 						&D[0][0],2,2,0,
-						Bprime,2,numdof,0,
+						Bprime,2,numnodes,0,
 						&Ke->values[0],1);
 		}
@@ -7058,5 +7056,5 @@
 
 	/*Constants*/
-	const int  numdof=NDOF1*NUMVERTICES;
+	const int  numnodes=NUMVERTICES;
 
 	/*Intermediaries*/
@@ -7097,7 +7095,7 @@
 		DL[1][1]=DL_scalar*vy;
 
-		TripleMultiply( &B[0][0],2,numdof,1,
+		TripleMultiply( &B[0][0],2,numnodes,1,
 					&DL[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
+					&Bprime[0][0],2,numnodes,0,
 					&Ke->values[0],1);
 	}
@@ -7131,5 +7129,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element vector and other vectors*/
@@ -7156,5 +7153,5 @@
 		GetNodalFunctions(basis,gauss);
 
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
 	}
 
@@ -7169,5 +7166,5 @@
 
 	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
+	const int    numnodes=NUMVERTICES;
 
 	/*Intermediaries */
@@ -7200,5 +7197,5 @@
 		GetNodalFunctions(&basis[0],gauss);
 
-		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
+		for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
 	}