Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12421)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12422)
@@ -276,4 +276,53 @@
 	*pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::InterpolationIDW{{{*/
+void Observations::InterpolationIDW(double *pprediction,double x_interp,double y_interp,double radius,int mindata,int maxdata,double power){
+
+	/*Intermediaries*/
+	int    i,n_obs;
+	double prediction;
+	double numerator,denominator,h,weight;
+	double *x   = NULL;
+	double *y   = NULL;
+	double *obs = NULL;
+
+	/*Some checks*/
+	_assert_(maxdata>0);
+	_assert_(pprediction);
+	_assert_(power>0);
+
+	/*If radius is not provided or is 0, return all observations*/
+	if(radius==0) radius=this->quadtree->root->length;
+
+	/*Get list of observations for current point*/
+	this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+	/*If we have less observations than mindata, return UNDEF*/
+	if(n_obs<mindata){
+		prediction = UNDEF; 
+	}
+	else{
+		numerator   = 0.;
+		denominator = 0.;
+		for(i=0;i<n_obs;i++){
+			h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
+			if (h<0.0000001){
+				numerator   = obs[i];
+				denominator = 1.;
+				break;
+			}
+			weight = 1./pow(h,power);
+			numerator   += weight*obs[i];
+			denominator += weight;
+		}
+		prediction = numerator/denominator; 
+	}
+
+	/*clean-up*/
+	*pprediction = prediction;
+	xfree((void**)&x);
+	xfree((void**)&y);
+	xfree((void**)&obs);
+}/*}}}*/
 /*FUNCTION Observations::InterpolationKriging{{{*/
 void Observations::InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram){
@@ -363,4 +412,16 @@
 
 }/*}}}*/
+/*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
+void Observations::InterpolationNearestNeighbor(double *pprediction,double x_interp,double y_interp,double radius){
+
+	/*Intermediaries*/
+	double        x,y,obs;
+
+	/*Get clostest observation*/
+	this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
+
+	/*Assign output pointer*/
+	*pprediction = obs;
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){
Index: /issm/trunk-jpl/src/c/Container/Observations.h
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.h	(revision 12421)
+++ /issm/trunk-jpl/src/c/Container/Observations.h	(revision 12422)
@@ -25,5 +25,7 @@
 		/*Methods*/
 		void ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius);
+		void InterpolationIDW(double *pprediction,double x_interp,double y_interp,double radius,int mindata,int maxdata,double power);
 		void InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram);
+		void InterpolationNearestNeighbor(double *pprediction,double x_interp,double y_interp,double radius);
 		void ObservationList(double **px,double **py,double **pobs,int* pnobs);
 		void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata);
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12421)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12422)
@@ -228,6 +228,5 @@
 
 	/*Intermediaries*/
-	int           i,j,n_obs;
-	double        x,y,obs;
+	int           i;
 	double        localpercent;
 
@@ -242,7 +241,5 @@
 		printf("\r      interpolation progress: %5.2lf%%",localpercent);
 
-		/*Get closest observation to current point*/
-		observations->ClosestObservation(&x,&y,&obs,x_interp[idx],y_interp[idx],radius);
-		predictions[idx] = obs; 
+		observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
 	}
 
@@ -279,10 +276,5 @@
 
 	/*Intermediaries*/
-	int           i,j,n_obs;
-	double        localpercent;
-	double        numerator,denominator,h,weight;
-	double *x     = NULL;
-	double *y     = NULL;
-	double *obs   = NULL;
+	double localpercent;
 	double  power = 2.;
 
@@ -294,25 +286,9 @@
 		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
 		localpercent=percent[0];
-		for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
+		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
 		printf("\r      interpolation progress: %5.2lf%%",localpercent);
 
-		/*Get closest observation to current point*/
-		observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
-		if(n_obs){
-			numerator   = 0;
-			denominator = 0;
-			for(j=0;j<n_obs;j++){
-				h      = sqrt( (x[j]-x_interp[idx])*(x[j]-x_interp[idx]) + (y[j]-y_interp[idx])*(y[j]-y_interp[idx]));
-				weight = 1./pow(h,power);
-				numerator   += weight*obs[j];
-				denominator += weight;
-			}
-			predictions[idx] = numerator/denominator; 
-		}
-		else{
-			predictions[idx] = UNDEF; 
-		}
-	}
-
+		observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
+	}
 	return NULL;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp	(revision 12421)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp	(revision 12422)
@@ -9,9 +9,5 @@
 #include "../../objects/objects.h"
 #include "../../Container/Container.h"
-#include "../modules.h"
-
-#ifdef _HAVE_GSL_
-#include <gsl/gsl_linalg.h>
-#endif
+#include "../../io/io.h"
 
 /*FUNCTION pKrigingx{{{*/
@@ -32,6 +28,5 @@
 	Observations *observations = NULL;
 
-	/*Get Variogram from Options*/
-	ProcessVariogram(&variogram,options);
+	/*Get some Options*/
 	options->Get(&radius,"searchradius",0.);
 	options->Get(&mindata,"mindata",1);
@@ -56,4 +51,7 @@
 	else if(strcmp(output,"prediction")==0){
 
+		/*Process Variogram*/
+		ProcessVariogram(&variogram,options);
+
 		/*partition loop across threads: */
 		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
@@ -72,8 +70,39 @@
 #endif
 	}
+	else if(strcmp(output,"nearestneighbor")==0){
+
+		/*partition loop across threads: */
+		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
+			_printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
+			observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
+		}
+		_printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
+
+#ifdef _HAVE_MPI_
+		double *sumpredictions =(double*)xmalloc(n_interp*sizeof(double));
+		MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+		xfree((void**)&predictions); predictions=sumpredictions;
+#endif
+	}
+	else if(strcmp(output,"idw")==0){
+		double power;
+		options->Get(&power,"power",2.);
+
+		/*partition loop across threads: */
+		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
+			_printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
+			observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
+		}
+		_printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
+
+#ifdef _HAVE_MPI_
+		double *sumpredictions =(double*)xmalloc(n_interp*sizeof(double));
+		MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+		xfree((void**)&predictions); predictions=sumpredictions;
+#endif
+	}
 	else{
 		_error_("output '%s' not supported yet",output);
 	}
-	_printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
 
 	/*clean-up and Assign output pointer*/
@@ -105,42 +134,2 @@
 	*pvariogram = variogram;
 }/*}}}*/
-void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/
-#ifdef _HAVE_GSL_
-
-		/*GSL Matrices and vectors: */
-		int              s;
-		gsl_matrix_view  a;
-		gsl_vector_view  b;
-		gsl_vector      *x = NULL;
-		gsl_permutation *p = NULL;
-
-		/*A will be modified by LU decomposition. Use copy*/
-		double* Acopy = (double*)xmalloc(n*n*sizeof(double));
-		memcpy(Acopy,A,n*n*sizeof(double));
-
-		/*Initialize gsl matrices and vectors: */
-		a = gsl_matrix_view_array (Acopy,n,n);
-		b = gsl_vector_view_array (B,n);
-		x = gsl_vector_alloc (n);
-
-		/*Run LU and solve: */
-		p = gsl_permutation_alloc (n);
-		gsl_linalg_LU_decomp (&a.matrix, p, &s);
-		gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
-
-		//printf ("x = \n");
-		//gsl_vector_fprintf (stdout, x, "%g");
-
-		/*Copy result*/
-		double* X = (double*)xmalloc(n*sizeof(double));
-		memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
-
-		/*Clean up and assign output pointer*/
-		xfree((void**)&Acopy);
-		gsl_permutation_free(p);
-		gsl_vector_free(x);
-		*pX=X;
-#else
-		_error_("GSL support required");
-#endif
-	}/*}}}*/