Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 14047)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 14048)
@@ -424,4 +424,74 @@
 	*pprediction = obs;
 }/*}}}*/
+/*FUNCTION Observations::InterpolationV4{{{*/
+void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
+
+	/*Intermediaries*/
+	int         i,j,n_obs;
+	IssmPDouble prediction,h;
+	IssmPDouble *x       = NULL;
+	IssmPDouble *y       = NULL;
+	IssmPDouble *obs     = NULL;
+	IssmPDouble *Green   = NULL;
+	IssmPDouble *weights = NULL;
+	IssmPDouble *d       = NULL;
+
+	/*Some checks*/
+	_assert_(maxdata>0);
+	_assert_(pprediction);
+
+	/*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{
+
+		/*Allocate intermediary matrix and vectors*/
+		Green = xNew<IssmPDouble>(n_obs*n_obs);
+		d     = xNew<IssmPDouble>(n_obs);
+
+		/*First: distance vector*/
+		for(i=0;i<n_obs;i++){
+			d[i] = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+		}
+
+		/*Build Green function matrix*/
+		for(i=0;i<n_obs;i++){
+			for(j=0;j<=i;j++){
+				h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
+				if(h>0){
+					Green[j*n_obs+i] = h*h*(log(h)-1.);
+				}
+				else{
+					Green[j*n_obs+i] = 0.;
+				}
+				Green[i*n_obs+j] = Green[j*n_obs+i];
+			}
+		}
+
+		/*Compute weights*/
+#if _HAVE_GSL_
+		SolverxSeq(&weights,Green,obs,n_obs); // Green^-1 obs
+#else
+		_error_("GSL is required");
+#endif
+
+		/*ok now we can interpolate with these weights*/
+
+		_error_("Not implemented yet");
+	}
+
+	/*clean-up*/
+	*pprediction = prediction;
+	xDelete<IssmPDouble>(x);
+	xDelete<IssmPDouble>(y);
+	xDelete<IssmPDouble>(obs);
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
Index: /issm/trunk-jpl/src/c/Container/Observations.h
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.h	(revision 14047)
+++ /issm/trunk-jpl/src/c/Container/Observations.h	(revision 14048)
@@ -28,4 +28,5 @@
 		void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
 		void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
+		void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
 		void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
 		void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 14047)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 14048)
@@ -117,4 +117,23 @@
 		xDelete<double>(gate.percent);
 	}
+	else if(strcmp(output,"v4")==0){ //Inverse distance weighting
+		/*initialize thread parameters: */
+		gate.n_interp     = n_interp;
+		gate.x_interp     = x_interp;
+		gate.y_interp     = y_interp;
+		gate.radius       = radius;
+		gate.mindata      = mindata;
+		gate.maxdata      = maxdata;
+		gate.variogram    = variogram;
+		gate.observations = observations;
+		gate.predictions  = predictions;
+		gate.error        = error;
+		gate.percent      = xNewZeroInit<double>(num);
+
+		/*launch the thread manager with Krigingxt as a core: */
+		LaunchThread(v4t,(void*)&gate,num);
+		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
+		xDelete<double>(gate.percent);
+	}
 	else if(strcmp(output,"prediction")==0){
 
@@ -292,4 +311,50 @@
 	return NULL;
 }/*}}}*/
+/*FUNCTION v4t{{{*/
+void* v4t(void* vpthread_handle){
+
+	/*gate variables :*/
+	KrigingxThreadStruct *gate        = NULL;
+	pthread_handle       *handle      = NULL;
+	int my_thread;
+	int num_threads;
+	int i0,i1;
+
+	/*recover handle and gate: */
+	handle      = (pthread_handle*)vpthread_handle;
+	gate        = (KrigingxThreadStruct*)handle->gate;
+	my_thread   = handle->id;
+	num_threads = handle->num;
+
+	/*recover parameters :*/
+	int           n_interp     = gate->n_interp;
+	double       *x_interp     = gate->x_interp;
+	double       *y_interp     = gate->y_interp;
+	double        radius       = gate->radius;
+	int           mindata      = gate->mindata;
+	int           maxdata      = gate->maxdata;
+	Variogram    *variogram    = gate->variogram;
+	Observations *observations = gate->observations;
+	double       *predictions  = gate->predictions;
+	double       *error        = gate->error;
+	double       *percent      = gate->percent;
+
+	/*Intermediaries*/
+	double localpercent;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+	for(int idx=i0;idx<i1;idx++){
+
+		/*Print info*/
+		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
+		localpercent=percent[0];
+		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
+		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+
+		observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
+	}
+	return NULL;
+}/*}}}*/
 
 void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 14047)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 14048)
@@ -35,3 +35,4 @@
 void* NearestNeighbort(void*);
 void* idwt(void*);
+void* v4t(void*);
 #endif /* _KRIGINGX_H */