Index: /issm/trunk-jpl/src/c/solutions/control_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/control_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/control_core.cpp	(revision 12477)
@@ -20,5 +20,5 @@
 	int     num_controls,num_responses;
 	int     nsteps;
-	double  tol_cm;
+	IssmDouble  tol_cm;
 	bool    cm_gradient;
 	int     dim;
@@ -28,11 +28,11 @@
 
 	int*    control_type = NULL;
-	double* responses=NULL;
+	IssmDouble* responses=NULL;
 	int*    step_responses=NULL;
-	double* maxiter=NULL;
-	double* cm_jump=NULL;
+	IssmDouble* maxiter=NULL;
+	IssmDouble* cm_jump=NULL;
 		
 	/*intermediary: */
-	double  search_scalar=1;
+	IssmDouble  search_scalar=1;
 	OptArgs optargs;
 	OptPars optpars;
@@ -43,5 +43,5 @@
 
 	/*output: */
-	double* J=NULL;
+	IssmDouble* J=NULL;
 
 	/*Recover parameters used throughout the solution*/
@@ -70,5 +70,5 @@
 
 	/*Initialize responses: */
-	J=xNew<double>(nsteps);
+	J=xNew<IssmDouble>(nsteps);
 	step_responses=xNew<int>(num_responses);
 		
@@ -122,7 +122,7 @@
 	xDelete<int>(control_type);
 	xDelete<int>(step_responses);
-	xDelete<double>(maxiter);
-	xDelete<double>(responses);
-	xDelete<double>(cm_jump);
-	xDelete<double>(J);
+	xDelete<IssmDouble>(maxiter);
+	xDelete<IssmDouble>(responses);
+	xDelete<IssmDouble>(cm_jump);
+	xDelete<IssmDouble>(J);
 }
Index: /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp	(revision 12477)
@@ -17,5 +17,5 @@
 #include "./solutions.h"
 
-bool controlconvergence(double J, double tol_cm){
+bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
 
 	int i;
Index: /issm/trunk-jpl/src/c/solutions/controlrestart.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controlrestart.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/controlrestart.cpp	(revision 12477)
@@ -7,5 +7,5 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-void controlrestart(FemModel* femmodel,double* J){
+void controlrestart(FemModel* femmodel,IssmDouble* J){
 
 	int      num_controls;
Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 12477)
@@ -17,5 +17,5 @@
 
 /*Local prototype*/
-int FormFunctionGradient(TaoSolver,Vec,double*,Vec,void*);
+int FormFunctionGradient(TaoSolver,Vec,IssmDouble*,Vec,void*);
 int IssmMonitor(TaoSolver,void*);
 typedef struct {
@@ -31,5 +31,5 @@
 	AppCtx     user;
 	TaoSolver  tao;
-	double    *dummy          = NULL;
+	IssmDouble    *dummy          = NULL;
 	int       *control_list   = NULL;
 	Vector    *X              = NULL;
@@ -50,5 +50,5 @@
 	femmodel->parameters->FindParam(&dummy,NULL,NULL,InversionMaxiterPerStepEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
-	maxiter=nsteps*(int)dummy[0]; xDelete<double>(dummy);
+	maxiter=nsteps*(int)dummy[0]; xDelete<IssmDouble>(dummy);
 
 	/*Initialize TAO*/
@@ -100,5 +100,5 @@
 	TaoFinalize();
 }
-int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, double *fcn,Vec G,void *userCtx){
+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
 
 	/*Retreive arguments*/
@@ -107,5 +107,5 @@
 	FemModel *femmodel       = user->femmodel;
 	int      *cost_functions = NULL;
-	double   *cost_functionsd= NULL;
+	IssmDouble   *cost_functionsd= NULL;
 	Vector   *gradient       = NULL;
 	Vector   *X              = NULL;
@@ -147,5 +147,5 @@
 	/*Clean-up and return*/
 	xDelete<int>(cost_functions);
-	xDelete<double>(cost_functionsd);
+	xDelete<IssmDouble>(cost_functionsd);
 	return 0;
 }
@@ -153,5 +153,5 @@
 
 	int       i,its,num_responses;
-	double    f,gnorm,cnorm,xdiff;
+	IssmDouble    f,gnorm,cnorm,xdiff;
 	AppCtx   *user      = (AppCtx *)userCtx;
 	FemModel *femmodel  = user->femmodel;
Index: /issm/trunk-jpl/src/c/solutions/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/convergence.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/convergence.cpp	(revision 12477)
@@ -19,15 +19,15 @@
 	Vector* KUoldF=NULL;
 	Vector* duf=NULL;
-	double ndu,nduinf,nu;
-	double nKUF;
-	double nKUoldF;
-	double nF;
-	double solver_residue,res;
+	IssmDouble ndu,nduinf,nu;
+	IssmDouble nKUF;
+	IssmDouble nKUoldF;
+	IssmDouble nF;
+	IssmDouble solver_residue,res;
 
 	/*convergence options*/
-	double eps_res;
-	double eps_rel;
-	double eps_abs;
-	double yts;
+	IssmDouble eps_res;
+	IssmDouble eps_rel;
+	IssmDouble eps_abs;
+	IssmDouble yts;
 
 	/*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from 
Index: /issm/trunk-jpl/src/c/solutions/gradient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/gradient_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/gradient_core.cpp	(revision 12477)
@@ -16,6 +16,6 @@
 
 	/*Intermediaries*/
-	double  norm_inf;
-	double *norm_list    = NULL;
+	IssmDouble  norm_inf;
+	IssmDouble *norm_list    = NULL;
 	Vector*     new_gradient = NULL;
 	Vector*     gradient     = NULL;
@@ -48,4 +48,4 @@
 
 	/*Clean up and return*/
-	xDelete<double>(norm_list);
+	xDelete<IssmDouble>(norm_list);
 }
Index: /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp	(revision 12477)
@@ -17,8 +17,8 @@
 
 	/*intermediary*/
-	double time;
+	IssmDouble time;
 	int    nsteps;
-	double starttime,final_time;
-	double dt;
+	IssmDouble starttime,final_time;
+	IssmDouble dt;
 	bool   save_results;
 	int    output_frequency;
Index: /issm/trunk-jpl/src/c/solutions/hydrology_core_step.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/hydrology_core_step.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/hydrology_core_step.cpp	(revision 12477)
@@ -13,5 +13,5 @@
 #include "../solvers/solvers.h"
 
-void hydrology_core_step(FemModel* femmodel,int step, double time){
+void hydrology_core_step(FemModel* femmodel,int step, IssmDouble time){
 	
 	bool modify_loads=true;
Index: /issm/trunk-jpl/src/c/solutions/issm.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/issm.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/issm.cpp	(revision 12477)
@@ -30,7 +30,7 @@
 
 	/*time*/
-	double   start, finish;
-	double   start_core, finish_core;
-	double   start_init, finish_init;
+	IssmPDouble   start, finish;
+	IssmPDouble   start_core, finish_core;
+	IssmPDouble   start_init, finish_init;
 	int      ierr;
 
@@ -50,5 +50,5 @@
 	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
 	#else
-	start=(double)clock();
+	start=(IssmPDouble)clock();
 	#endif
 
@@ -73,5 +73,5 @@
 	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
 	#else
-	start_init=(double)clock();
+	start_init=(IssmPDouble)clock();
 	#endif
 	femmodel=new FemModel(binfilename,outbinfilename,solution_type,analyses,numanalyses);
@@ -97,5 +97,5 @@
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
 	#else
-	finish_init=(double)clock();
+	finish_init=(IssmPDouble)clock();
 	#endif
 
@@ -104,5 +104,5 @@
 	MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
 	#else
-	start_core=(double)clock();
+	start_core=(IssmPDouble)clock();
 	#endif
 	
@@ -130,5 +130,5 @@
 	MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	#else
-	finish_core=(double)clock();
+	finish_core=(IssmPDouble)clock();
 	#endif
 	
@@ -158,5 +158,5 @@
 	_printf_(true,"\n   %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
 	#else
-	finish=(double)clock();
+	finish=(IssmPDouble)clock();
 	_printf_(true,"\n   %-34s %f seconds  \n","FemModel initialization elapsed time:",(finish_init-start_init)/CLOCKS_PER_SEC);
 	_printf_(true,"   %-34s %f seconds  \n","Core solution elapsed time:",(finish_core-start_core)/CLOCKS_PER_SEC);
Index: /issm/trunk-jpl/src/c/solutions/kriging.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/kriging.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/kriging.cpp	(revision 12477)
@@ -8,5 +8,5 @@
 /*Local prototypes*/
 void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv);
-void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid);
+void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);
 
 int main(int argc,char **argv){
@@ -24,14 +24,14 @@
 	/*Input*/
 	int      ninterp,nobs;
-	double  *x        = NULL;
-	double  *y        = NULL;
-	double  *data     = NULL;
-	double  *x_interp = NULL;
-	double  *y_interp = NULL;
+	IssmDouble  *x        = NULL;
+	IssmDouble  *y        = NULL;
+	IssmDouble  *data     = NULL;
+	IssmDouble  *x_interp = NULL;
+	IssmDouble  *y_interp = NULL;
 	Options *options  = NULL;
 
 	/*Output*/
-	double *predictions = NULL;
-	double *error       = NULL;
+	IssmDouble *predictions = NULL;
+	IssmDouble *error       = NULL;
 
 	ISSMBOOT();
@@ -89,11 +89,11 @@
 	xDelete<char>(binfilename);
 	xDelete<char>(outbinfilename);
-	xDelete<double>(x);
-	xDelete<double>(y);
-	xDelete<double>(data);
-	xDelete<double>(x_interp);
-	xDelete<double>(y_interp);
-	xDelete<double>(predictions);
-	xDelete<double>(error);
+	xDelete<IssmDouble>(x);
+	xDelete<IssmDouble>(y);
+	xDelete<IssmDouble>(data);
+	xDelete<IssmDouble>(x_interp);
+	xDelete<IssmDouble>(y_interp);
+	xDelete<IssmDouble>(predictions);
+	xDelete<IssmDouble>(error);
 	delete options;
 	delete results;
@@ -134,12 +134,12 @@
 }
 
-void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){
+void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){
 
 	int      ninterp,nobs,numoptions;
-	double  *x        = NULL;
-	double  *y        = NULL;
-	double  *data     = NULL;
-	double  *x_interp = NULL;
-	double  *y_interp = NULL;
+	IssmDouble  *x        = NULL;
+	IssmDouble  *y        = NULL;
+	IssmDouble  *data     = NULL;
+	IssmDouble  *x_interp = NULL;
+	IssmDouble  *y_interp = NULL;
 	Options *options  = NULL;
 	Option  *option   = NULL;
Index: /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/objectivefunction.cpp	(revision 12477)
@@ -20,10 +20,10 @@
 /*}}}*/
 
-double objectivefunction(double search_scalar,OptArgs* optargs){
+IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs){
 
 	int i;  
 	
 	/*output: */
-	double J;
+	IssmDouble J;
 	
 	/*parameters: */
Index: /issm/trunk-jpl/src/c/solutions/solutions.h
===================================================================
--- /issm/trunk-jpl/src/c/solutions/solutions.h	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/solutions.h	(revision 12477)
@@ -18,5 +18,5 @@
 void diagnostic_core(FemModel* femmodel);
 void hydrology_core(FemModel* femmodel);
-void hydrology_core_step(FemModel* femmodel,int step, double time);
+void hydrology_core_step(FemModel* femmodel,int step, IssmDouble time);
 void thermal_core(FemModel* femmodel);
 void enthalpy_core(FemModel* femmodel);
@@ -30,18 +30,18 @@
 void steadystate_core(FemModel* femmodel);
 void transient_core(FemModel* femmodel);
-double objectivefunction(double search_scalar,OptArgs* optargs);
+IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs);
 
 //convergence:
 void convergence(bool* pconverged, Matrix* K_ff,Vector* p_f,Vector* u_f,Vector* u_f_old,Parameters* parameters);
-bool controlconvergence(double J,double tol_cm);
+bool controlconvergence(IssmDouble J,IssmDouble tol_cm);
 bool steadystateconvergence(FemModel* femmodel);
 
 //optimization
-int GradJSearch(double* search_vector,FemModel* femmodel,int step);
+int GradJSearch(IssmDouble* search_vector,FemModel* femmodel,int step);
 
 //diverse
 void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ppetscname,char** plockname,int argc,char **argv);
 void WriteLockFile(char* filename);
-void controlrestart(FemModel* femmodel,double* J);
+void controlrestart(FemModel* femmodel,IssmDouble* J);
 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
 
Index: /issm/trunk-jpl/src/c/solutions/steadystateconvergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/steadystateconvergence.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/steadystateconvergence.cpp	(revision 12477)
@@ -28,5 +28,5 @@
 	int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum};
 	int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum 
-	double convergencecriterionvalue[1]; //value of criterion to be respected
+	IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected
 
 	/*retrieve parameters: */
Index: /issm/trunk-jpl/src/c/solutions/thermal_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/thermal_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/thermal_core.cpp	(revision 12477)
@@ -16,5 +16,5 @@
 
 	/*intermediary*/
-	double melting_offset;
+	IssmDouble melting_offset;
 	bool   save_results;
 	bool   dakota_analysis  = false;
Index: /issm/trunk-jpl/src/c/solutions/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/transient_core.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solutions/transient_core.cpp	(revision 12477)
@@ -23,5 +23,5 @@
 
 	/*parameters: */
-	double starttime,finaltime,dt,yts;
+	IssmDouble starttime,finaltime,dt,yts;
 	bool   isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy;
 	bool   save_results,dakota_analysis;
@@ -34,5 +34,5 @@
 	/*intermediary: */
 	int    step;
-	double time;
+	IssmDouble time;
 
 	//first recover parameters common to all solutions
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 12477)
@@ -17,5 +17,5 @@
 	int    num_unstable_constraints;
 	int    count;
-	double kmax;
+	IssmDouble kmax;
 	Matrix* Kff = NULL;
 	Matrix* Kfs    = NULL;
@@ -70,5 +70,5 @@
 			bool max_iteration_state=false;
 			int tempStep=1;
-			double tempTime=1.0;
+			IssmDouble tempTime=1.0;
 			femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
 			break;
@@ -78,5 +78,5 @@
 			bool max_iteration_state=true;
 			int tempStep=1;
-			double tempTime=1.0;
+			IssmDouble tempTime=1.0;
 			femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
 			break;
Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 12477)
@@ -88,5 +88,5 @@
 			bool max_iteration_state=false;
 			int tempStep=1;
-			double tempTime=1.0;
+			IssmDouble tempTime=1.0;
 			femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
 			break;
@@ -99,5 +99,5 @@
 			bool max_iteration_state=true;
 			int tempStep=1;
-			double tempTime=1.0;
+			IssmDouble tempTime=1.0;
 			femmodel->results->AddObject(new BoolExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
 			break;
Index: /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 12477)
@@ -16,5 +16,5 @@
 	Vector* tf_old=NULL; 
 	Vector* ys=NULL; 
-	double melting_offset;
+	IssmDouble melting_offset;
 
 	/*intermediary: */
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp	(revision 12477)
@@ -33,18 +33,18 @@
 	this->N=pN;
 	this->matrix=NULL;
-	if(M*N) this->matrix=xNewZeroInit<double>(pM*pN);
-}
-/*}}}*/
-/*FUNCTION SeqMat::SeqMat(int M,int N, double sparsity){{{*/
-SeqMat::SeqMat(int pM,int pN, double sparsity){
-
-	this->M=pM;
-	this->N=pN;
-	this->matrix=NULL;
-	if(M*N) this->matrix=xNewZeroInit<double>(pM*pN);
-}
-/*}}}*/
-/*FUNCTION SeqMat(double* serial_mat,int M,int N,double sparsity){{{*/
-SeqMat::SeqMat(double* serial_mat,int pM,int pN,double sparsity){
+	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
+}
+/*}}}*/
+/*FUNCTION SeqMat::SeqMat(int M,int N, IssmDouble sparsity){{{*/
+SeqMat::SeqMat(int pM,int pN, IssmDouble sparsity){
+
+	this->M=pM;
+	this->N=pN;
+	this->matrix=NULL;
+	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
+}
+/*}}}*/
+/*FUNCTION SeqMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity){{{*/
+SeqMat::SeqMat(IssmDouble* serial_mat,int pM,int pN,IssmDouble sparsity){
 
 	int i,j;
@@ -54,6 +54,6 @@
 	this->matrix=NULL;
 	if(M*N){
-		this->matrix=xNewZeroInit<double>(pM*pN);
-		memcpy(this->matrix,serial_mat,pM*pN*sizeof(double));
+		this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
+		xMemCpy<IssmDouble>(this->matrix,serial_mat,pM*pN);
 	}
 
@@ -66,5 +66,5 @@
 	this->N=pN;
 	this->matrix=NULL;
-	if(M*N) this->matrix=xNewZeroInit<double>(pM*pN);
+	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
 }
 /*}}}*/
@@ -72,5 +72,5 @@
 SeqMat::~SeqMat(){
 
-	xDelete<double>(this->matrix);
+	xDelete<IssmDouble>(this->matrix);
 	M=0;
 	N=0;
@@ -100,8 +100,8 @@
 /*}}}*/
 /*FUNCTION SeqMat::Norm{{{*/
-double SeqMat::Norm(NormMode mode){
-
-	double norm;
-	double absolute;
+IssmDouble SeqMat::Norm(NormMode mode){
+
+	IssmDouble norm;
+	IssmDouble absolute;
 	int i,j;
 
@@ -145,5 +145,5 @@
 	int i,j;
 	int XM,AXM;
-	double dummy;
+	IssmDouble dummy;
 
 	X->GetSize(&XM);
@@ -166,5 +166,5 @@
 SeqMat* SeqMat::Duplicate(void){
 
-	double dummy=0;
+	IssmDouble dummy=0;
 
 	return new SeqMat(this->matrix,this->M,this->N,dummy);
@@ -173,11 +173,11 @@
 /*}}}*/
 /*FUNCTION SeqMat::ToSerial{{{*/
-double* SeqMat::ToSerial(void){
-
-	double* buffer=NULL;
+IssmDouble* SeqMat::ToSerial(void){
+
+	IssmDouble* buffer=NULL;
 
 	if(this->M*this->N){
-		buffer=xNew<double>(this->M*this->N);
-		memcpy(buffer,this->matrix,this->M*this->N*sizeof(double));
+		buffer=xNew<IssmDouble>(this->M*this->N);
+		xMemCpy<IssmDouble>(buffer,this->matrix,this->M*this->N);
 	}
 	return buffer;
@@ -186,5 +186,5 @@
 /*}}}*/
 /*FUNCTION SeqMat::SetValues{{{*/
-void SeqMat::SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode){
+void SeqMat::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
 	
 	int i,j;
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h	(revision 12476)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h	(revision 12477)
@@ -1,4 +1,4 @@
 /*!\file:  SeqMat.h
- * \brief wrapper to SeqMat objects, which are just wrappers to a simple double* buffer.
+ * \brief wrapper to SeqMat objects, which are just wrappers to a simple IssmDouble* buffer.
  */ 
 
@@ -24,11 +24,11 @@
 	
 		int M,N; 
-		double* matrix; 
+		IssmDouble* matrix; 
 
 		/*SeqMat constructors, destructors {{{*/
 		SeqMat();
 		SeqMat(int M,int N);
-		SeqMat(int M,int N,double sparsity);
-		SeqMat(double* serial_mat,int M,int N,double sparsity);
+		SeqMat(int M,int N,IssmDouble sparsity);
+		SeqMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity);
 		SeqMat(int M,int N,int connectivity,int numberofdofspernode);
 		~SeqMat();
@@ -37,11 +37,11 @@
 		void Echo(void);
 		void Assemble(void);
-		double Norm(NormMode norm_type);
+		IssmDouble Norm(NormMode norm_type);
 		void GetSize(int* pM,int* pN);
 		void GetLocalSize(int* pM,int* pN);
 		void MatMult(SeqVec* X,SeqVec* AX);
 		SeqMat* Duplicate(void);
-		double* ToSerial(void);
-		void SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode);
+		IssmDouble* ToSerial(void);
+		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode);
 		void Convert(MatrixType type);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp	(revision 12476)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp	(revision 12477)
@@ -31,9 +31,9 @@
 	this->M=pM;
 	this->vector=NULL;
-	if(this->M) this->vector=xNewZeroInit<double>(pM);
-}
-/*}}}*/
-/*FUNCTION SeqVec::SeqVec(double* serial_vec,int M){{{*/
-SeqVec::SeqVec(double* buffer,int pM){
+	if(this->M) this->vector=xNewZeroInit<IssmDouble>(pM);
+}
+/*}}}*/
+/*FUNCTION SeqVec::SeqVec(IssmDouble* serial_vec,int M){{{*/
+SeqVec::SeqVec(IssmDouble* buffer,int pM){
 
 	int i,j;
@@ -42,6 +42,6 @@
 	this->vector=NULL;
 	if(this->M){
-		this->vector=xNewZeroInit<double>(pM);
-		memcpy(this->vector,buffer,pM*sizeof(double));
+		this->vector=xNewZeroInit<IssmDouble>(pM);
+		xMemCpy<IssmDouble>(this->vector,buffer,pM);
 	}
 }
@@ -49,5 +49,5 @@
 		/*FUNCTION SeqVec::~SeqVec(){{{*/
 SeqVec::~SeqVec(){
-	xDelete<double>(this->vector);
+	xDelete<IssmDouble>(this->vector);
 	M=0;
 }
@@ -74,5 +74,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::SetValues{{{*/
-void SeqVec::SetValues(int ssize, int* list, double* values, InsMode mode){
+void SeqVec::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){
 	
 	int i;
@@ -92,5 +92,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::SetValue{{{*/
-void SeqVec::SetValue(int dof, double value, InsMode mode){
+void SeqVec::SetValue(int dof, IssmDouble value, InsMode mode){
 
 	switch(mode){
@@ -108,5 +108,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::GetValue{{{*/
-void SeqVec::GetValue(double* pvalue,int dof){
+void SeqVec::GetValue(IssmDouble* pvalue,int dof){
 
 	*pvalue=this->vector[dof];
@@ -137,5 +137,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::Set{{{*/
-void SeqVec::Set(double value){
+void SeqVec::Set(IssmDouble value){
 
 	int i;
@@ -145,5 +145,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::AXPY{{{*/
-void SeqVec::AXPY(SeqVec* X, double a){
+void SeqVec::AXPY(SeqVec* X, IssmDouble a){
 
 	int i;
@@ -155,5 +155,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::AYPX{{{*/
-void SeqVec::AYPX(SeqVec* X, double a){
+void SeqVec::AYPX(SeqVec* X, IssmDouble a){
 	
 	int i;
@@ -165,11 +165,11 @@
 /*}}}*/
 /*FUNCTION SeqVec::ToMPISerial{{{*/
-double* SeqVec::ToMPISerial(void){
-
-	double* buffer=NULL;
+IssmDouble* SeqVec::ToMPISerial(void){
+
+	IssmDouble* buffer=NULL;
 
 	if(this->M){
-		buffer=xNew<double>(this->M);
-		memcpy(buffer,this->vector,this->M*sizeof(double));
+		buffer=xNew<IssmDouble>(this->M);
+		xMemCpy<IssmDouble>(buffer,this->vector,this->M);
 	}
 	return buffer;
@@ -188,7 +188,7 @@
 /*}}}*/
 /*FUNCTION SeqVec::Norm{{{*/
-double SeqVec::Norm(NormMode mode){
-
-	double norm;
+IssmDouble SeqVec::Norm(NormMode mode){
+
+	IssmDouble norm;
 	int i;
 
@@ -210,5 +210,5 @@
 /*}}}*/
 /*FUNCTION SeqVec::Scale{{{*/
-void SeqVec::Scale(double scale_factor){
+void SeqVec::Scale(IssmDouble scale_factor){
 
 	int i;
@@ -218,9 +218,9 @@
 /*}}}*/
 /*FUNCTION SeqVec::Dot{{{*/
-double SeqVec::Dot(SeqVec* input){
-
-	int i;
-
-	double dot=0;
+IssmDouble SeqVec::Dot(SeqVec* input){
+
+	int i;
+
+	IssmDouble dot=0;
 	for(i=0;i<this->M;i++)dot+=this->vector[i]*input->vector[i];
 	return dot;
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h	(revision 12476)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h	(revision 12477)
@@ -1,4 +1,4 @@
 /*!\file:  SeqVec.h
- * \brief wrapper to our SeqVec object, which is just a wrapper to a double* 
+ * \brief wrapper to our SeqVec object, which is just a wrapper to a IssmDouble* 
  */ 
 
@@ -22,5 +22,5 @@
 	public:
 	
-		double* vector;
+		IssmDouble* vector;
 		int M;
 
@@ -28,5 +28,5 @@
 		SeqVec();
 		SeqVec(int M,bool fromlocalsize=false);
-		SeqVec(double* buffer, int M);
+		SeqVec(IssmDouble* buffer, int M);
 		~SeqVec();
 		/*}}}*/
@@ -34,19 +34,19 @@
 		void Echo(void);
 		void Assemble(void);
-		void SetValues(int ssize, int* list, double* values, InsMode mode);
-		void SetValue(int dof, double value, InsMode  mode);
-		void GetValue(double* pvalue, int dof);
+		void SetValues(int ssize, int* list, IssmDouble* values, InsMode mode);
+		void SetValue(int dof, IssmDouble value, InsMode  mode);
+		void GetValue(IssmDouble* pvalue, int dof);
 		void GetSize(int* pM);
 		void GetLocalSize(int* pM);
 		SeqVec* Duplicate(void);
-		void Set(double value);
-		void AXPY(SeqVec* X, double a);
-		void AYPX(SeqVec* X, double a);
-		double* ToMPISerial(void);
+		void Set(IssmDouble value);
+		void AXPY(SeqVec* X, IssmDouble a);
+		void AYPX(SeqVec* X, IssmDouble a);
+		IssmDouble* ToMPISerial(void);
 		void Copy(SeqVec* to);
-		double Norm(NormMode norm_type);
-		void Scale(double scale_factor);
+		IssmDouble Norm(NormMode norm_type);
+		void Scale(IssmDouble scale_factor);
 		void PointwiseDivide(SeqVec* x,SeqVec* y);
-		double Dot(SeqVec* vector);
+		IssmDouble Dot(SeqVec* vector);
 		/*}}}*/
 };
Index: /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 12476)
+++ /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 12477)
@@ -6,4 +6,6 @@
 #define _ISSM_TOOLKIT_H_
 
+#include "../../include/include.h"
+
 #include "./SeqMat.h"
 #include "./SeqVec.h"