Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 25904)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 25905)
@@ -651,24 +651,20 @@
 	const int NUM_VERTICES = this->GetNumberOfVertices();
 
-	int        i;
-	IssmDouble accuref, runoffref; //reference values at given altitude
-	IssmDouble accualti, runoffalti; //reference altitudes
+	IssmDouble accuref, runoffref;  //reference values at given altitude
 	IssmDouble accugrad, runoffgrad; //gradients from reference altitude
-	IssmDouble rho_water, rho_ice;
-	IssmDouble time;
-
-	IssmDouble*		smb	 = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble*		surf	 = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble*		accu	 = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble*		runoff = xNew<IssmDouble>(NUM_VERTICES);
+
+	IssmDouble smb[MAXVERTICES];
+	IssmDouble surf[MAXVERTICES];
+	IssmDouble accu[MAXVERTICES];
+	IssmDouble runoff[MAXVERTICES];
 
 	/*Get material parameters :*/
-	rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
-	rho_ice=this->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_water = this->FindParam(MaterialsRhoSeawaterEnum);
+	IssmDouble rho_ice   = this->FindParam(MaterialsRhoIceEnum);
 
 	/*Recover parameters*/
-	parameters->FindParam(&time,TimeEnum);
-	parameters->FindParam(&accualti,SmbAccualtiEnum);
-	parameters->FindParam(&runoffalti,SmbRunoffaltiEnum);
+	IssmDouble time       = this->FindParam(TimeEnum);
+	IssmDouble accualti   = this->FindParam(SmbAccualtiEnum);
+	IssmDouble runoffalti = this->FindParam(SmbRunoffaltiEnum);
 
 	/*Recover reference values at current time*/
@@ -683,7 +679,7 @@
 	/*Compute the temperature and precipitation*/
 	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));
-		smb[iv]=(accu[iv]-runoff[iv])*rho_ice/rho_water;
+		accu[iv]   = max(0.,(accuref+(surf[iv]-accualti)*accugrad));
+		runoff[iv] = max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
+		smb[iv]    = (accu[iv]-runoff[iv])*rho_ice/rho_water;
 	}
 
@@ -701,9 +697,4 @@
 	default: _error_("Not implemented yet");
 	}
-	/*clean-up*/
-	xDelete<IssmDouble>(surf);
-	xDelete<IssmDouble>(accu);
-	xDelete<IssmDouble>(runoff);
-	xDelete<IssmDouble>(smb);
 }
 /*}}}*/
@@ -1293,70 +1284,58 @@
 }
 /*}}}*/
-/* void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/\*{{{*\/ */
-
-/* 	/\*Fetch number vertices for this element and allocate arrays*\/ */
-/* 	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(NUM_VERTICES,vertexpidlist,values,INS_VAL); */
-
-/* 	/\*Clean up*\/ */
-/* 	xDelete<int>(vertexpidlist); */
-/* 	xDelete<IssmDouble>(values); */
-
-/* } */
-/* /\*}}}*\/ */
 void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/
 
-	/*Fetch number vertices for this element and allocate arrays*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	int         numnodes    = this->GetNumberOfNodes();
-	int*        doflist     = NULL;
-	IssmDouble  value;
-	IssmDouble* values      = NULL;
-	Input*     input       = NULL;
 
 	switch(type){
-		case ElementSIdEnum:
-			input=this->GetInput(input_enum); _assert_(input);
+		case ElementSIdEnum:{
+         IssmDouble  value;
+			Input* input=this->GetInput(input_enum); _assert_(input);
 			input->GetInputAverage(&value);
 			vector->SetValue(this->sid,value,INS_VAL);
+                          }
 			break;
-		case VertexPIdEnum:
-			doflist = xNew<int>(NUM_VERTICES);
-			values = xNew<IssmDouble>(NUM_VERTICES);
+		case VertexPIdEnum:{
+         int        doflist[MAXVERTICES];
+         IssmDouble values[MAXVERTICES];
+			const int  NUM_VERTICES = this->GetNumberOfVertices();
 			/*Fill in values*/
-			this->GetVerticesPidList(doflist);
-			this->GetInputListOnVertices(values,input_enum);
+			this->GetVerticesPidList(&doflist[0]);
+			this->GetInputListOnVertices(&values[0],input_enum);
 			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
+                         }
 			break;
-		case VertexSIdEnum:
-			doflist = xNew<int>(NUM_VERTICES);
-			values = xNew<IssmDouble>(NUM_VERTICES);
+		case VertexSIdEnum:{
+         int        doflist[MAXVERTICES];
+         IssmDouble values[MAXVERTICES];
+			const int  NUM_VERTICES = this->GetNumberOfVertices();
 			/*Fill in values*/
 			this->GetVerticesSidList(doflist);
 			this->GetInputListOnVertices(values,input_enum);
 			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
+                         }
 			break;
-		case NodesEnum:
-			doflist = xNew<int>(numnodes);
-			values = xNew<IssmDouble>(numnodes);
+		case NodesEnum:{
+         int         numnodes = this->GetNumberOfNodes();
+         int        *doflist  = xNew<int>(numnodes);
+         IssmDouble *values   = xNew<IssmDouble>(numnodes);
 			/*Fill in values*/
 			this->GetInputListOnNodes(values,input_enum);
 			this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 			vector->SetValues(numnodes,doflist,values,INS_VAL);
+         xDelete<int>(doflist);
+         xDelete<IssmDouble>(values);
+                     }
 			break;
-		case NodeSIdEnum:
-			doflist = xNew<int>(numnodes);
-			values = xNew<IssmDouble>(numnodes);
+		case NodeSIdEnum:{
+         int         numnodes = this->GetNumberOfNodes();
+         int        *doflist  = xNew<int>(numnodes);
+         IssmDouble *values   = xNew<IssmDouble>(numnodes);
 			/*Fill in values*/
 			this->GetNodesSidList(doflist);
 			this->GetInputListOnNodes(values,input_enum);
 			vector->SetValues(numnodes,doflist,values,INS_VAL);
+         xDelete<int>(doflist);
+         xDelete<IssmDouble>(values);
+                       }
 			break;
 		default:
@@ -1364,44 +1343,31 @@
 	}
 
-	/*Clean up*/
-	xDelete<int>(doflist);
-	xDelete<IssmDouble>(values);
-
-}
-/*}}}*/
+}/*}}}*/
 void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type, IssmDouble time){/*{{{*/
 
-	/*Fetch number vertices for this element and allocate arrays*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	int         numnodes    = this->GetNumberOfNodes();
-	int*        doflist     = NULL;
-	IssmDouble* values      = NULL;
-
 	switch(type){
-		case VertexPIdEnum:
-			doflist = xNew<int>(NUM_VERTICES);
-			values = xNew<IssmDouble>(NUM_VERTICES);
+		case VertexPIdEnum:{
+			int        doflist[MAXVERTICES];
+			IssmDouble values[MAXVERTICES];
+			const int  NUM_VERTICES = this->GetNumberOfVertices();
 			/*Fill in values*/
 			this->GetVerticesPidList(doflist);
 			this->GetInputListOnVerticesAtTime(values,input_enum,time);
-			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
+			vector->SetValues(NUM_VERTICES,doflist,&values[0],INS_VAL);
+								 }
 			break;
-		case VertexSIdEnum:
-			doflist = xNew<int>(NUM_VERTICES);
-			values = xNew<IssmDouble>(NUM_VERTICES);
+		case VertexSIdEnum:{
+			int        doflist[MAXVERTICES];
+			IssmDouble values[MAXVERTICES];
+			const int  NUM_VERTICES = this->GetNumberOfVertices();
 			/*Fill in values*/
 			this->GetVerticesSidList(doflist);
 			this->GetInputListOnVerticesAtTime(values,input_enum,time);
-			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
+			vector->SetValues(NUM_VERTICES,doflist,&values[0],INS_VAL);
+								 }
 			break;
 		default:
 			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
 	}
-
-	/*Clean up*/
-	xDelete<int>(doflist);
-	xDelete<IssmDouble>(values);
-
 }
 /*}}}*/
@@ -1449,11 +1415,9 @@
 	/*Create list of x*/
 	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	IssmDouble* x_list = xNew<IssmDouble>(NUM_VERTICES);
-
+	IssmDouble x_list[MAXVERTICES];
 	for(int i=0;i<NUM_VERTICES;i++) x_list[i]=xyz_list[i*3+0];
-	ValueP1OnGauss(&x,x_list,gauss);
-
-	xDelete<IssmDouble>(x_list);
+
+	/*Get value at gauss point*/
+	ValueP1OnGauss(&x,&x_list[0],gauss);
 	return x;
 }/*}}}*/
@@ -1465,11 +1429,9 @@
 	/*Create list of y*/
 	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	IssmDouble* y_list      = xNew<IssmDouble>(NUM_VERTICES);
-
+	IssmDouble y_list[MAXVERTICES];
 	for(int i=0;i<NUM_VERTICES;i++) y_list[i]=xyz_list[i*3+1];
-	ValueP1OnGauss(&y,y_list,gauss);
-
-	xDelete<IssmDouble>(y_list);
+
+	/*Get value at gauss point*/
+	ValueP1OnGauss(&y,&y_list[0],gauss);
 	return y;
 }/*}}}*/
@@ -1481,11 +1443,9 @@
 	/*Create list of z*/
 	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	IssmDouble* z_list      = xNew<IssmDouble>(NUM_VERTICES);
-
+	IssmDouble z_list[MAXVERTICES];
 	for(int i=0;i<NUM_VERTICES;i++) z_list[i]=xyz_list[i*3+2];
-	ValueP1OnGauss(&z,z_list,gauss);
-
-	xDelete<IssmDouble>(z_list);
+
+	/*Get value at gauss point*/
+	ValueP1OnGauss(&z,&z_list[0],gauss);
 	return z;
 }/*}}}*/
@@ -1612,7 +1572,4 @@
 void       Element::InputCreate(IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/
 
-	/*Intermediaries*/
-	int i,t;
-
 	/*Branch on type of vector: nodal or elementary: */
 	if(vector_type==1){ //nodal vector
@@ -1620,11 +1577,11 @@
 		const int NUM_VERTICES = this->GetNumberOfVertices();
 
-		int        *vertexids  = xNew<int>(NUM_VERTICES);
-		int        *vertexlids = xNew<int>(NUM_VERTICES);
-		IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);
+		int        vertexids[MAXVERTICES];
+		int        vertexlids[MAXVERTICES];
+		IssmDouble values[MAXVERTICES];
 
 		/*Recover vertices ids needed to initialize inputs*/
 		_assert_(iomodel->elements);
-		for(i=0;i<NUM_VERTICES;i++){
+		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
 			vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
@@ -1637,5 +1594,5 @@
 		}
 		else if(M==iomodel->numberofvertices){
-			for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
+			for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
 			this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
 		}
@@ -1643,12 +1600,12 @@
 			/*create transient input: */
 			IssmDouble* times = xNew<IssmDouble>(N);
-			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
+			for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
 			inputs->SetTransientInput(vector_enum,times,N);
 			TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
-			for(t=0;t<N;t++){
-				for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
+			for(int t=0;t<N;t++){
+				for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
 				switch(this->ObjectEnum()){
-					case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
-					case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,values,P1Enum); break;
+					case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,&values[0],P1Enum); break;
+					case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,&values[0],P1Enum); break;
 					default: _error_("Not implemented yet");
 				}
@@ -1659,23 +1616,23 @@
 
 			/*This is a Patch!*/
-			xDelete<IssmDouble>(values);
-			values = xNew<IssmDouble>(N);
-			for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
+			IssmDouble* evalues = xNew<IssmDouble>(N);
+			for(int j=0;j<N;j++) evalues[j]=vector[this->Sid()*N+j];
 
 			if (N==this->GetNumberOfNodes(P1Enum)){
-				this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
+				this->SetElementInput(inputs,NUM_VERTICES,vertexlids,evalues,vector_enum);
 			}
 			else if(N==this->GetNumberOfNodes(P0Enum)){
-				this->SetElementInput(inputs,vector_enum,values[0]);
+				this->SetElementInput(inputs,vector_enum,evalues[0]);
 			}
 			else if(N==this->GetNumberOfNodes(P1xP2Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
-				inputs->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,values);
+				inputs->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,evalues);
 			}
 			else if(N==this->GetNumberOfNodes(P1xP3Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
-				inputs->SetPentaInput(vector_enum,P1xP3Enum,this->lid,N,values);
+				inputs->SetPentaInput(vector_enum,P1xP3Enum,this->lid,N,evalues);
 			}
 			else{
 				_error_("Patch interpolation not supported yet");
 			}
+			xDelete<IssmDouble>(evalues);
 
 		}
@@ -1683,8 +1640,4 @@
 			_error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
 		}
-
-		xDelete<IssmDouble>(values);
-		xDelete<int>(vertexids);
-		xDelete<int>(vertexlids);
 	}
 	else if(vector_type==2){ //element vector
@@ -1708,8 +1661,8 @@
 			/*create transient input: */
 			IssmDouble* times = xNew<IssmDouble>(N);
-			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
+			for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
 			inputs->SetTransientInput(vector_enum,times,N);
 			TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
-			for(t=0;t<N;t++){
+			for(int t=0;t<N;t++){
 				value=vector[N*this->Sid()+t];
 				switch(this->ObjectEnum()){
@@ -1728,5 +1681,5 @@
 		if(M==iomodel->numberofelements){
 			IssmDouble* layers = xNewZeroInit<IssmDouble>(N);
-			for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
+			for(int t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
 			inputs->SetArrayInput(vector_enum,this->lid,layers,N);
 			xDelete<IssmDouble>(layers);
@@ -1850,7 +1803,4 @@
 	 */
 
-	/*Intermediaries*/
-	int i,t;
-
 	/*Branch on type of vector: nodal or elementary: */
 	if(vector_type==1){ //nodal vector
@@ -1858,11 +1808,11 @@
 		const int NUM_VERTICES = this->GetNumberOfVertices();
 
-		int        *vertexids  = xNew<int>(NUM_VERTICES);
-		int        *vertexlids = xNew<int>(NUM_VERTICES);
-		IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);
+		int        vertexids[MAXVERTICES];
+		int        vertexlids[MAXVERTICES];
+		IssmDouble values[MAXVERTICES];
 
 		/*Recover vertices ids needed to initialize inputs*/
 		_assert_(iomodel->elements);
-		for(i=0;i<NUM_VERTICES;i++){
+		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
 			vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
@@ -1876,5 +1826,5 @@
 		}
 		else if(M==iomodel->numberofvertices){
-			for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
+			for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
 			switch(this->ObjectEnum()){
 				case TriaEnum:  inputs->SetTriaDatasetInput(enum_type,input_id,P1Enum,NUM_VERTICES,vertexlids,values); break;
@@ -1886,8 +1836,8 @@
 			/*create transient input: */
 			IssmDouble* times = xNew<IssmDouble>(N);
-			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
+			for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
 			TransientInput* transientinput = inputs->SetDatasetTransientInput(enum_type,input_id,times,N);
-			for(t=0;t<N;t++){
-				for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
+			for(int t=0;t<N;t++){
+				for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
 				switch(this->ObjectEnum()){
 					case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
@@ -1901,8 +1851,4 @@
 			_error_("not implemented yet (M="<<M<<")");
 		}
-
-		xDelete<IssmDouble>(values);
-		xDelete<int>(vertexids);
-		xDelete<int>(vertexlids);
 	}
 	else if(vector_type==2){ //element vector
@@ -2065,5 +2011,5 @@
 	/*Allocate some arrays*/
 	const int numvertices = this->GetNumberOfVertices();
-	IssmDouble* basalmeltrate = xNew<IssmDouble>(numvertices);
+	IssmDouble basalmeltrate[MAXVERTICES];
 
 	/*Get variables*/
@@ -2114,5 +2060,4 @@
 	xDelete<IssmDouble>(mean_tf);
 	xDelete<IssmDouble>(depths);
-	xDelete<IssmDouble>(basalmeltrate);
 
 }/*}}}*/
@@ -2121,9 +2066,9 @@
 	const int NUM_VERTICES = this->GetNumberOfVertices();
 
-	IssmDouble  deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
-	IssmDouble *base         = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble *values       = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble *perturbation = xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble 	time;
+	IssmDouble deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
+	IssmDouble base[MAXVERTICES];
+	IssmDouble values[MAXVERTICES];
+	IssmDouble perturbation[MAXVERTICES];
+	IssmDouble time;
 
 	parameters->FindParam(&time,TimeEnum);
@@ -2134,6 +2079,6 @@
 	_assert_(upperwaterel>deepwaterel);
 
-	this->GetInputListOnVertices(base,BaseEnum);
-   this->GetInputListOnVertices(perturbation,BasalforcingsPerturbationMeltingRateEnum);
+	this->GetInputListOnVertices(&base[0],BaseEnum);
+   this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum);
 	for(int i=0;i<NUM_VERTICES;i++){
 		if(base[i]>=upperwaterel){
@@ -2152,7 +2097,4 @@
 
 	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
-	xDelete<IssmDouble>(base);
-   xDelete<IssmDouble>(perturbation);
-	xDelete<IssmDouble>(values);
 
 }/*}}}*/
@@ -2161,14 +2103,14 @@
 	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);
-	this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
-	this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
-	this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
+	IssmDouble deepwatermelt[MAXVERTICES];
+	IssmDouble deepwaterel[MAXVERTICES];;
+	IssmDouble upperwaterel[MAXVERTICES];
+	IssmDouble base[MAXVERTICES];
+	IssmDouble values[MAXVERTICES];
+
+	this->GetInputListOnVertices(&base[0],BaseEnum);
+	this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsDeepwaterMeltingRateEnum);
+	this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsDeepwaterElevationEnum);
+	this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsUpperwaterElevationEnum);
 
 	for(int i=0;i<NUM_VERTICES;i++){
@@ -2178,11 +2120,5 @@
 	}
 
-	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
-	xDelete<IssmDouble>(base);
-	xDelete<IssmDouble>(deepwaterel);
-	xDelete<IssmDouble>(deepwatermelt);
-	xDelete<IssmDouble>(upperwaterel);
-	xDelete<IssmDouble>(values);
-
+	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&values[0],P1Enum);
 }/*}}}*/
 void       Element::MantlePlumeGeothermalFlux(){/*{{{*/
@@ -2193,5 +2129,5 @@
 	IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
 	IssmDouble  x,y,z,c;
-	IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble  values[MAXVERTICES];
 	IssmDouble *xyz_list = NULL;
 
@@ -2228,7 +2164,6 @@
 	}
 
-	this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum);
+	this->AddInput(BasalforcingsGeothermalfluxEnum,&values[0],P1Enum);
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(values);
 
 }/*}}}*/
@@ -2253,14 +2188,13 @@
 
 	const int  NUM_VERTICES = this->GetNumberOfVertices();
-	int        i,migration_style;
+	int        migration_style;
 	IssmDouble bed_hydro,yts;
-	IssmDouble rho_water,rho_ice,density;
-	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);
+	IssmDouble melting[MAXVERTICES];
+	IssmDouble phi[MAXVERTICES];
+	IssmDouble h[MAXVERTICES];
+	IssmDouble s[MAXVERTICES];
+	IssmDouble b[MAXVERTICES];
+	IssmDouble r[MAXVERTICES];
+	IssmDouble sl[MAXVERTICES];
 
 	/*Recover info at the vertices: */
@@ -2273,10 +2207,10 @@
 	GetInputListOnVertices(&sl[0],SealevelEnum);
 	GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
-	rho_water   = FindParam(MaterialsRhoSeawaterEnum);
-	rho_ice     = FindParam(MaterialsRhoIceEnum);
-	density     = rho_ice/rho_water;
+	IssmDouble rho_water   = FindParam(MaterialsRhoSeawaterEnum);
+	IssmDouble rho_ice     = FindParam(MaterialsRhoIceEnum);
+	IssmDouble density     = rho_ice/rho_water;
 
 	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
-	for(i=0;i<NUM_VERTICES;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		/* Contact FS*/
 		if(migration_style == ContactEnum){
@@ -2284,6 +2218,8 @@
 			if(phi[i]>=0.) b[i]=r[i];
 		}
-		else if(migration_style == GroundingOnlyEnum && b[i]<r[i]) b[i]=r[i];
-		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
+		else if(migration_style == GroundingOnlyEnum && b[i]<r[i]){
+			/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
+			b[i]=r[i];
+		}
 		else if(phi[i]<=0.){
 			if(b[i]<=r[i]){
@@ -2299,10 +2235,10 @@
 				/*Unground only if the element is connected to the ice shelf*/
 				if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum){
-					s[i]        = (1-density)*h[i]+sl[i];
-					b[i]        = -density*h[i]+sl[i];
+					s[i] = (1.-density)*h[i]+sl[i];
+					b[i] = -density*h[i]+sl[i];
 				}
 				else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
-					s[i]        = (1-density)*h[i]+sl[i];
-					b[i]        = -density*h[i]+sl[i];
+					s[i] = (1.-density)*h[i]+sl[i];
+					b[i] = -density*h[i]+sl[i];
 				}
 				else{
@@ -2314,5 +2250,5 @@
 
 	/*Recalculate phi*/
-	for(i=0;i<NUM_VERTICES;i++){
+	for(int i=0;i<NUM_VERTICES;i++){
 		if(migration_style==SoftMigrationEnum){
 			bed_hydro=-density*h[i]+sl[i];
@@ -2321,5 +2257,7 @@
 			}
 		}
-		else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density;
+		else if(migration_style!=ContactEnum){
+			phi[i]=h[i]+(r[i]-sl[i])/density;
+		}
 		else{
 			/*do nothing*/
@@ -2332,15 +2270,5 @@
 	this->AddInput(BaseEnum,&b[0],P1Enum);
 
-	/*Delete*/
-	xDelete<IssmDouble>(melting);
-	xDelete<IssmDouble>(phi);
-	xDelete<IssmDouble>(r);
-	xDelete<IssmDouble>(b);
-	xDelete<IssmDouble>(s);
-	xDelete<IssmDouble>(sl);
-	xDelete<IssmDouble>(h);
-
-}
-/*}}}*/
+}/*}}}*/
 void       Element::MismipFloatingiceMeltingRate(){/*{{{*/
 
@@ -2566,5 +2494,4 @@
 
 	const int NUM_VERTICES = this->GetNumberOfVertices();
-
 	int        basinid, maxbox, num_basins, numnodes, M;
 	IssmDouble gamma_T, overturning_coeff, thickness;
@@ -2593,6 +2520,7 @@
 	this->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
 	this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
-	Input *thickness_input    = this->GetInput(ThicknessEnum);                         _assert_(thickness_input);
-	_assert_(basinid<=num_basins);
+	Input *thickness_input = this->GetInput(ThicknessEnum); 
+   _assert_(basinid<=num_basins);
+   _assert_(thickness_input);
 
 	IssmDouble* boxareas = NULL;
@@ -2603,8 +2531,8 @@
 	IssmDouble g1               = area_boxi*gamma_T;
 
-	IssmDouble* basalmeltrates_shelf 				= xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble* potential_pressure_melting_point	= xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble* Tocs 								= xNew<IssmDouble>(NUM_VERTICES);
-	IssmDouble* Socs								= xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble basalmeltrates_shelf[MAXVERTICES];
+	IssmDouble potential_pressure_melting_point[MAXVERTICES];
+	IssmDouble Tocs[MAXVERTICES];
+	IssmDouble Socs[MAXVERTICES];
 
 	/* First box calculations */
@@ -2612,9 +2540,9 @@
 		/* Get box1 parameters and inputs */
 		IssmDouble time, toc_farocean, soc_farocean;
+      IssmDouble overturnings[MAXVERTICES];
 		this->parameters->FindParam(&time,TimeEnum);
 		this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
 		this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
-		IssmDouble 	s1 				= soc_farocean/(nu*lambda);
-		IssmDouble* overturnings 	= xNew<IssmDouble>(NUM_VERTICES);
+		IssmDouble s1 = soc_farocean/(nu*lambda);
 		Input *overturningC_input = this->GetInput(BasalforcingsPicoOverturningCoeffEnum); _assert_(overturningC_input);
 
@@ -2640,8 +2568,8 @@
 		}
 
-		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1DGEnum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1DGEnum);
+		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&basalmeltrates_shelf[0],P1DGEnum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,&Tocs[0],P1DGEnum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,&Socs[0],P1DGEnum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,&overturnings[0],P1DGEnum);
 
 		/*Cleanup and return*/
@@ -2716,5 +2644,5 @@
 
 	/*Define arrays*/
-	IssmDouble* basalmeltrates_shelf = xNew<IssmDouble>(NUM_VERTICES); //Basal melt-rate
+	IssmDouble basalmeltrates_shelf[MAXVERTICES];
 
 	/*Polynomial coefficients*/
@@ -2789,9 +2717,8 @@
 
 	/*Save computed melt-rate*/
-	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
+	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&basalmeltrates_shelf[0],P1DGEnum);
 
 	/*Cleanup and return*/
 	delete gauss;
-
 }/*}}}*/
 void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
@@ -4232,6 +4159,5 @@
 	/*Fetch number vertices and allocate memory*/
 	const int NUM_VERTICES  = this->GetNumberOfVertices();
-
-	IssmDouble* maxprincipal = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble maxprincipal[MAXVERTICES];
 
 	/*Retrieve all inputs and parameters*/
@@ -4251,5 +4177,5 @@
 	/*loop over vertices: */
 	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<NUM_VERTICES;iv++){
+	for(int iv=0;iv<NUM_VERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -4296,8 +4222,7 @@
 
 	/*Create input*/
-	this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
+	this->AddInput(StressMaxPrincipalEnum,&maxprincipal[0],P1Enum);
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(maxprincipal);
 	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
@@ -4575,6 +4500,5 @@
 	/*Fetch number vertices and allocate memory*/
 	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	IssmDouble* viscousheating = xNew<IssmDouble>(NUM_VERTICES);
+	IssmDouble viscousheating[MAXVERTICES];
 
 	/*Retrieve all inputs and parameters*/
@@ -4595,8 +4519,7 @@
 
 	/*Create PentaVertex input, which will hold the basal friction:*/
-	this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
+	this->AddInput(ViscousHeatingEnum,&viscousheating[0],P1Enum);
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(viscousheating);
 	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
@@ -4683,6 +4606,6 @@
 	IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
 	IssmDouble  Hc,Ht;
-	IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
-	IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
+	IssmDouble PIE[MAXVERTICES];
+	IssmDouble dHpmp[MAXVERTICES];
 
 	for(int iv=0; iv<numvertices; iv++){
@@ -4720,6 +4643,4 @@
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(PIE);
-	xDelete<IssmDouble>(dHpmp);
 	return kappa;
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25904)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25905)
@@ -2222,5 +2222,5 @@
 
 	/*Intermediaries*/
-	int i, numiceverts, numnoiceverts;
+	int numiceverts, numnoiceverts;
 	int ind0, ind1, lastindex;
 	int indices_ice[NUMVERTICES],indices_noice[NUMVERTICES];
@@ -2234,5 +2234,5 @@
 	 * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
 	lastindex=0;
-	for(i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
+	for(int i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
 		ind0=(NUMVERTICES-i)%NUMVERTICES;
 		ind1=(ind0-1+NUMVERTICES)%NUMVERTICES;
@@ -2248,5 +2248,5 @@
 	numiceverts=0;
 	numnoiceverts=0;
-	for(i=0;i<NUMVERTICES;i++){
+	for(int i=0;i<NUMVERTICES;i++){
 		ind0=(lastindex+i)%NUMVERTICES;
 		if(lsf[i]<=level){
@@ -2260,24 +2260,28 @@
 	}
 	//merge indices
-	for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
-	for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
-
-	switch (numiceverts){
+	for(int i=0;i<numiceverts;  i++){
+		indices[i]=indices_ice[i];
+	}
+	for(int i=0;i<numnoiceverts;i++){
+		indices[numiceverts+i]=indices_noice[i];
+	}
+
+	switch(numiceverts){
 		case 0: // no vertex has ice: element is ice free, no intersection
-			for(i=0;i<2;i++)
+			for(int i=0;i<2;i++)
 				fraction[i]=0.;
 			break;
 		case 1: // one vertex has ice:
-			for(i=0;i<2;i++){
+			for(int i=0;i<2;i++){
 				fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
 			}
 			break;
 		case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
-			for(i=0;i<2;i++){
+			for(int i=0;i<2;i++){
 				fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
 			}
 			break;
 		case NUMVERTICES: // all vertices have ice: return triangle area
-			for(i=0;i<2;i++)
+			for(int i=0;i<2;i++)
 				fraction[i]=1.;
 			break;
@@ -3257,18 +3261,10 @@
 
 	/*intermediary: */
-	IssmDouble* values=NULL;
-	Input*     thickness_input=NULL;
 	IssmDouble  thickness;
-	IssmDouble  weight;
 	IssmDouble  Jdet;
-	IssmDouble  volume;
-	IssmDouble  rho_ice;
 	int         point1;
 	IssmDouble  fraction1,fraction2;
 	bool        mainlynegative=true;
 
-	/*Output:*/
-	volume=0;
-
 	/* Get node coordinates and dof list: */
 	IssmDouble  xyz_list[NUMVERTICES][3];
@@ -3276,11 +3272,11 @@
 
 	/*Retrieve inputs required:*/
-	thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
 
 	/*Retrieve material parameters: */
-	rho_ice=FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
 
 	/*Retrieve values of the levelset defining the masscon: */
-	values = xNew<IssmDouble>(NUMVERTICES);
+	IssmDouble values[NUMVERTICES];
 	for(int i=0;i<NUMVERTICES;i++){
 		values[i]=levelset[this->vertices[i]->Sid()];
@@ -3288,19 +3284,15 @@
 
 	/*Ok, use the level set values to figure out where we put our gaussian points:*/
-	this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
+	this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,&values[0]);
 	Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
 
-	volume=0;
-
+	IssmDouble volume=0.;
 	while(gauss->next()){
-
 		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
 		thickness_input->GetInputValue(&thickness, gauss);
-
 		volume+=thickness*gauss->weight*Jdet;
 	}
 
 	/* clean up and Return: */
-	xDelete<IssmDouble>(values);
 	delete gauss;
 	return rho_ice*volume;