Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 7446)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 7447)
@@ -8,5 +8,4 @@
 
 enum definitions{
-
 	/*Datasets {{{1*/
 	ConstraintsEnum,
@@ -411,5 +410,10 @@
 	SsetEnum,
 	GroundingLineMigrationEnum,
-	YtsEnum
+	YtsEnum,
+	/*}}}*/
+	/*Interpolation {{{1*/
+	TriangleInterpEnum,
+	BilinearInterpEnum,
+	NearestInterpEnum
 	/*}}}*/
 };
Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 7446)
+++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 7447)
@@ -366,4 +366,7 @@
 		case GroundingLineMigrationEnum : return "GroundingLineMigration";
 		case YtsEnum : return "Yts";
+		case TriangleInterpEnum : return "TriangleInterp";
+		case BilinearInterpEnum : return "BilinearInterp";
+		case NearestInterpEnum : return "NearestInterp";
 		default : return "unknown";
 
Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 7446)
+++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 7447)
@@ -364,4 +364,7 @@
 	else if (strcmp(name,"GroundingLineMigration")==0) return GroundingLineMigrationEnum;
 	else if (strcmp(name,"Yts")==0) return YtsEnum;
+	else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
+	else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
+	else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
 	else _error_("Enum %s not found",name);
 
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 7446)
+++ /issm/trunk/src/c/Makefile.am	(revision 7447)
@@ -505,5 +505,4 @@
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
-					./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
 					./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
 					./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\
@@ -1084,5 +1083,4 @@
 					./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\
 					./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
-					./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
 					./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
 					./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\
Index: /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 7446)
+++ /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 7447)
@@ -3,4 +3,5 @@
  */ 
 
+/*Include {{{*/
 #ifdef HAVE_CONFIG_H
 	#include "config.h"
@@ -12,7 +13,8 @@
 #include "../../shared/shared.h"
 #include "../../include/include.h"
-
-int InterpFromGridToMeshx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value) {
-
+/*}}}*/
+
+/*InterpFromGridToMeshx{{{*/
+int InterpFromGridToMeshx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value, int interpenum){
 
 	/*output: */
@@ -28,5 +30,4 @@
 	InterpFromGridToMeshxThreadStruct gate;
 	int num=1;
-
 	#ifdef _MULTITHREADING_
 	num=_NUMTHREADS_;
@@ -82,4 +83,5 @@
 	gate.M=M;
 	gate.N=N;
+	gate.interp=interpenum;
 
 	/*launch the thread manager with InterpFromGridToMeshxt as a core: */
@@ -89,2 +91,225 @@
 	*pdata_mesh=data_mesh;
 }
+/*}}}*/
+/*InterpFromGridToMeshxt {{{1*/
+void* InterpFromGridToMeshxt(void* vpthread_handle){
+
+	/*gate variables :*/
+	InterpFromGridToMeshxThreadStruct *gate    = NULL;
+	pthread_handle                    *handle  = NULL;
+	int my_thread;
+	int num_threads;
+	int i0,i1;
+
+	/*intermediary: */
+	int    i,m,n;
+	double x_grid;
+	double y_grid;
+	double data_value;
+	double x1,x2,y1,y2;
+	double Q11,Q12,Q21,Q22;
+
+	/*recover handle and gate: */
+	handle=(pthread_handle*)vpthread_handle;
+	gate=(InterpFromGridToMeshxThreadStruct*)handle->gate;
+	my_thread=handle->id;
+	num_threads=handle->num;
+
+	/*recover parameters :*/
+	double *x_mesh        = gate->x_mesh;
+	double *y_mesh        = gate->y_mesh;
+	int     x_rows        = gate->x_rows;
+	int     y_rows        = gate->y_rows;
+	double *x             = gate->x;
+	double *y             = gate->y;
+	int     nods          = gate->nods;
+	Vec     data_mesh     = gate->data_mesh;
+	double *data          = gate->data;
+	double  default_value = gate->default_value;
+	int     interpenum    = gate->interp;
+	int     M             = gate->M;
+	int     N             = gate->N;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,nods,num_threads,my_thread);
+	for (i=i0;i<i1;i++) {
+
+		x_grid=*(x_mesh+i);
+		y_grid=*(y_mesh+i);
+
+		/*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
+		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
+
+			/*    Q12             Q22
+			 * y2 x---------+-----x
+			 *    |         |     |
+			 *    |         |P    |
+			 *    |---------+-----|
+			 *    |         |     |
+			 *    |         |     |
+			 * y1 x---------+-----x Q21
+			 *    x1                 x2       
+			 *
+			 */
+			x1=x[n]; x2=x[n+1];
+			y1=y[m]; y2=y[m+1];
+			Q11=data[m*N+n];
+			Q12=data[(m+1)*N+n];
+			Q21=data[m*N+n+1];
+			Q22=data[(m+1)*N+n+1];
+
+			switch(interpenum){
+				case TriangleInterpEnum:
+					data_value=triangleinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
+					break;
+				case BilinearInterpEnum:
+					data_value=bilinearinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
+					break;
+				case NearestInterpEnum:
+					data_value=nearestinterp(x1,x2,y1,y2, Q11,Q12,Q21,Q22,x_grid,y_grid);
+					break;
+				default:
+					printf("Interpolation %s not supported yet\n",EnumToString(interpenum));
+					return NULL; /*WARNING: no error because it would blow up the multithreading!*/
+			}
+			if(isnan(data_value)) data_value=default_value;
+		}
+		else{
+			data_value=default_value;
+		}
+
+		VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
+	}
+}/*}}}*/
+
+/*findindices {{{1*/
+bool findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){
+
+	bool foundx=false,foundy=false;
+	int m=-1,n=-1;
+	int i;
+
+	for (i=0;i<x_rows-1;i++){
+		if ((x[i]<=xgrid) && (xgrid<x[i+1])){
+			n=i;
+			foundx=true;
+			break;
+		}
+	}
+	if(xgrid==x[x_rows-1]){
+		n=x_rows-2;
+		foundx=true;
+	}
+
+	for (i=0;i<y_rows-1;i++){
+		if ((y[i]<=ygrid) && (ygrid<y[i+1])){
+			m=i;
+			foundy=true;
+			break;
+		}
+	}
+	if(ygrid==y[y_rows-1]){
+		m=y_rows-2;
+		foundy=true;
+	}
+
+	/*Assign output pointers:*/
+	*pm=m; *pn=n;
+	return (foundx && foundy);
+}/*}}}*/
+/*triangleinterp{{{1*/
+double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+	/*split the rectangle in 2 triangle and
+	 * use Lagrange P1 interpolation
+	 *       
+	 *   +3----+2,3' Q12----Q22
+	 *   |    /|     |    /|
+	 *   |   / |     |   / |
+	 *   |  /  |     |  /  |
+	 *   | /   |     | /   |
+	 *   |/    |     |/    |
+	 *   1-----2'    Q11---Q21        */
+
+	/*Intermediaries*/
+	double area,area_1,area_2,area_3;
+
+	/*Checks*/
+	_assert_(x2>x1 && y2>y1);
+	_assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+
+	/*area of the rectangle*/
+	area=(x2-x1)*(y2-y1);
+
+	/*is it the upper left triangle?*/
+	if ((x-x1)/(x2-x1)<(y-y1)/(y2-y1)){
+
+		area_1=((y2-y)*(x2-x1))/area;
+		area_2=((x-x1)*(y2-y1))/area;
+		area_3=1-area_1-area_2;
+
+		return area_1*Q11+area_2*Q22+area_3*Q12;
+	}
+	else {
+
+		area_1=((y-y1)*(x2-x1))/area;
+		area_2=((x2-x)*(y2-y1))/area;
+		area_3=1-area_1-area_2;
+
+		return area_1*Q22+area_2*Q11+area_3*Q21;
+	}
+}/*}}}*/
+/*bilinearinterp{{{1*/
+double bilinearinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+	/*Bilinear  interpolation: (http://en.wikipedia.org/wiki/Bilinear_interpolation) */
+
+	/*    Q12    R2        Q22
+	 * y2 x------x---------x
+	 *    |      |         |
+	 *    |      |         |
+	 *    |      +P        |
+	 *    |      |         |
+	 *    |Q11   R1        Q21
+	 * y1 x------x---------x
+	 *    x1               x2
+	 *
+	 */
+
+	/*Checks*/
+	_assert_(x2>x1 && y2>y1);
+	_assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+
+	return 
+	  +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
+	  +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
+	  +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1))
+	  +Q22*(x-x1)*(y-y1)/((x2-x1)*(y2-y1));
+}
+/*}}}*/
+/*nearestinterp{{{1*/
+double nearestinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+	/*Nearest neighbor interpolation*/
+
+	/*    Q12             Q22
+	 * y2 x--------x---------x
+	 *    |        |         |
+	 *    |        |  xP     |
+	 * ym |--------+---------|
+	 *    |        |         |
+	 *    |        |         |
+	 * y1 x--------x---------x Q21
+	 *    x1       xm        x2 
+	 *
+	 */
+	/*Checks*/
+	_assert_(x2>x1 && y2>y1);
+	_assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+	
+	double xm=(x2-x1)/2;
+	double ym=(y2-y1)/2;
+
+	if (x<xm && y<ym) return Q11;
+	if (x<xm && y>ym) return Q12;
+	if (x>xm && y<ym) return Q21;
+	else return Q22;
+}
+/*}}}*/
Index: /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h
===================================================================
--- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h	(revision 7446)
+++ /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h	(revision 7447)
@@ -7,32 +7,29 @@
 
 #include "../../toolkits/toolkits.h"
-
-
+#include "../../EnumDefinitions/EnumDefinitions.h"
 
 /*threading: */
 typedef struct{
-
+	double* x;
+	int     x_rows;
+	double* y;
+	int     y_rows;
+	double* data;
+	double  default_value;
+	int     interp;
+	int     M;
+	int     N;
+	int     nods;
 	double* x_mesh;
 	double* y_mesh;
-	int     x_rows;
-	int     y_rows;
-	double* x;
-	double* y;
-	int     nods;
-	double  default_value;
 	Vec     data_mesh;
-	double* data;
-	int     M;
-	int     N;
-
 } InterpFromGridToMeshxThreadStruct;
 
-int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
-
-int InterpFromGridToMeshx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value);
-
-void* InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
-
+int    InterpFromGridToMeshx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value, int interpenum=BilinearInterpEnum);
+void*  InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
+bool   findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
+double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
+double bilinearinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
+double nearestinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
 
 #endif /* _INTERPFROMGRIDTOMESHX_H */
-
Index: sm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp	(revision 7446)
+++ 	(revision )
@@ -1,208 +1,0 @@
-/*!\file:  InterpFromGridToMeshxt.cpp
- * \brief  thread core for InterpFromGridToMeshxt code
- */ 
-
-#include "./InterpFromGridToMeshx.h"
-#include "../../shared/shared.h"
-
-void* InterpFromGridToMeshxt(void* vpthread_handle){
-
-	/*gate variables :*/
-	InterpFromGridToMeshxThreadStruct* gate=NULL;
-	pthread_handle* handle=NULL;
-	int     my_thread;
-	int     num_threads;
-	
-	Vec     data_mesh=NULL;
-	double* x_mesh=NULL;
-	double* y_mesh=NULL;
-	int     x_rows,y_rows;
-	double* x=NULL;
-	double* y=NULL;
-	int     nods;
-	double  default_value;
-	double* data=NULL;
-	int     M,N;
-
-
-	/*intermediary: */
-	int     i0;
-	int     i1;
-	int     i,m,n;
-	double  x_grid;
-	double  y_grid;
-	double  area;
-	double G1,G2,G3,G4,data_value;
-	double area_1,area_2,area_3;
-	double xi,eta;
-	bool   triangle=true;
-
-
-	/*recover handle and gate: */
-	handle=(pthread_handle*)vpthread_handle;
-	gate=(InterpFromGridToMeshxThreadStruct*)handle->gate;
-	my_thread=handle->id;
-	num_threads=handle->num;
-	
-	/*recover parameters :*/
-	x_mesh=gate->x_mesh;
-	y_mesh=gate->y_mesh;
-	x_rows=gate->x_rows;
-	y_rows=gate->y_rows;
-	x=gate->x;
-	y=gate->y;
-	nods=gate->nods;
-	data_mesh=gate->data_mesh;
-	data=gate->data;
-	default_value=gate->default_value;
-	M=gate->M;
-	N=gate->N;
-
-	/*partition loop across threads: */
-	PartitionRange(&i0,&i1,nods,num_threads,my_thread);
-
-	for ( i=i0; i<i1; i++) {
-
-		x_grid=*(x_mesh+i);
-		y_grid=*(y_mesh+i);
-
-		/*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
-		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
-			
-			if(triangle){
-				/*Linear (triangle) interpolation: {{{1*/
-				/*Get area*/
-				area=(x[n+1]-x[n])*(y[m+1]-y[m]);
-
-				/*is it the upper right triangle?*/
-				/*2'     3'
-				 *+-----+
-				 *1\    |
-				 *| \   |
-				 *|  \  |
-				 *|   \ |
-				 *|    \|
-				 *2----3+1' */
-				if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
-
-					/*Then find values of data at each summit*/
-					G1=*(data+m*N+n);
-					G2=*(data+(m+1)*N+n+1);
-					G3=*(data+(m+1)*N+n);
-
-					/*Get first area coordinate*/
-					area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
-					/*Get second area coordinate*/
-					area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
-					/*Get third area coordinate = 1-area1-area2*/
-					area_3=1-area_1-area_2;
-
-					/*interpolate*/
-					data_value=area_1*G1+area_2*G2+area_3*G3;
-				}
-				else {
-
-					/*Then find values of data at each summit*/
-					G1=*(data+(m+1)*N+n+1);
-					G2=*(data+m*N+n);
-					G3=*(data+m*N+n+1);
-
-					/*Get first area coordinate*/
-					area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
-					/*Get second area coordinate*/
-					area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
-					/*Get third area coordinate = 1-area1-area2*/
-					area_3=1-area_1-area_2;
-
-					/*interpolate*/
-					data_value=area_1*G1+area_2*G2+area_3*G3;
-				}
-				/*}}}*/
-			}
-			else{
-				/*Bilinear  interpolation: (http://en.wikipedia.org/wiki/Bilinear_interpolation) {{{1*/
-
-				/*    Q12    R2        Q22
-				 * y2 x------x---------x
-				 *    |      |         |
-				 *    |      |         |
-				 *    |      +P        |
-				 *    |      |         |
-				 *    |Q11   R1        Q21
-				 * y1 x------x---------x
-				 *    x1               x2
-				 *
-				 */
-
-				/*Coordinates*/
-				double x1=x[n];
-				double x2=x[n+1];
-				double y1=y[m];
-				double y2=y[m+1];
-
-				/*Values on each summit*/
-				double Q11=*(data+m    *N+n);
-				double Q21=*(data+m    *N+n+1);
-				double Q12=*(data+(m+1)*N+n);
-				double Q22=*(data+(m+1)*N+n+1);
-
-				/*Interpolate*/
-				data_value=
-				  +Q11*(x2-x_grid)*(y2-y_grid)/((x2-x1)*(y2-y1))
-				  +Q21*(x_grid-x1)*(y2-y_grid)/((x2-x1)*(y2-y1))
-				  +Q12*(x2-x_grid)*(y_grid-y1)/((x2-x1)*(y2-y1))
-				  +Q22*(x_grid-x1)*(y_grid-y1)/((x2-x1)*(y2-y1));
-				/*}}}*/
-			}
-
-			/*Treat NANs*/
-			if (isnan(data_value)){
-				data_value=default_value;
-			}
-		}
-		else{
-			data_value=default_value;
-		}
-		VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
-	}
-
-}
-
-
-int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){
-
-	int foundx=0;
-	int foundy=0;
-	int i;
-	int m=-1;
-	int n=-1;
-
-	for (i=0;i<x_rows-1;i++){
-		if ( (*(x+i)<=xgrid) && (xgrid<*(x+i+1)) ){
-			m=i;
-			foundx= 1;
-			break;
-		}
-	}
-	if(*(x+x_rows-1)==xgrid){
-		m=x_rows-2;
-		foundx=1;
-	}
-	
-	for (i=0;i<y_rows-1;i++){
-		if ( (*(y+i)<=ygrid) && (ygrid<*(y+i+1)) ){
-			n=i;
-			foundy= 1;
-			break;
-		}
-	}
-	if(*(y+y_rows-1)==ygrid){
-		m=y_rows-2;
-		foundy=1;
-	}
-
-	/*Assign output pointers:*/
-	*pm=m;
-	*pn=n;
-	return foundx*foundy;
-}