Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 25424)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 25425)
@@ -1263,4 +1263,89 @@
 }
 /*}}}*/
+#ifdef _HAVE_AD_
+void  IoModel::FetchData(IssmPDouble** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
+
+	/*First, look if has already been loaded (might be an independent variable)*/
+	vector<IoData*>::iterator iter;
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(strcmp(iodata->name,data_name)==0){
+			_error_(data_name <<" has already been loaded as in IssmDouble and cannot be converted to IssmPDouble");
+		}
+	}
+
+	/*output: */
+	int          M,N;
+	IssmPDouble *matrix = NULL;
+	int          code   = 0;
+
+	/*recover my_rank:*/
+	int my_rank=IssmComm::GetRank();
+
+	/*Set file pointer to beginning of the data: */
+	fid=this->SetFilePointerToData(&code,NULL,data_name);
+	if(code!=5 && code!=6 && code!=7 && code!=10)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
+
+	/*Now fetch: */
+
+	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+	/*numberofelements: */
+	if(my_rank==0){
+		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+	}
+	ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	if(my_rank==0){
+		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+	}
+	ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	/*Now allocate matrix: */
+	if(M*N){
+		if(code==10){
+			/*Special case for Compressed mat*/
+			IssmPDouble offset,range;
+			if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset");
+			ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+			if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range");
+			ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+			*pmatrix=xNew<IssmPDouble>(M*N);
+
+			/*Read matrix*/
+			uint8_t* rawmatrix=xNew<uint8_t>((M-1)*N);
+			if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix ");
+			ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm());
+
+			for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast<IssmPDouble>(rawmatrix[i])/255.;
+			xDelete<uint8_t>(rawmatrix);
+
+			/*read time now*/
+			IssmPDouble* timematrix=xNew<IssmPDouble>(N);
+			if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix");
+			ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+			for(int i=0;i<N;++i) (*pmatrix)[(M-1)*N+i]=timematrix[i];
+			xDelete<IssmPDouble>(timematrix);
+
+		}
+		else{
+			/*Read matrix on node 0, then broadcast: */
+			matrix=xNew<IssmPDouble>(M*N);
+			if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+			ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+			*pmatrix=matrix;
+		}
+	}
+	else{
+		*pmatrix=NULL;
+	}
+	/*Assign output pointers: */
+	if(pM) *pM=M;
+	if(pN) *pN=N;
+}
+/*}}}*/
+#endif
 void  IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,const char* data_name){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 25424)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 25425)
@@ -136,4 +136,7 @@
 		void        FetchData(int** pmatrix,int* pM,int* pN,const char* data_name);
 		void        FetchData(IssmDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
+#if _HAVE_AD_  && !defined(_WRAPPERS_)
+		void        FetchData(IssmPDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
+#endif
 		void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
 		void        FetchData(Options *options,const char* data_name);
Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 25424)
+++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 25425)
@@ -44,5 +44,5 @@
 
 	/*Get tree type (FIXME)*/
-	IssmDouble dtree = 0.;
+	double dtree = 0.;
 	options->Get(&dtree,"treetype",1.);
 	this->treetype = reCast<int>(dtree);
@@ -552,5 +552,5 @@
 	IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
 
-	IssmDouble unbias = variogram->Covariance(0.,0.);
+	IssmPDouble unbias = variogram->Covariance(0.,0.);
 	/*First: Create semivariogram matrix for observations*/
 	for(i=0;i<n_obs;i++){
Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Observations.h	(revision 25424)
+++ /issm/trunk-jpl/src/c/classes/kriging/Observations.h	(revision 25425)
@@ -25,26 +25,26 @@
 		/*constructors, destructors*/
 		Observations();
-		Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
+		Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
 		~Observations();
 
 		/*Initialize data structures*/
-		void InitCovertree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
-		void InitQuadtree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
+		void InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
+		void InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
 
 		/*Methods*/
-		void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
-		void ClosestObservationCovertree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
-		void ClosestObservationQuadtree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
+		void ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+		void ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+		void ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
 		void Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble 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);
-		void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs);
-		void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-		void ObservationListCovertree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-		void ObservationListQuadtree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-		void QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n);
-		void Variomap(IssmDouble* gamma,IssmDouble *x,int n);
+		void InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power);
+		void InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata);
+		void InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram);
+		void InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+		void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs);
+		void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+		void ObservationListCovertree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+		void ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+		void QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n);
+		void Variomap(IssmPDouble* gamma,IssmPDouble *x,int n);
 
 };
Index: /issm/trunk-jpl/src/c/main/kriging.cpp
===================================================================
--- /issm/trunk-jpl/src/c/main/kriging.cpp	(revision 25424)
+++ /issm/trunk-jpl/src/c/main/kriging.cpp	(revision 25425)
@@ -7,5 +7,5 @@
 /*Local prototypes*/
 void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv);
-void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);
+void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid);
 
 int main(int argc,char **argv){
@@ -23,14 +23,14 @@
 	/*Input*/
 	int         ninterp,nobs;
-	IssmDouble *x        = NULL;
-	IssmDouble *y        = NULL;
-	IssmDouble *data     = NULL;
-	IssmDouble *x_interp = NULL;
-	IssmDouble *y_interp = NULL;
+	double *x        = NULL;
+	double *y        = NULL;
+	double *data     = NULL;
+	double *x_interp = NULL;
+	double *y_interp = NULL;
 	Options    *options  = NULL;
 
 	/*Output*/
-	IssmDouble *predictions = NULL;
-	IssmDouble *error       = NULL;
+	double *predictions = NULL;
+	double *error       = NULL;
 
 	/*Initialize exception trapping: */
@@ -73,11 +73,11 @@
 	xDelete<char>(outbinfilename);
 	xDelete<char>(rootpath);
-	xDelete<IssmDouble>(x);
-	xDelete<IssmDouble>(y);
-	xDelete<IssmDouble>(data);
-	xDelete<IssmDouble>(x_interp);
-	xDelete<IssmDouble>(y_interp);
-	xDelete<IssmDouble>(predictions);
-	xDelete<IssmDouble>(error);
+	xDelete<double>(x);
+	xDelete<double>(y);
+	xDelete<double>(data);
+	xDelete<double>(x_interp);
+	xDelete<double>(y_interp);
+	xDelete<double>(predictions);
+	xDelete<double>(error);
 	delete options;
 	delete results;
@@ -130,12 +130,12 @@
 }
 
-void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){
+void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){
 
-	int         ninterp    ,nobs;
-	IssmDouble *x        = NULL;
-	IssmDouble *y        = NULL;
-	IssmDouble *data     = NULL;
-	IssmDouble *x_interp = NULL;
-	IssmDouble *y_interp = NULL;
+	int     ninterp,nobs;
+	double *x        = NULL;
+	double *y        = NULL;
+	double *data     = NULL;
+	double *x_interp = NULL;
+	double *y_interp = NULL;
 	Options    *options  = NULL;
 
@@ -144,7 +144,7 @@
 	iomodel->fid=fid;
 	iomodel->CheckFile();
-	iomodel->FetchData(&x,&M,&N,"md.x");        nobs=M*N;
-	iomodel->FetchData(&y,&M,&N,"md.y");        _assert_(M*N==nobs);
-	iomodel->FetchData(&data,&M,&N,"md.data");     _assert_(M*N==nobs);
+	iomodel->FetchData(&x,&M,&N,"md.x");               nobs=M*N;
+	iomodel->FetchData(&y,&M,&N,"md.y");               _assert_(M*N==nobs);
+	iomodel->FetchData(&data,&M,&N,"md.data");         _assert_(M*N==nobs);
 	iomodel->FetchData(&x_interp,&M,&N,"md.x_interp"); ninterp=M*N;
 	iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);