Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 22728)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 22729)
@@ -6,8 +6,7 @@
 #include "../../shared/shared.h"
 
-void InterpFromMeshToGridx(double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value){
+void InterpFromMeshToGridx(double** pgriddata,int* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value){
 
 	/*Intermediary*/
-	int    i,j;
 	int    i1,i2,j1,j2;
 	int    interpolation_type;
@@ -28,18 +27,16 @@
 	else _error_("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
 
-	/*First, allocate pointers: */
-	double* griddata=xNewZeroInit<double>(nlines*ncols);
+	/*First, prepare output*/
+	double* griddata=xNew<double>(nlines*ncols);
+	for(int i=0;i<nlines;i++) for(int j=0;j<ncols; j++) griddata[i*ncols+j]=default_value;
 
-	/*Set debug to 1 if there are lots of elements*/
+	/*Set debug to "true" if there are lots of elements*/
 	bool debug=(bool)((double)ncols*nlines*nels >= 5.e10);
-
-	/*Initialize coordintes and griddata*/
-	for(i=0;i<nlines;i++) for(j=0;j<ncols; j++) griddata[i*ncols+j]=default_value;
 
 	/*figure out if x or y are flipped*/
 	bool xflip = false;
 	bool yflip = false;
-	if (x_grid[1]-x_grid[0]<0) xflip=true;
-	if (y_grid[1]-y_grid[0]<0) yflip=true;
+	if(x_grid[1]-x_grid[0]<0) xflip=true;
+	if(y_grid[1]-y_grid[0]<0) yflip=true;
 
 	/*Get min/max coordinates of the grid*/
@@ -67,15 +64,15 @@
 
 		/*display current iteration*/
-		if(debug && fmod((double)n,(double)100)==0)
+		if(debug && fmod(double(n),1000)==0)
 		 _printf_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(n)/double(nels)*100<<"%   ");
 
 		/*Get extrema coordinates of current elements*/
-		x_tria_min=x_mesh[(int)index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
-		y_tria_min=y_mesh[(int)index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
-		for(i=1;i<3;i++){
-			if(x_mesh[(int)index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[(int)index_mesh[3*n+i]-1];
-			if(x_mesh[(int)index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[(int)index_mesh[3*n+i]-1];
-			if(y_mesh[(int)index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[(int)index_mesh[3*n+i]-1];
-			if(y_mesh[(int)index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[(int)index_mesh[3*n+i]-1];
+		x_tria_min=x_mesh[index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
+		y_tria_min=y_mesh[index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
+		for(int i=1;i<3;i++){
+			if(x_mesh[index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[index_mesh[3*n+i]-1];
+			if(x_mesh[index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[index_mesh[3*n+i]-1];
+			if(y_mesh[index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[index_mesh[3*n+i]-1];
+			if(y_mesh[index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[index_mesh[3*n+i]-1];
 		}
 
@@ -83,5 +80,5 @@
 		if( (x_tria_min>x_grid_max) || (x_tria_max<x_grid_min) || (y_tria_min>y_grid_max) || (y_tria_max<y_grid_min) ) continue;
 
-		/*Get indices i and j that form a square around the currant triangle*/
+		/*Get indices i and j that form a square around the current triangle*/
 		if(yflip){
 			i1=max(0,       (int)floor((y_tria_max-y_grid_max)/yposting)-1);
@@ -103,15 +100,15 @@
 		/*get area of the current element (Jacobian = 2 * area)*/
 		//area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
-		area=x_mesh[(int)index_mesh[3*n+1]-1]*y_mesh[(int)index_mesh[3*n+2]-1]-y_mesh[(int)index_mesh[3*n+1]-1]*x_mesh[(int)index_mesh[3*n+2]-1]
-		  +  x_mesh[(int)index_mesh[3*n+0]-1]*y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+0]-1]*x_mesh[(int)index_mesh[3*n+1]-1]
-		  +  x_mesh[(int)index_mesh[3*n+2]-1]*y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1]*x_mesh[(int)index_mesh[3*n+0]-1];
+		area=x_mesh[index_mesh[3*n+1]-1]*y_mesh[index_mesh[3*n+2]-1]-y_mesh[index_mesh[3*n+1]-1]*x_mesh[index_mesh[3*n+2]-1]
+		  +  x_mesh[index_mesh[3*n+0]-1]*y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+0]-1]*x_mesh[index_mesh[3*n+1]-1]
+		  +  x_mesh[index_mesh[3*n+2]-1]*y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1]*x_mesh[index_mesh[3*n+0]-1];
 
 		/*Go through x_grid and y_grid and interpolate if necessary*/
-		for(i=i1;i<=i2;i++){
+		for(int i=i1;i<=i2;i++){
 
 			//exit if y not between y_tria_min and y_tria_max
 			if((y_grid[i]>y_tria_max) || (y_grid[i]<y_tria_min)) continue;
 
-			for(j=j1;j<=j2; j++){
+			for(int j=j1;j<=j2; j++){
 
 				//exit if x not between x_tria_min and x_tria_max
@@ -119,19 +116,19 @@
 
 				/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
-				area_1=((x_grid[j]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+2]-1]) 
-							-  (y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_mesh[(int)index_mesh[3*n+1]-1]-x_mesh[(int)index_mesh[3*n+2]-1]))/area;
+				area_1=((x_grid[j]-x_mesh[index_mesh[3*n+2]-1])*(y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+2]-1]) 
+							-  (y_grid[i]-y_mesh[index_mesh[3*n+2]-1])*(x_mesh[index_mesh[3*n+1]-1]-x_mesh[index_mesh[3*n+2]-1]))/area;
 				/*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
-				area_2=((x_mesh[(int)index_mesh[3*n+0]-1]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1]) 
-							- (y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[(int)index_mesh[3*n+2]-1]))/area;
+				area_2=((x_mesh[index_mesh[3*n+0]-1]-x_mesh[index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[index_mesh[3*n+2]-1]) 
+							- (y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[index_mesh[3*n+2]-1]))/area;
 				/*Get third area coordinate = 1-area1-area2*/
 				area_3=1.-area_1-area_2;
 
 				/*is the current point in the current element?*/
-				if(area_1>-10e-12 && area_2>-10e-12 && area_3>-10e-12){
+				if(area_1>-1e-12 && area_2>-1e-12 && area_3>-1e-12){
 
 					/*Yes ! compute the value on the point*/
 					if(interpolation_type==1){
 						/*nodal interpolation*/
-						data_value=area_1*data_mesh[(int)index_mesh[3*n+0]-1]+area_2*data_mesh[(int)index_mesh[3*n+1]-1]+area_3*data_mesh[(int)index_mesh[3*n+2]-1];
+						data_value=area_1*data_mesh[index_mesh[3*n+0]-1]+area_2*data_mesh[index_mesh[3*n+1]-1]+area_3*data_mesh[index_mesh[3*n+2]-1];
 					}
 					else{
@@ -147,5 +144,5 @@
 		}
 	}
-	if(debug) _printf_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%  \n");
+	if(debug) _printf_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<fixed<<100.<<"%   \n");
 
 	/*Assign output pointers:*/
Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h	(revision 22728)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h	(revision 22729)
@@ -8,5 +8,5 @@
 #include "../../toolkits/toolkits.h"
 
-void InterpFromMeshToGridx(double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value);
+void InterpFromMeshToGridx(double** pgriddata,int* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value);
 
 #endif /* _INTERPFROMMESHTOGRIDX_H */
Index: /issm/trunk-jpl/src/wrappers/InterpFromMeshToGrid/InterpFromMeshToGrid.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/InterpFromMeshToGrid/InterpFromMeshToGrid.cpp	(revision 22728)
+++ /issm/trunk-jpl/src/wrappers/InterpFromMeshToGrid/InterpFromMeshToGrid.cpp	(revision 22729)
@@ -23,5 +23,5 @@
 
 	/*inputs */
-	double* index=NULL;
+	int*    index=NULL;
 	double* x=NULL;
 	double* y=NULL;
@@ -63,5 +63,5 @@
 
 	/*Free ressources: */
-	xDelete<double>(index);
+	xDelete<int>(index);
 	xDelete<double>(x);
 	xDelete<double>(y);