Index: /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 12476)
@@ -16,5 +16,5 @@
 #include <float.h>
 
-void BrentSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs){
+void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs){
 
 	/* This routine is optimizing a given function using Brent's method
@@ -22,13 +22,13 @@
 
 	/*Intermediary*/
-	double si,gold,intervalgold,oldintervalgold;
-	double parab_num,parab_den;
-	double distance,cm_jump;
-	double fxmax,fxmin,fxbest;
-	double fx,fx1,fx2;
-	double xmax,xmin,xbest;
-	double x,x1,x2,xm;
-	double tol1,tol2,seps;
-	double tolerance=1.e-4;
+	IssmDouble si,gold,intervalgold,oldintervalgold;
+	IssmDouble parab_num,parab_den;
+	IssmDouble distance,cm_jump;
+	IssmDouble fxmax,fxmin,fxbest;
+	IssmDouble fx,fx1,fx2;
+	IssmDouble xmax,xmin,xbest;
+	IssmDouble x,x1,x2,xm;
+	IssmDouble tol1,tol2,seps;
+	IssmDouble tolerance=1.e-4;
 	int    maxiter,iter;
 	bool   loop=true,goldenflag;
@@ -58,5 +58,5 @@
 
 	/*initialize optimization variables*/
-	seps=sqrt(DBL_EPSILON);    //precision of a double
+	seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
 	distance=0.0;              //new_x=old_x + distance
 	gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
Index: /issm/trunk-jpl/src/c/shared/Numerics/IsInputConverged.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/IsInputConverged.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/IsInputConverged.cpp	(revision 12476)
@@ -14,17 +14,17 @@
 #include "../../objects/objects.h"
 
-void IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum){
+void IsInputConverged(IssmDouble* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum){
 
 	int i,j;
 
 	/*output: */
-	double eps;
+	IssmDouble eps;
 	
 	/*intermediary: */
-	double *newvalues     = NULL;
-	double *oldvalues     = NULL;
+	IssmDouble *newvalues     = NULL;
+	IssmDouble *oldvalues     = NULL;
 	int     num_values;
-	double  ndu        = 0;
-	double  nu         = 0;
+	IssmDouble  ndu        = 0;
+	IssmDouble  nu         = 0;
 
 	if(criterion_enum==RelativeEnum){
Index: /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp	(revision 12476)
@@ -16,15 +16,15 @@
 #include <float.h>
 
-void OptimalSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs){
+void OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs){
 
 	/* This routine is optimizing a given function*/
 
 	/*function values: */
-	double fx1,fx2,fxbest;
-	double x1,x2,xmin,xbest;
+	IssmDouble fx1,fx2,fxbest;
+	IssmDouble x1,x2,xmin,xbest;
 
 	/*tolerances: */
-	double seps;
-	double tolerance=1.e-4;
+	IssmDouble seps;
+	IssmDouble tolerance=1.e-4;
 	int    maxiter;
 
Index: /issm/trunk-jpl/src/c/shared/Numerics/OptionsFromAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/OptionsFromAnalysis.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/OptionsFromAnalysis.cpp	(revision 12476)
@@ -22,5 +22,5 @@
 	/*intermediary: */
 	int     dummy;
-	double* analyses=NULL;
+	IssmDouble* analyses=NULL;
 	char**  strings=NULL;
 	char*   string=NULL;
@@ -65,5 +65,5 @@
 
 	/*Free ressources*/
-	xDelete<double>(analyses);
+	xDelete<IssmDouble>(analyses);
 	for(i=0;i<numanalyses;i++){
 		string=strings[i];
Index: /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp	(revision 12476)
@@ -15,10 +15,10 @@
 #include "../../shared/shared.h"
 
-double UnitConversionScaleFactor(int type_enum);
+IssmDouble UnitConversionScaleFactor(int type_enum);
 /*}}}*/
 
-void UnitConversion(double* values, int numvalues,int direction_enum, int type_enum){
+void UnitConversion(IssmDouble* values, int numvalues,int direction_enum, int type_enum){
 
-	double scale;
+	IssmDouble scale;
 	int i;
 
@@ -33,5 +33,5 @@
 }
 
-double UnitConversion(double value, int direction_enum, int type_enum){
+IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum){
 
 	UnitConversion(&value,1,direction_enum,type_enum);
@@ -41,9 +41,9 @@
 
 
-double UnitConversionScaleFactor(int type_enum){
+IssmDouble UnitConversionScaleFactor(int type_enum){
 
-	double yts=365.0*24.0*3600.0;
+	IssmDouble yts=365.0*24.0*3600.0;
 	
-	double scale;
+	IssmDouble scale;
 	switch(type_enum){
 		case TimesteppingStartTimeEnum:              scale=1.0/yts;break; //yr
@@ -67,5 +67,5 @@
 		case SurfaceforcingsMassBalanceEnum:         scale=yts;break; //m/yr
 		case MisfitEnum:                             scale=pow(yts,2);break; //(m/yr)^2
-		case MassFluxEnum:                           scale=pow((double)10,-12)*yts;break; // (GigaTon/year)
+		case MassFluxEnum:                           scale=pow((IssmDouble)10,-12)*yts;break; // (GigaTon/year)
 		default: scale=1.0; break;
 	}
Index: /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp	(revision 12476)
@@ -5,9 +5,9 @@
 #include <math.h>
 
-void XZvectorsToCoordinateSystem(double* T,double* xzvectors){
+void XZvectorsToCoordinateSystem(IssmDouble* T,IssmDouble* xzvectors){
 
 	int		i,j;
-	double	x[3],y[3],z[3];
-	double	x_norm, y_norm, z_norm;
+	IssmDouble	x[3],y[3],z[3];
+	IssmDouble	x_norm, y_norm, z_norm;
 
 	for(i=0;i<6;i++){
Index: /issm/trunk-jpl/src/c/shared/Numerics/cross.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/cross.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/cross.cpp	(revision 12476)
@@ -9,5 +9,7 @@
 #endif
 
-void cross(double* result,double* vector1,double* vector2){
+#include "../../include/include.h"
+
+void cross(IssmDouble* result,IssmDouble* vector1,IssmDouble* vector2){
 
 	/*result,vector1 and vector2 are all assumed to be of size 3: */
Index: /issm/trunk-jpl/src/c/shared/Numerics/extrema.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/extrema.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/extrema.cpp	(revision 12476)
@@ -9,6 +9,7 @@
 #endif
 
+#include "../../include/include.h"
 
-double min(double a,double b){
+IssmDouble min(IssmDouble a,IssmDouble b){
 	if (a<b)return a;
 	else return b;
@@ -18,5 +19,5 @@
 	else return b;
 }
-double max(double a,double b){
+IssmDouble max(IssmDouble a,IssmDouble b){
 	if (a>b)return a;
 	else return b;
Index: /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp	(revision 12476)
@@ -2,5 +2,5 @@
 #ifdef _INTEL_WIN_
 
-int isnan(double x){
+int isnan(IssmDouble x){
 	if (x!=x)return 1;
 	else return 0;
Index: /issm/trunk-jpl/src/c/shared/Numerics/isnan.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/isnan.h	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/isnan.h	(revision 12476)
@@ -7,5 +7,5 @@
 
 #ifdef _INTEL_WIN_
-int isnan(double X);
+int isnan(IssmDouble X);
 #endif
 
Index: /issm/trunk-jpl/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 12475)
+++ /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 12476)
@@ -16,17 +16,17 @@
 struct OptPars;
 
-double min(double a,double b);
-double max(double a,double b);
+IssmDouble min(IssmDouble a,IssmDouble b);
+IssmDouble max(IssmDouble a,IssmDouble b);
 int    min(int a,int b);
 int    max(int a,int b);
-double OptFunc(double scalar, OptArgs* optargs);
-void   BrentSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
-void   OptimalSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
-void   cross(double* result,double* vector1,double* vector2);
-void   IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
-void   UnitConversion(double* values, int numvalues,int direction_enum, int type_enum);
-double UnitConversion(double value, int direction_enum, int type_enum);
+IssmDouble OptFunc(IssmDouble scalar, OptArgs* optargs);
+void   BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
+void   OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
+void   cross(IssmDouble* result,IssmDouble* vector1,IssmDouble* vector2);
+void   IsInputConverged(IssmDouble* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
+void   UnitConversion(IssmDouble* values, int numvalues,int direction_enum, int type_enum);
+IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum);
 char*  OptionsFromAnalysis(Parameters* parameters,int analysis_type);
-void   XZvectorsToCoordinateSystem(double* T,double* xzvectors);
+void   XZvectorsToCoordinateSystem(IssmDouble* T,IssmDouble* xzvectors);
 #ifdef _HAVE_PETSC_
 void   PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type);