Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24012)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24013)
@@ -66,10 +66,11 @@
 
 	/*Allocate arrays*/
-	int numvertices = this->GetNumberOfVertices();
-	IssmDouble* lambdas = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* lambdas = xNew<IssmDouble>(NUM_VERTICES);
 
 	/* Start looping on the number of vertices: */
 	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -245,16 +246,17 @@
 
 	/*Allocate arrays*/
-	int numvertices = this->GetNumberOfVertices();
-	IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);
-	IssmDouble* eps_ef = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* eps_xx = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_yy = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_zz = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_xy = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_xz = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_yz = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* eps_ef = xNew<IssmDouble>(NUM_VERTICES);
 
 	/* Start looping on the number of vertices: */
 	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -388,6 +390,6 @@
 	_printf_("   sid: "<<this->sid<<"\n");
 	if(vertices){
-		int numvertices = this->GetNumberOfVertices();
-		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
+		const int NUM_VERTICES = this->GetNumberOfVertices();
+		for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
 	}
 	else _printf_("vertices = NULL\n");
@@ -427,13 +429,14 @@
 	if(!IsOnBase()) return;
 
-	int        numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES 					= this->GetNumberOfVertices();
+	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;
 
 	int        i;
-	IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* tmp=xNew<IssmDouble>(numvertices);
+	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
 	IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
 	IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
@@ -453,5 +456,5 @@
 	Gauss* gauss=this->NewGauss();
 	for(int month=0;month<12;month++){
-		for(int iv=0;iv<numvertices;iv++){
+		for(int iv=0;iv<NUM_VERTICES;iv++){
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
@@ -472,5 +475,5 @@
 
 	/*Compute the temperature and precipitation*/
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
 					Delta18oPresent, Delta18oLgm, Delta18oTime,
@@ -484,5 +487,5 @@
 	TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
 	for (int imonth=0;imonth<12;imonth++) {
-		for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -491,5 +494,5 @@
 			default: _error_("Not implemented yet");
 		}
-		for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -530,14 +533,15 @@
 	if(!IsOnBase()) return;
 
-	int        numvertices = this->GetNumberOfVertices();
-
-	int        i;
-	IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* tmp=xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES 					= this->GetNumberOfVertices();
+	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;
+
+	int        	i;
+	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
 	IssmDouble Delta18oTime;
 	IssmDouble dpermil,f;
@@ -583,5 +587,5 @@
 	Gauss* gauss=this->NewGauss();
 	for(int month=0;month<12;month++) {
-		for(int iv=0;iv<numvertices;iv++) {
+		for(int iv=0;iv<NUM_VERTICES;iv++) {
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
@@ -604,5 +608,5 @@
 
 	/*Compute the temperature and precipitation*/
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,isPrecipScaled,
 					f,&PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
@@ -615,5 +619,5 @@
 	TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
 	for (int imonth=0;imonth<12;imonth++) {
-		for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -622,5 +626,5 @@
 			default: _error_("Not implemented yet");
 		}
-		for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -662,5 +666,5 @@
 	/*Are we on the base? If not, return*/
 	if(!IsOnBase()) return;
-	int        numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 
 	int        i;
@@ -671,8 +675,8 @@
 	IssmDouble time;
 
-	IssmDouble*		smb		 = xNew<IssmDouble>(numvertices);
-	IssmDouble*		surf	 = xNew<IssmDouble>(numvertices);
-	IssmDouble*		accu	 = xNew<IssmDouble>(numvertices);
-	IssmDouble*		runoff = xNew<IssmDouble>(numvertices);
+	IssmDouble*		smb		= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble*		surf	= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble*		accu	= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble*		runoff 	= xNew<IssmDouble>(NUM_VERTICES);
 
 	/*Get material parameters :*/
@@ -695,5 +699,5 @@
 
 	/*Compute the temperature and precipitation*/
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		accu[iv]=max(0.,(accuref+(surf[iv]-accualti)*accugrad));
 		runoff[iv]=max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
@@ -914,6 +918,6 @@
 	_printf_("   sid: "<<this->sid<<"\n");
 	if(vertices){
-		int numvertices = this->GetNumberOfVertices();
-		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
+		const int NUM_VERTICES = this->GetNumberOfVertices();
+		for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
 	}
 	else _printf_("vertices = NULL\n");
@@ -1173,5 +1177,5 @@
 
 	/*Fetch number vertices for this element*/
-	int numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 
 	/*Checks in debugging mode*/
@@ -1180,5 +1184,5 @@
 	/* Start looping on the number of vertices: */
 	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 		input->GetInputValue(&pvalue[iv],gauss);
@@ -1196,5 +1200,5 @@
 
 	/*Fetch number vertices for this element*/
-	int numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 
 	/*Checks in debugging mode*/
@@ -1203,5 +1207,5 @@
 	/* Start looping on the number of vertices: */
 	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 		input->GetInputValue(&pvalue[iv],gauss,time);
@@ -1221,10 +1225,10 @@
 
 	/*Fetch number vertices for this element*/
-	int numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 
 	/* Start looping on the number of vertices: */
 	if (input){
 		Gauss* gauss=this->NewGauss();
-		for (int iv=0;iv<numvertices;iv++){
+		for (int iv=0;iv<NUM_VERTICES;iv++){
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&pvalue[iv],gauss);
@@ -1233,5 +1237,5 @@
 	}
 	else{
-		for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+		for(int iv=0;iv<NUM_VERTICES;iv++) pvalue[iv]=defaultvalue;
 	}
 }
@@ -1406,12 +1410,12 @@
 
 /* 	/\*Fetch number vertices for this element and allocate arrays*\/ */
-/* 	int numvertices = this->GetNumberOfVertices(); */
-/* 	int*        vertexpidlist = xNew<int>(numvertices); */
-/* 	IssmDouble* values        = xNew<IssmDouble>(numvertices); */
+/* 	const int NUM_VERTICES = this->GetNumberOfVertices(); */
+/* 	int*        vertexpidlist = xNew<int>(NUM_VERTICES); */
+/* 	IssmDouble* values        = xNew<IssmDouble>(NUM_VERTICES); */
 
 /* 	/\*Fill in values*\/ */
 /* 	this->GetVerticesPidList(vertexpidlist); */
 /* 	this->GetInputListOnVertices(values,input_enum); */
-/* 	vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); */
+/* 	vector->SetValues(NUM_VERTICES,vertexpidlist,values,INS_VAL); */
 
 /* 	/\*Clean up*\/ */
@@ -1424,5 +1428,6 @@
 
 	/*Fetch number vertices for this element and allocate arrays*/
-	int         numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	int         numnodes    = this->GetNumberOfNodes();
 	int*        doflist     = NULL;
@@ -1438,18 +1443,18 @@
 			break;
 		case VertexPIdEnum:
-			doflist = xNew<int>(numvertices);
-			values = xNew<IssmDouble>(numvertices);
+			doflist = xNew<int>(NUM_VERTICES);
+			values = xNew<IssmDouble>(NUM_VERTICES);
 			/*Fill in values*/
 			this->GetVerticesPidList(doflist);
 			this->GetInputListOnVertices(values,input_enum);
-			vector->SetValues(numvertices,doflist,values,INS_VAL);
+			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
 			break;
 		case VertexSIdEnum:
-			doflist = xNew<int>(numvertices);
-			values = xNew<IssmDouble>(numvertices);
+			doflist = xNew<int>(NUM_VERTICES);
+			values = xNew<IssmDouble>(NUM_VERTICES);
 			/*Fill in values*/
 			this->GetVerticesSidList(doflist);
 			this->GetInputListOnVertices(values,input_enum);
-			vector->SetValues(numvertices,doflist,values,INS_VAL);
+			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
 			break;
 		case NodesEnum:
@@ -1482,5 +1487,6 @@
 
 	/*Fetch number vertices for this element and allocate arrays*/
-	int         numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	int         numnodes    = this->GetNumberOfNodes();
 	int*        doflist     = NULL;
@@ -1489,18 +1495,18 @@
 	switch(type){
 	case VertexPIdEnum:
-		doflist = xNew<int>(numvertices);
-		values = xNew<IssmDouble>(numvertices);
+		doflist = xNew<int>(NUM_VERTICES);
+		values = xNew<IssmDouble>(NUM_VERTICES);
 		/*Fill in values*/
 		this->GetVerticesPidList(doflist);
 		this->GetInputListOnVerticesAtTime(values,input_enum,time);
-		vector->SetValues(numvertices,doflist,values,INS_VAL);
+		vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
 		break;
 	case VertexSIdEnum:
-		doflist = xNew<int>(numvertices);
-		values = xNew<IssmDouble>(numvertices);
+		doflist = xNew<int>(NUM_VERTICES);
+		values = xNew<IssmDouble>(NUM_VERTICES);
 		/*Fill in values*/
 		this->GetVerticesSidList(doflist);
 		this->GetInputListOnVerticesAtTime(values,input_enum,time);
-		vector->SetValues(numvertices,doflist,values,INS_VAL);
+		vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
 		break;
 	default:
@@ -1516,6 +1522,6 @@
 void       Element::GetVerticesLidList(int* lidlist){/*{{{*/
 
-	int numvertices = this->GetNumberOfVertices();
-	for(int i=0;i<numvertices;i++) lidlist[i]=vertices[i]->Lid();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+	for(int i=0;i<NUM_VERTICES;i++) lidlist[i]=vertices[i]->Lid();
 
 }
@@ -1523,6 +1529,6 @@
 void       Element::GetVerticesPidList(int* pidlist){/*{{{*/
 
-	int numvertices = this->GetNumberOfVertices();
-	for(int i=0;i<numvertices;i++) pidlist[i]=vertices[i]->Pid();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+	for(int i=0;i<NUM_VERTICES;i++) pidlist[i]=vertices[i]->Pid();
 
 }
@@ -1530,13 +1536,14 @@
 void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
 
-	int numvertices = this->GetNumberOfVertices();
-	for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+	for(int i=0;i<NUM_VERTICES;i++) connectivity[i]=this->vertices[i]->Connectivity();
 }
 /*}}}*/
 void       Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/
 
-	int         numvertices = this->GetNumberOfVertices();
-	IssmDouble* xyz_list    = xNew<IssmDouble>(numvertices*3);
-	::GetVerticesCoordinates(xyz_list,this->vertices,numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* xyz_list    = xNew<IssmDouble>(NUM_VERTICES*3);
+	::GetVerticesCoordinates(xyz_list,this->vertices,NUM_VERTICES);
 
 	*pxyz_list = xyz_list;
@@ -1545,6 +1552,6 @@
 void       Element::GetVerticesSidList(int* sidlist){/*{{{*/
 
-	int numvertices = this->GetNumberOfVertices();
-	for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+	for(int i=0;i<NUM_VERTICES;i++) sidlist[i]=this->vertices[i]->Sid();
 }
 /*}}}*/
@@ -1555,8 +1562,9 @@
 
 	/*Create list of x*/
-	int         numvertices = this->GetNumberOfVertices();
-	IssmDouble* x_list      = xNew<IssmDouble>(numvertices);
-
-	for(int i=0;i<numvertices;i++) x_list[i]=xyz_list[i*3+0];
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* x_list = xNew<IssmDouble>(NUM_VERTICES);
+
+	for(int i=0;i<NUM_VERTICES;i++) x_list[i]=xyz_list[i*3+0];
 	ValueP1OnGauss(&x,x_list,gauss);
 
@@ -1570,8 +1578,9 @@
 
 	/*Create list of y*/
-	int         numvertices = this->GetNumberOfVertices();
-	IssmDouble* y_list      = xNew<IssmDouble>(numvertices);
-
-	for(int i=0;i<numvertices;i++) y_list[i]=xyz_list[i*3+1];
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* y_list      = xNew<IssmDouble>(NUM_VERTICES);
+
+	for(int i=0;i<NUM_VERTICES;i++) y_list[i]=xyz_list[i*3+1];
 	ValueP1OnGauss(&y,y_list,gauss);
 
@@ -1585,8 +1594,9 @@
 
 	/*Create list of z*/
-	int         numvertices = this->GetNumberOfVertices();
-	IssmDouble* z_list      = xNew<IssmDouble>(numvertices);
-
-	for(int i=0;i<numvertices;i++) z_list[i]=xyz_list[i*3+2];
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* z_list      = xNew<IssmDouble>(NUM_VERTICES);
+
+	for(int i=0;i<NUM_VERTICES;i++) z_list[i]=xyz_list[i*3+2];
 	ValueP1OnGauss(&z,z_list,gauss);
 
@@ -1603,5 +1613,5 @@
 
 	/*Get number of vertices*/
-	int numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 	if(isautodiff){
 		int* N=NULL;
@@ -1617,6 +1627,6 @@
 
 		for(int n=0;n<N[control_index];n++){
-			for(int i=0;i<numvertices;i++){
-				indexing[i+n*numvertices]=this->vertices[i]->Sid() + start + n*M[control_index];
+			for(int i=0;i<NUM_VERTICES;i++){
+				indexing[i+n*NUM_VERTICES]=this->vertices[i]->Sid() + start + n*M[control_index];
 			}
 		}
@@ -1626,5 +1636,5 @@
 		parameters->FindParam(&M,ControlInputSizeMEnum);
 		/*get gradient indices*/
-		for(int i=0;i<numvertices;i++){
+		for(int i=0;i<NUM_VERTICES;i++){
 			indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
 		}
@@ -1720,12 +1730,13 @@
     if(vector_type==1){ //nodal vector
 
-        int         numvertices = this->GetNumberOfVertices();
-        int        *vertexids   = xNew<int>(numvertices);
-        IssmDouble *values      = xNew<IssmDouble>(numvertices);
+        const int NUM_VERTICES = this->GetNumberOfVertices();
+
+        int        *vertexids   = xNew<int>(NUM_VERTICES);
+        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
 
         /*Recover vertices ids needed to initialize inputs*/
         _assert_(iomodel->elements);
-        for(i=0;i<numvertices;i++){
-            vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
+        for(i=0;i<NUM_VERTICES;i++){
+            vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
         }
 
@@ -1736,5 +1747,5 @@
 		  }
 		  else if(M==iomodel->numberofvertices){
-            for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
+            for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
             this->AddInput(vector_enum,values,P1Enum);
         }
@@ -1745,5 +1756,5 @@
             TransientInput* transientinput=new TransientInput(vector_enum,times,N);
             for(t=0;t<N;t++){
-                for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
+                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
                 switch(this->ObjectEnum()){
                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
@@ -1833,9 +1844,10 @@
 
 	/*Intermediaries*/
-	int         numvertices = this->GetNumberOfVertices();
-	int        *vertexids   = xNew<int>(numvertices);
-	IssmDouble *values      = xNew<IssmDouble>(numvertices);
-	IssmDouble *values_min  = xNew<IssmDouble>(numvertices);
-	IssmDouble *values_max  = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	int        *vertexids   = xNew<int>(NUM_VERTICES);
+	IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble *values_min  = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble *values_max  = xNew<IssmDouble>(NUM_VERTICES);
 
 	/*Some sanity checks*/
@@ -1850,11 +1862,11 @@
 	/*Recover vertices ids needed to initialize inputs*/
 	_assert_(iomodel->elements);
-	for(int i=0;i<numvertices;i++){
-		vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
+	for(int i=0;i<NUM_VERTICES;i++){
+		vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
 	}
 
 	/*Are we in transient or static? */
 	if(M==iomodel->numberofvertices){
-		for(int i=0;i<numvertices;i++){
+		for(int i=0;i<NUM_VERTICES;i++){
 			values[i]=vector[vertexids[i]-1];
 			values_min[i] = min_vector[vertexids[i]-1];
@@ -1874,5 +1886,5 @@
 			   TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
 				for(int t=0;t<N;t++){
-                for(int i=0;i<numvertices;i++){
+                for(int i=0;i<NUM_VERTICES;i++){
 						values[i]=vector[N*(vertexids[i]-1)+t];
 						values_min[i] = min_vector[N*(vertexids[i]-1)+t];
@@ -1935,12 +1947,13 @@
     if(vector_type==1){ //nodal vector
 
-        int         numvertices = this->GetNumberOfVertices();
-        int        *vertexids   = xNew<int>(numvertices);
-        IssmDouble *values      = xNew<IssmDouble>(numvertices);
+        const int NUM_VERTICES = this->GetNumberOfVertices();
+
+        int        *vertexids   = xNew<int>(NUM_VERTICES);
+        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
 
         /*Recover vertices ids needed to initialize inputs*/
         _assert_(iomodel->elements);
-        for(i=0;i<numvertices;i++){
-            vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
+        for(i=0;i<NUM_VERTICES;i++){
+            vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
         }
 
@@ -1956,5 +1969,5 @@
 		  }
 		  else if(M==iomodel->numberofvertices){
-            for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
+            for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
 				switch(this->ObjectEnum()){
                     case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
@@ -1969,5 +1982,5 @@
             TransientInput* transientinput=new TransientInput(input_enum,times,N);
             for(t=0;t<N;t++){
-                for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
+                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
                 switch(this->ObjectEnum()){
                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break;
@@ -2161,8 +2174,9 @@
 	if(!this->IsIceInElement() || !this->IsFloating()) return;
 
-	int         numvertices = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	int         basinid,num_basins,M,N;
 	IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
-	IssmDouble  basalmeltrate[numvertices];
+	IssmDouble  basalmeltrate[NUM_VERTICES];
 	bool        islocal;
 	IssmDouble* delta_t = NULL;
@@ -2194,5 +2208,5 @@
 	/*Compute melt rate for Local and Nonlocal parameterizations*/
 	Gauss* gauss=this->NewGauss();
-	for(int i=0;i<numvertices;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		gauss->GaussVertex(i);
 		tf_input->GetInputValue(&tf,gauss);
@@ -2220,9 +2234,10 @@
 void       Element::LinearFloatingiceMeltingRate(){/*{{{*/
 
-	int numvertices      = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	IssmDouble  deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
-	IssmDouble* base     = xNew<IssmDouble>(numvertices);
-	IssmDouble* values   = xNew<IssmDouble>(numvertices);
-	IssmDouble time;
+	IssmDouble* base     		= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* values   		= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble 	time;
 
 	parameters->FindParam(&time,TimeEnum);
@@ -2231,8 +2246,8 @@
 	parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum,time);
 	parameters->FindParam(&upperwatermelt,BasalforcingsUpperwaterMeltingRateEnum,time);
-	_assert_(upperwaterel>deepwaterel); 
+	_assert_(upperwaterel>deepwaterel);
 
 	this->GetInputListOnVertices(base,BaseEnum);
-	for(int i=0;i<numvertices;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		if(base[i]>=upperwaterel){
 			values[i]=upperwatermelt;
@@ -2253,10 +2268,11 @@
 void       Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/
 
-	int numvertices      = this->GetNumberOfVertices();
-	IssmDouble* deepwatermelt     = xNew<IssmDouble>(numvertices);
-	IssmDouble* deepwaterel     = xNew<IssmDouble>(numvertices);
-	IssmDouble* upperwaterel     = xNew<IssmDouble>(numvertices);
-	IssmDouble* base     = xNew<IssmDouble>(numvertices);
-	IssmDouble* values   = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* deepwatermelt	= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* deepwaterel     = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* upperwaterel	= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* base			= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* values			= xNew<IssmDouble>(NUM_VERTICES);
 
 	this->GetInputListOnVertices(base,BaseEnum);
@@ -2265,5 +2281,5 @@
 	this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
 
-	for(int i=0;i<numvertices;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		if(base[i]>upperwaterel[i])      values[i]=0;
 		else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
@@ -2281,10 +2297,10 @@
 void       Element::MantlePlumeGeothermalFlux(){/*{{{*/
 
-	int numvertices      = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
 	IssmDouble  mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey;
 	IssmDouble  crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat;
 	IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
 	IssmDouble  x,y,z,c;
-	IssmDouble* values   = xNew<IssmDouble>(numvertices);
+	IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
 	IssmDouble *xyz_list = NULL;
 
@@ -2307,5 +2323,5 @@
 	e=pow(a*a-c*c,1./2.)/a;
 	A0=(1-pow(e,2.))/pow(e,3.)*(1./2.*log((1+e)/(1-e))-e);
-	for(int i=0;i<numvertices;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		y=xyz_list[i*3+0]-plumex;
 		z=xyz_list[i*3+1]-plumey;
@@ -2346,15 +2362,15 @@
 void       Element::MigrateGroundingLine(IssmDouble* phi_ungrounding){/*{{{*/
 
-	int         numvertices = this->GetNumberOfVertices();
+	const int  NUM_VERTICES = this->GetNumberOfVertices();
 	int        i,migration_style;
 	IssmDouble bed_hydro,yts;
 	IssmDouble rho_water,rho_ice,density;
-	IssmDouble* melting = xNew<IssmDouble>(numvertices);
-	IssmDouble* phi     = xNew<IssmDouble>(numvertices);
-	IssmDouble* h       = xNew<IssmDouble>(numvertices);
-	IssmDouble* s       = xNew<IssmDouble>(numvertices);
-	IssmDouble* b       = xNew<IssmDouble>(numvertices);
-	IssmDouble* r       = xNew<IssmDouble>(numvertices);
-	IssmDouble* sl      = xNew<IssmDouble>(numvertices);
+	IssmDouble* melting = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* phi     = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* h       = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* s       = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* b       = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* r       = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* sl      = xNew<IssmDouble>(NUM_VERTICES);
 
 	/*Recover info at the vertices: */
@@ -2372,5 +2388,5 @@
 
 	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUM_VERTICES;i++){
 		/* Contact FS*/
 		if(migration_style == ContactEnum){
@@ -2408,5 +2424,5 @@
 
 	/*Recalculate phi*/
-	for(i=0;i<numvertices;i++){
+	for(i=0;i<NUM_VERTICES;i++){
 		if(migration_style==SoftMigrationEnum){
 			bed_hydro=-density*h[i]+sl[i];
@@ -2438,10 +2454,10 @@
 /*}}}*/
 void       Element::MismipFloatingiceMeltingRate(){/*{{{*/
-
-	int numvertices      = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	IssmDouble  meltratefactor,thresholdthickness,upperdepthmelt;
-	IssmDouble* base     = xNew<IssmDouble>(numvertices);
-	IssmDouble* bed      = xNew<IssmDouble>(numvertices);
-	IssmDouble* values   = xNew<IssmDouble>(numvertices);
+	IssmDouble* base     = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* bed      = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
 
 	parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum);
@@ -2451,5 +2467,5 @@
 	this->GetInputListOnVertices(base,BaseEnum);
 	this->GetInputListOnVertices(bed,BedEnum);
-	for(int i=0;i<numvertices;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		if(base[i]>upperdepthmelt){
 			values[i]=0;
@@ -2470,14 +2486,15 @@
 	if(!IsOnBase()) return;
 
-	int        numvertices = this->GetNumberOfVertices();
-
-	int        i;
-	IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* tmp=xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES 					= this->GetNumberOfVertices();
+	const int NUM_VERTICES_MONTHS_PER_YEAR 	= NUM_VERTICES * 12;
+
+	int i;
+	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
 	IssmDouble TdiffTime,PfacTime;
 	IssmDouble time,yts,time_yr;
@@ -2494,5 +2511,5 @@
 	Gauss* gauss=this->NewGauss();
 	for(int month=0;month<12;month++) {
-		for(int iv=0;iv<numvertices;iv++) {
+		for(int iv=0;iv<NUM_VERTICES;iv++) {
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
@@ -2511,5 +2528,5 @@
 
 	/*Compute the temperature and precipitation*/
-	for(int iv=0;iv<numvertices;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
 					&PrecipitationsLgm[iv*12],&PrecipitationsPresentday[iv*12],
@@ -2522,5 +2539,5 @@
 	TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
 	for (int imonth=0;imonth<12;imonth++) {
-		for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -2529,5 +2546,5 @@
 			default: _error_("Not implemented yet");
 		}
-		for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
+		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
 		switch(this->ObjectEnum()){
 			case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
@@ -2616,5 +2633,5 @@
 
 	this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));
-	
+
 }/*}}}*/
 void       Element::PicoUpdateBox(int loopboxid){/*{{{*/
@@ -2626,5 +2643,6 @@
 	if(loopboxid!=boxid) return;
 
-	int        NUMVERTICES = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	int        basinid, maxbox, num_basins, numnodes, M;
 	IssmDouble gamma_T, overturning_coeff, thickness;
@@ -2664,8 +2682,8 @@
 	IssmDouble g1               = area_boxi*gamma_T;
 
-	IssmDouble basalmeltrates_shelf[NUMVERTICES];
-	IssmDouble potential_pressure_melting_point[NUMVERTICES];
-	IssmDouble Tocs[NUMVERTICES];
-	IssmDouble Socs[NUMVERTICES];
+	IssmDouble basalmeltrates_shelf[NUM_VERTICES];
+	IssmDouble potential_pressure_melting_point[NUM_VERTICES];
+	IssmDouble Tocs[NUM_VERTICES];
+	IssmDouble Socs[NUM_VERTICES];
 
 	/* First box calculations */
@@ -2677,9 +2695,9 @@
 		this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
 		IssmDouble s1 = soc_farocean/(nu*lambda);
-		IssmDouble overturnings[NUMVERTICES];
+		IssmDouble overturnings[NUM_VERTICES];
 
 		/* Start looping on the number of verticies and calculate ocean vars */
 		Gauss* gauss=this->NewGauss();
-		for(int i=0;i<NUMVERTICES;i++){
+		for(int i=0;i<NUM_VERTICES;i++){
 			gauss->GaussVertex(i);
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -2724,5 +2742,5 @@
 		/* Start looping on the number of verticies and calculate ocean vars */
 		Gauss* gauss=this->NewGauss();
-		for(int i=0;i<NUMVERTICES;i++){
+		for(int i=0;i<NUM_VERTICES;i++){
 			gauss->GaussVertex(i);
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -2748,5 +2766,5 @@
 	/*Cleanup and return*/
 	xDelete<IssmDouble>(boxareas);
-	
+
 }/*}}}*/
 void       Element::PicoComputeBasalMelt(){/*{{{*/
@@ -2754,5 +2772,6 @@
 	if(!this->IsIceInElement() || !this->IsFloating()) return;
 
-	int        NUMVERTICES = this->GetNumberOfVertices();
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
 	IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
 	IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
@@ -2774,5 +2793,5 @@
 
 	/*Define arrays*/
-	IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
+	IssmDouble basalmeltrates_shelf[NUM_VERTICES];  //Basal melt-rate
 
 	/*Polynomial coefficients*/
@@ -2802,5 +2821,5 @@
 	/*Loop over the number of vertices in this element*/
 	Gauss* gauss=this->NewGauss();
-	for(int i=0;i<NUMVERTICES;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		gauss->GaussVertex(i);
 
@@ -2851,22 +2870,23 @@
 	/*Cleanup and return*/
 	delete gauss;
-		
+
 }/*}}}*/
 void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
 
-	int  numvertices = this->GetNumberOfVertices();
-
-	int        i;
-	IssmDouble* agd=xNew<IssmDouble>(numvertices); // surface mass balance
-	IssmDouble* melt=xNew<IssmDouble>(numvertices); // surface mass balance
-	IssmDouble* accu=xNew<IssmDouble>(numvertices); // surface mass balance
-	IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
-	IssmDouble* tmp=xNew<IssmDouble>(numvertices);
-	IssmDouble* h=xNew<IssmDouble>(numvertices);
-	IssmDouble* s=xNew<IssmDouble>(numvertices);
-	IssmDouble* s0p=xNew<IssmDouble>(numvertices);
-	IssmDouble* s0t=xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES 		= this->GetNumberOfVertices();
+	const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 12;
+
+	int  		i;
+	IssmDouble* agd=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
+	IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
+	IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
+	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
+	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* h=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);
 	IssmDouble pddsnowfac = -1.;
 	IssmDouble pddicefac = -1.;
@@ -2895,5 +2915,5 @@
 	Gauss* gauss=this->NewGauss();
 	for(int month=0;month<12;month++) {
-		for(int iv=0;iv<numvertices;iv++) {
+		for(int iv=0;iv<NUM_VERTICES;iv++) {
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,time_yr+month/12.*yts);
@@ -2932,5 +2952,5 @@
 
 	/*measure the surface mass balance*/
-	for (int iv = 0; iv<numvertices; iv++){
+	for (int iv = 0; iv<NUM_VERTICES; iv++){
 		gauss->GaussVertex(iv);
 		pddsnowfac=-1.;
@@ -2952,9 +2972,9 @@
 	// TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
 	// for (int imonth=0;imonth<12;imonth++) {
-	//   for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
+	//   for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
 	//   TriaInput* newmonthinput1 = new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum);
 	//   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
 	//
-	//   for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
+	//   for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
 	//   TriaInput* newmonthinput2 = new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum);
 	//   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
@@ -3042,20 +3062,21 @@
 	/* General FIXMEs: get Tmelting point, pddicefactor, pddsnowfactor, sigma from parameters/user input */
 
-	int  numvertices = this->GetNumberOfVertices();
-
-	int        i;
-	IssmDouble* smb=xNew<IssmDouble>(numvertices);		// surface mass balance
-	IssmDouble* melt=xNew<IssmDouble>(numvertices);		// melting comp. of surface mass balance
-	IssmDouble* accu=xNew<IssmDouble>(numvertices);		// accuumulation comp. of surface mass balance
-	IssmDouble* melt_star=xNew<IssmDouble>(numvertices);
-	IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
-	IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
-	IssmDouble* s=xNew<IssmDouble>(numvertices);			// actual surface height
-	IssmDouble* s0p=xNew<IssmDouble>(numvertices);		// reference elevation for precip.
-	IssmDouble* s0t=xNew<IssmDouble>(numvertices);		// reference elevation for temperature
-	IssmDouble* smbcorr=xNew<IssmDouble>(numvertices); // surface mass balance correction; will be added after pdd call
-	IssmDouble* p_ampl=xNew<IssmDouble>(numvertices);	// precip anomaly
-	IssmDouble* t_ampl=xNew<IssmDouble>(numvertices);	// remperature anomaly
+	const int NUM_VERTICES 		= this->GetNumberOfVertices();
+	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;
+
+	int        	i;
+	IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);		// surface mass balance
+	IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);		// melting comp. of surface mass balance
+	IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES);		// accuumulation comp. of surface mass balance
+	IssmDouble* melt_star=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
+	IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
+	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);			// actual surface height
+	IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);		// reference elevation for precip.
+	IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);		// reference elevation for temperature
+	IssmDouble* smbcorr=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance correction; will be added after pdd call
+	IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);	// precip anomaly
+	IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);	// remperature anomaly
 	IssmDouble rho_water,rho_ice,desfac,rlaps;
 	IssmDouble inv_twelve=1./12.;								//factor for monthly average
@@ -3088,5 +3109,5 @@
 	Gauss* gauss=this->NewGauss();
 	for(int month=0;month<12;month++){
-		for(int iv=0;iv<numvertices;iv++){
+		for(int iv=0;iv<NUM_VERTICES;iv++){
 			gauss->GaussVertex(iv);
 			input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,(month+1)/12.*yts);
@@ -3106,5 +3127,5 @@
 
 	/*measure the surface mass balance*/
-	for (int iv = 0; iv<numvertices; iv++){
+	for (int iv = 0; iv<NUM_VERTICES; iv++){
 		smb[iv]=PddSurfaceMassBalanceSicopolis(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
 					&melt[iv], &accu[iv], &melt_star[iv], &t_ampl[iv], &p_ampl[iv], yts, s[iv],
@@ -3324,15 +3345,16 @@
 		}
 		case P1Enum:{
-			int         numvertices = this->GetNumberOfVertices();
-			IssmDouble *values      = xNew<IssmDouble>(numvertices);
-			int        *connectivity= xNew<int>(numvertices);
-			int        *sidlist     = xNew<int>(numvertices);
+			const int NUM_VERTICES = this->GetNumberOfVertices();
+
+			IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
+			int        *connectivity= xNew<int>(NUM_VERTICES);
+			int        *sidlist     = xNew<int>(NUM_VERTICES);
 
 			this->GetVerticesSidList(sidlist);
 			this->GetVerticesConnectivityList(connectivity);
 			this->GetInputListOnVertices(values,output_enum);
-			for(int i=0;i<numvertices;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
-
-			vector->SetValues(numvertices,sidlist,values,ADD_VAL);
+			for(int i=0;i<NUM_VERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
+
+			vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL);
 
 			xDelete<IssmDouble>(values);
@@ -3421,25 +3443,26 @@
 	if (!IsOnSurface()) return;
 
-	int  numvertices = this->GetNumberOfVertices();
-
-	IssmDouble* s=xNew<IssmDouble>(numvertices);
-	IssmDouble* s0gcm=xNew<IssmDouble>(numvertices);
-	IssmDouble* st=xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES 					= this->GetNumberOfVertices();
+	const int NUM_VERTICES_DAYS_PER_YEAR 	= NUM_VERTICES * 365;
+
+	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* s0gcm=xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble* st=xNew<IssmDouble>(NUM_VERTICES);
 
 	// daily forcing inputs
-	IssmDouble* dailyrainfall=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailysnowfall=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailydlradiation=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailydsradiation=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailywindspeed=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailypressure=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailyairdensity=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailyairhumidity=xNew<IssmDouble>(365*numvertices);
-	IssmDouble* dailytemperature=xNew<IssmDouble>(365*numvertices);
+	IssmDouble* dailyrainfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailysnowfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailydlradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailydsradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailywindspeed=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailypressure=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailyairdensity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailyairhumidity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
+	IssmDouble* dailytemperature=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
 	// daily outputs
-	IssmDouble* tsurf_out=xNew<IssmDouble>(numvertices); memset(tsurf_out, 0., numvertices*sizeof(IssmDouble));
-	IssmDouble* smb_out=xNew<IssmDouble>(numvertices); memset(smb_out, 0., numvertices*sizeof(IssmDouble));
-	IssmDouble* saccu_out=xNew<IssmDouble>(numvertices); memset(saccu_out, 0., numvertices*sizeof(IssmDouble));
-	IssmDouble* smelt_out=xNew<IssmDouble>(numvertices); memset(smelt_out, 0., numvertices*sizeof(IssmDouble));
+	IssmDouble* tsurf_out=xNew<IssmDouble>(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble));
+	IssmDouble* smb_out=xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
+	IssmDouble* saccu_out=xNew<IssmDouble>(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble));
+	IssmDouble* smelt_out=xNew<IssmDouble>(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble));
 
 	IssmDouble rho_water,rho_ice,desfac,rlaps,rdl;
@@ -3476,5 +3499,5 @@
 	Gauss* gauss=this->NewGauss();
 	for (int iday = 0; iday < 365; iday++){
-		for(int iv=0;iv<numvertices;iv++) {
+		for(int iv=0;iv<NUM_VERTICES;iv++) {
 			gauss->GaussVertex(iv);
 			/* get forcing */
@@ -3508,7 +3531,7 @@
 	}
 
-	for (int iv = 0; iv<numvertices; iv++){
+	for (int iv = 0; iv<NUM_VERTICES; iv++){
 		/* call semic */
-		run_semic_(&dailysnowfall[iv*365], &dailyrainfall[iv*365], &dailydsradiation[iv*365], &dailydlradiation[iv*365], 
+		run_semic_(&dailysnowfall[iv*365], &dailyrainfall[iv*365], &dailydsradiation[iv*365], &dailydlradiation[iv*365],
 					&dailywindspeed[iv*365], &dailypressure[iv*365], &dailyairdensity[iv*365], &dailyairhumidity[iv*365], &dailytemperature[iv*365],
 					&tsurf_out[iv], &smb_out[iv], &saccu_out[iv], &smelt_out[iv]);
@@ -3516,5 +3539,5 @@
 
 	switch(this->ObjectEnum()){
-		case TriaEnum:  
+		case TriaEnum:
 			this->inputs->AddInput(new TriaInput(TemperatureSEMICEnum,&tsurf_out[0],P1Enum)); // TODO add TemperatureSEMICEnum to EnumDefinitions
 			this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb_out[0],P1Enum));
@@ -3525,5 +3548,5 @@
 			// TODO
 			break;
-		case TetraEnum: 
+		case TetraEnum:
 			// TODO
 			break;
@@ -3540,5 +3563,5 @@
 	xDelete<IssmDouble>(dailypressure);
 	xDelete<IssmDouble>(dailyairdensity);
-	xDelete<IssmDouble>(dailyairhumidity);	
+	xDelete<IssmDouble>(dailyairhumidity);
 	xDelete<IssmDouble>(dailypressure);
 	xDelete<IssmDouble>(dailytemperature);
@@ -3548,5 +3571,5 @@
 	xDelete<IssmDouble>(tsurf_out);
 	xDelete<IssmDouble>(s);
-	xDelete<IssmDouble>(st);	
+	xDelete<IssmDouble>(st);
 	xDelete<IssmDouble>(s0gcm);
 }
@@ -4161,6 +4184,7 @@
 
 	/*Fetch number vertices and allocate memory*/
-	int         numvertices  = this->GetNumberOfVertices();
-	IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES  = this->GetNumberOfVertices();
+
+	IssmDouble* maxprincipal = xNew<IssmDouble>(NUM_VERTICES);
 
 	/*Retrieve all inputs and parameters*/
@@ -4180,5 +4204,5 @@
 	/*loop over vertices: */
 	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -4503,6 +4527,7 @@
 
 	/*Fetch number vertices and allocate memory*/
-	int         numvertices    = this->GetNumberOfVertices();
-	IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
+	const int NUM_VERTICES = this->GetNumberOfVertices();
+
+	IssmDouble* viscousheating = xNew<IssmDouble>(NUM_VERTICES);
 
 	/*Retrieve all inputs and parameters*/
@@ -4514,5 +4539,5 @@
 	/*loop over vertices: */
 	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
+	for (int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);