Index: /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 7438)
+++ /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 7439)
@@ -34,5 +34,5 @@
 
 	/*Some checks on arguments: */
-	if ((M<=2) || (N<=2) || (nods<=0)){
+	if ((M<2) || (N<2) || (nods<=0)){
 		_error_("nothing to be done according to the dimensions of input matrices and vectors.");
 	}
Index: /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp	(revision 7438)
+++ /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp	(revision 7439)
@@ -36,4 +36,5 @@
 	double area_1,area_2,area_3;
 	double xi,eta;
+	bool   triangle=true;
 
 
@@ -61,5 +62,4 @@
 	PartitionRange(&i0,&i1,nods,num_threads,my_thread);
 
-	/*Linear (triangle) interpolation: */
 	for ( i=i0; i<i1; i++) {
 
@@ -70,49 +70,88 @@
 		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
 			
-			/*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;
+			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 {
-
-				/*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));
+				/*}}}*/
 			}