Index: /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp
===================================================================
--- /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp	(revision 1173)
+++ /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp	(revision 1174)
@@ -11,5 +11,5 @@
 int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
 
-int InterpFromGridx( 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) {
+int InterpFromGridx( 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) {
 
 
@@ -21,6 +21,4 @@
 	double* y=NULL;
 	int i,m,n;
-	int ind;
-	double dis;
 	double x_grid,y_grid;
 	double xi,eta;
@@ -116,27 +114,13 @@
 			}
 
-			/*Treat NANs: take the closest non nan value avomg neighbors*/
+			/*Treat NANs*/
 			if isnan(data_value){
-				ind=m*N+n;
-				dis=pow(x[n+1]-x[n],2);
-				if ( ((pow(y_grid-y[m],2)+pow(x_grid-x[n+1],2))<dis) && (!isnan(data[m*N+n+1]))){
-					ind=m*N+n+1;
-					dis=pow(y_grid-y[m],2)+pow(x_grid-x[n+1],2);
-				}
-				if ( ((pow(y_grid-y[m+1],2)+pow(x_grid-x[n],2))<dis) && (!isnan(data[(m+1)*N+n]))){
-					ind=(m+1)*N+n;
-					dis=pow(y_grid-y[m+1],2)+pow(x_grid-x[n],2);
-				}
-				if ( ((pow(y_grid-y[m+1],2)+pow(x_grid-x[n+1],2))<dis) && (!isnan(data[(m+1)*N+n+1]))){
-					ind=(m+1)*N+n+1;
-					dis=pow(y_grid-y[m+1],2)+pow(x_grid-x[n+1],2);
-				}
-				data_value=*(data+ind);
+				data_value=default_value;
 			}
 		}
 		else{
-			data_value=-9999;
+			data_value=default_value;
 		}
-		VecSetValues(data_mesh,1,&i,&data_value,INSERT_VALUES);
+		VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
 	}
 
Index: /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.h
===================================================================
--- /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.h	(revision 1173)
+++ /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.h	(revision 1174)
@@ -8,5 +8,5 @@
 #include "../toolkits/toolkits.h"
 
-int InterpFromGridx( 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);
+int InterpFromGridx( 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);
 
 #endif /* _INTERPFROMGRIDX_H */
Index: /issm/trunk/src/m/classes/public/plugvelocities.m
===================================================================
--- /issm/trunk/src/m/classes/public/plugvelocities.m	(revision 1173)
+++ /issm/trunk/src/m/classes/public/plugvelocities.m	(revision 1174)
@@ -1,3 +1,3 @@
-function md=plugvelocities(md,filename)
+function md=plugvelocities(md,filename,default_value)
 %PLUGVELOCITIES - load velocities on a model
 %
@@ -7,13 +7,13 @@
 %
 %   Usage:
-%      md=plugvelocities(md,filename)
+%      md=plugvelocities(md,filename,default_value)
 %
 %   Example:
-%      md=plugvelocities(md,'velocityfile.mat');
+%      md=plugvelocities(md,'velocityfile.mat',0);
 %
-%   See also: PLUGDATA, GRIDDATA, GRIDDATA_MESH_TO_MESH, DATAINTERP
+%   See also: INTERPFROMFILE, GRIDDATA
 
 %some checks
-if nargin~=2 | nargout~=1
+if nargin~=3 | nargout~=1
 	help plugvelocities
 	error('plugvelocities error message: bad usage');
@@ -94,12 +94,6 @@
 
 %interpolate
-if length(x)==size(vx,2)
-	md.vx_obs=griddata(x,y,vx,md.x,md.y);
-	md.vy_obs=griddata(x,y,vy,md.x,md.y);
-
-else
-	md.vx_obs=DataInterp(x,y,vx,md.x,md.y);
-	md.vy_obs=DataInterp(x,y,vy,md.x,md.y);
-end
+md.vx_obs=InterpFromGrid(x,y,vx,md.x,md.y,default_value);
+md.vy_obs=InterpFromGrid(x,y,vy,md.x,md.y,default_value);
 md.vx=md.vx_obs;
 md.vy=md.vy_obs;
Index: /issm/trunk/src/m/utils/Interp/InterpFromFile.m
===================================================================
--- /issm/trunk/src/m/utils/Interp/InterpFromFile.m	(revision 1173)
+++ /issm/trunk/src/m/utils/Interp/InterpFromFile.m	(revision 1174)
@@ -1,3 +1,3 @@
-function data_out=InterpFromFile(x,y,filename)
+function data_out=InterpFromFile(x,y,filename,default_value)
 %INTERPFROMFILE - load data and interpolate on the given nodes
 %
@@ -16,13 +16,13 @@
 %
 %   Usage:
-%      data=InterpFromFile(x,y,filename);
+%      data=InterpFromFile(x,y,filename,default_value);
 %
 %   Example:
-%      md.surface=InterpFromFile(md.x,md.y,'surfacefile.mat');
+%      md.surface=InterpFromFile(md.x,md.y,'surfacefile.mat',0);
 %
 %   See also: PLUGVELOCITIES, INTERPFROMGRID, INTERPFROMMESH2D, INTERPFROMMESH3D
 
 %some checks
-if nargin~=3 | nargout~=1
+if nargin~=4 | nargout~=1
 	help InterpFromFile
 	error('plugdata error message: bad usage');
@@ -200,5 +200,5 @@
 
 	%interpolate
-	data_out=InterpFromMesh2d(index_data,x_data,y_data,data,x,y);
+	data_out=InterpFromMesh2d(index_data,x_data,y_data,data,x,y,default_value);
 
 else
@@ -218,4 +218,4 @@
 
 	%interpolate
-	data_out=InterpFromGrid(x_data,y_data,data,x,y);
-end
+	data_out=InterpFromGrid(x_data,y_data,data,x,y,default_value);
+end
Index: /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.cpp
===================================================================
--- /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.cpp	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.cpp	(revision 1174)
@@ -27,7 +27,7 @@
 	/*input: */
 	double* x=NULL;
+	double* y=NULL;
+
 	int     x_rows;
-	
-	double* y=NULL;
 	int     y_rows;
 
@@ -40,4 +40,6 @@
 	int     x_mesh_rows;
 	int     y_mesh_rows;
+
+	double default_value;
 
 	/* output: */
@@ -56,7 +58,8 @@
 	FetchData((void**)&x_mesh,&x_mesh_rows,NULL,XMESHHANDLE,"Matrix","Mat");
 	FetchData((void**)&y_mesh,&y_mesh_rows,NULL,YMESHHANDLE,"Matrix","Mat");
+	FetchData((void**)&default_value,NULL,NULL,DEFAULTHANDLE,"Scalar",NULL);
 
 	/* Run core computations: */
-	InterpFromGridx( &data_mesh, x, x_rows,  y, y_rows, data, data_rows,data_cols, x_mesh, y_mesh, x_mesh_rows);
+	InterpFromGridx( &data_mesh, x, x_rows,  y, y_rows, data, data_rows,data_cols, x_mesh, y_mesh, x_mesh_rows,default_value);
 
 	/*Write data: */
@@ -70,9 +73,10 @@
 {
 	_printf_("   usage:\n");
-	_printf_("   data_mesh=InterpFromGrid(x,y,data,x_mesh,y_mesh);\n\n");
+	_printf_("   data_mesh=InterpFromGrid(x,y,data,x_mesh,y_mesh,defult_value);\n\n");
 	_printf_("   where:\n");
 	_printf_("      x,y: coordinates of matrix data\n");
 	_printf_("      data - matrix holding the data to be interpolated onto the mesh.\n");
 	_printf_("      x_mesh,y_mesh: coordinates of the mesh grids onto which we interpolate.\n");
+	_printf_("      default_value: default value if no interpolation is found.\n");
 	_printf_("      data_mesh:  vector of mesh interpolated data.\n");
 	_printf_("\n");
Index: /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.h
===================================================================
--- /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.h	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromGrid/InterpFromGrid.h	(revision 1174)
@@ -24,4 +24,5 @@
 #define XMESHHANDLE prhs[3]
 #define YMESHHANDLE prhs[4]
+#define DEFAULTHANDLE prhs[5]
 
 /* serial output macros: */
@@ -32,5 +33,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  5
+#define NRHS  6
 
 #endif  /* _INTERPFROMGRId_H */
Index: /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp
===================================================================
--- /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp	(revision 1174)
@@ -43,4 +43,6 @@
 	int     y_prime_rows;
 
+	double default_value;
+
 	/*Intermediary*/
 	int i,j;
@@ -71,4 +73,5 @@
 	FetchData((void**)&x_prime,&x_prime_rows,NULL,XPRIMEHANDLE,"Matrix","Mat");
 	FetchData((void**)&y_prime,&y_prime_rows,NULL,YPRIMEHANDLE,"Matrix","Mat");
+	FetchData((void**)&default_value,NULL,NULL,DEFAULTHANDLE,"Scalar",NULL);
 
 	/* Run core computations: */
@@ -76,4 +79,7 @@
 
 	/*some checks*/
+	if (index_data_rows<1 || x_data_rows<3 || y_data_rows<3){
+		throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
+	}
 	if (x_data_rows!=y_data_rows){
 		throw ErrorException(__FUNCT__,"vectors x and y should have the same length!");
@@ -151,4 +157,5 @@
 					data_value=data[i];
 				}
+				if isnan(data_value) data_value=default_value;
 
 				/*insert value and go to the next point*/
@@ -169,5 +176,5 @@
 {
 	_printf_("   usage:\n");
-	_printf_("      data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);\n\n");
+	_printf_("      data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);\n\n");
 	_printf_("   where:\n");
 	_printf_("      x,y: coordinates of the nodes where data is defined\n");
@@ -175,4 +182,5 @@
 	_printf_("      data - vector holding the data to be interpolated onto the points.\n");
 	_printf_("      x_prime,y_prime: coordinates of the mesh grids onto which we interpolate.\n");
+	_printf_("      default_value - default value if no interpolation is found.\n");
 	_printf_("      data_prime:  vector of prime interpolated data.\n");
 	_printf_("\n");
Index: /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h
===================================================================
--- /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h	(revision 1174)
@@ -25,4 +25,5 @@
 #define XPRIMEHANDLE prhs[4]
 #define YPRIMEHANDLE prhs[5]
+#define DEFAULTHANDLE prhs[6]
 
 /* serial output macros: */
@@ -33,5 +34,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  6
+#define NRHS  7
 
-#endif  /* _DATAINTERP2_H */
+#endif  /* _INTERPFROMMESH2D_H */
Index: /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp
===================================================================
--- /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp	(revision 1174)
@@ -46,4 +46,6 @@
 	int     y_prime_rows;
 	int     z_prime_rows;
+
+	double  default_value;
 
 	/*Intermediary*/
@@ -78,4 +80,5 @@
 	FetchData((void**)&y_prime,&y_prime_rows,NULL,YPRIMEHANDLE,"Matrix","Mat");
 	FetchData((void**)&z_prime,&z_prime_rows,NULL,ZPRIMEHANDLE,"Matrix","Mat");
+	FetchData((void**)&default_value,NULL,NULL,DEFAULTHANDLE,"Scalar",NULL);
 
 	/*some checks*/
@@ -166,4 +169,5 @@
 						data_value=data[i];
 					}
+					if isnan(data_value) data_value=default_value;
 
 					/*insert value and go to the next point*/
@@ -185,5 +189,5 @@
 {
 	_printf_("   usage:\n");
-	_printf_("      data_prime=InterpFromMesh3d(index,x,y,z,data,x_prime,y_prime,z_prime);\n\n");
+	_printf_("      data_prime=InterpFromMesh3d(index,x,y,z,data,x_prime,y_prime,z_prime,default_value);\n\n");
 	_printf_("   where:\n");
 	_printf_("      x,y,z: coordinates of the nodes where data is defined\n");
@@ -191,4 +195,5 @@
 	_printf_("      data - vector holding the data to be interpolated onto the points.\n");
 	_printf_("      x_prime,y_prime,z_prime: coordinates of the mesh grids onto which we interpolate.\n");
+	_printf_("      default_value - default value if no interpolation is found.\n");
 	_printf_("      data_prime:  vector of prime interpolated data.\n");
 	_printf_("\n");
Index: /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h
===================================================================
--- /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h	(revision 1173)
+++ /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h	(revision 1174)
@@ -27,4 +27,5 @@
 #define YPRIMEHANDLE prhs[6]
 #define ZPRIMEHANDLE prhs[7]
+#define DEFAULTHANDLE prhs[8]
 
 /* serial output macros: */
@@ -35,5 +36,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  8
+#define NRHS  9
 
 #endif  /* _INTERPFROMMESH3D_H */