Index: /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp	(revision 12529)
@@ -119,5 +119,5 @@
 			if (verbosity>1) _printLine_("   Merging Metric with hVertices...");
 			for (i=0;i<BTh.nbv;i++){
-				if (!isnan(bamgopts->hVertices[i])){
+				if (!xIsNan<IssmDouble>(bamgopts->hVertices[i])){
 					BTh[i].m=Metric((float)bamgopts->hVertices[i]);
 				}
@@ -129,5 +129,5 @@
 			if (verbosity>1) _printLine_("   Merging Metric with hminVertices...");
 			for (i=0;i<BTh.nbv;i++){
-				if (!isnan(bamgopts->hminVertices[i])){
+				if (!xIsNan<IssmDouble>(bamgopts->hminVertices[i])){
 					Metric M=BTh.vertices[i].m;
 					EigenMetric Vp(M/coef);
@@ -142,5 +142,5 @@
 			if (verbosity>1) _printLine_("   Merging Metric with hmaxVertices...");
 			for (i=0;i<BTh.nbv;i++){
-				if (!isnan(bamgopts->hmaxVertices[i])){
+				if (!xIsNan<IssmDouble>(bamgopts->hmaxVertices[i])){
 					Metric M=BTh.vertices[i].m;
 					EigenMetric Vp(M/coef);
Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 12529)
@@ -55,5 +55,5 @@
 	norm_inf=gradient->Norm(NORM_INF);
 	if(norm_inf<=0)    _error2_("||∂J/∂α||∞ = 0    gradient norm is zero");
-	if(isnan(norm_inf))_error2_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
+	if(xIsNan<IssmDouble>(norm_inf))_error2_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
 
 	/*Clean-up and assign output pointer*/
Index: /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 12529)
@@ -13,4 +13,5 @@
 #include "../../shared/shared.h"
 #include "../../include/include.h"
+#include "../../io/io.h"
 /*}}}*/
 
@@ -140,5 +141,5 @@
 
 		if(debug && my_thread==0)
-		 _pprintString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(i-i0)/double(i1-i0)*100<<"%");
+		 _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(i-i0)/double(i1-i0)*100<<"%");
 		x_grid=*(x_mesh+i);
 		y_grid=*(y_mesh+i);
@@ -179,5 +180,5 @@
 					return NULL; /*WARNING: no error because it would blow up the multithreading!*/
 			}
-			if(isnan(data_value)) data_value=default_value;
+			if(xIsNan<IssmDouble>(data_value)) data_value=default_value;
 		}
 		else{
Index: /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 12529)
@@ -113,5 +113,5 @@
 						data_value=data[i];
 					}
-					if (isnan(data_value)){
+					if (xIsNan<IssmDouble>(data_value)){
 						if(num_default_values==1) data_value=default_values[0];
 						else data_value=default_values[j];
Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 12529)
@@ -160,5 +160,5 @@
 						data_value=data_mesh[n];
 					}
-					if (isnan(data_value)) data_value=default_value;
+					if (xIsNan<IssmDouble>(data_value)) data_value=default_value;
 
 					/*insert value and go to the next point*/
Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp	(revision 12529)
@@ -124,5 +124,5 @@
 						data_value=data[i];
 					}
-					if (isnan(data_value)) data_value=default_value;
+					if (xIsNan<IssmDouble>(data_value)) data_value=default_value;
 
 					/*insert value and go to the next point*/
Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 12529)
@@ -55,5 +55,5 @@
 			if((iomodel->my_vertices[i])){
 
-				if (!isnan(IssmDoublevector[i])){
+				if (!xIsNan<IssmDouble>(IssmDoublevector[i])){
 
 					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type));
@@ -87,5 +87,5 @@
 				for(j=0;j<N;j++){
 					values[j]=IssmDoublevector[i*N+j];
-					if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+					if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
 				}
 
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 12529)
@@ -120,9 +120,9 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -135,9 +135,9 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -157,9 +157,9 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -172,13 +172,13 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvz[i])){
+						if (!xIsNan<IssmDouble>(spcvz[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -197,9 +197,9 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -212,13 +212,13 @@
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 						count++;
-						if (!isnan(spcvx[i])){
+						if (!xIsNan<IssmDouble>(spcvx[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvy[i])){
+						if (!xIsNan<IssmDouble>(spcvy[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
 						}
-						if (!isnan(spcvz[i])){
+						if (!xIsNan<IssmDouble>(spcvz[i])){
 							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 							count++;
@@ -229,5 +229,5 @@
 			/*Now add the regular spcs*/
 			else{
-				if (Mx==numberofvertices && !isnan(spcvx[i])){
+				if (Mx==numberofvertices && !xIsNan<IssmDouble>(spcvx[i])){
 					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 					count++;
@@ -239,5 +239,5 @@
 					for(j=0;j<Nx;j++){
 						values[j]=spcvx[i*Nx+j]/yts;
-						if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+						if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
 					}
 
@@ -253,5 +253,5 @@
 				}
 
-				if (My==numberofvertices && !isnan(spcvy[i])){
+				if (My==numberofvertices && !xIsNan<IssmDouble>(spcvy[i])){
 					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
 					count++;
@@ -263,5 +263,5 @@
 					for(j=0;j<Ny;j++){
 						values[j]=spcvy[i*Ny+j]/yts;
-						if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+						if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
 					}
 					if(spcpresent){
@@ -277,5 +277,5 @@
 
 				if ((int)vertices_type[i]==StokesApproximationEnum ||  ((int)vertices_type[i]==NoneApproximationEnum)){
-					if (Mz==numberofvertices && !isnan(spcvz[i])){
+					if (Mz==numberofvertices && !xIsNan<IssmDouble>(spcvz[i])){
 						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
 						count++;
@@ -287,5 +287,5 @@
 						for(j=0;j<Nz;j++){
 							values[j]=spcvz[i*Nz+j]/yts;
-							if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+							if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
 						}
 						if(spcpresent){
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 12529)
@@ -60,10 +60,10 @@
 			}
 			else{
-				if (!isnan(iomodel->Data(DiagnosticSpcvxEnum)[i])){
+				if (!xIsNan<IssmDouble>(iomodel->Data(DiagnosticSpcvxEnum)[i])){
 					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(DiagnosticSpcvxEnum)[i]/yts,DiagnosticHutterAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
 					count++;
 				}
 
-				if (!isnan(iomodel->Data(DiagnosticSpcvyEnum)[i])){
+				if (!xIsNan<IssmDouble>(iomodel->Data(DiagnosticSpcvyEnum)[i])){
 					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,iomodel->Data(DiagnosticSpcvyEnum)[i]/yts,DiagnosticHutterAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
 					count++;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 12529)
@@ -56,5 +56,5 @@
 				count++;
 			}
-			else if (!isnan(iomodel->Data(DiagnosticSpcvzEnum)[i])){
+			else if (!xIsNan<IssmDouble>(iomodel->Data(DiagnosticSpcvzEnum)[i])){
 				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,
 								iomodel->Data(DiagnosticSpcvzEnum)[i]/yts,DiagnosticVertAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp	(revision 12529)
@@ -55,5 +55,5 @@
 		if((iomodel->my_vertices[i])){
 
-			if (!isnan(spctemperature[i])){
+			if (!xIsNan<IssmDouble>(spctemperature[i])){
 
 				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spctemperature[i]-referencetemperature),EnthalpyAnalysisEnum));
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp	(revision 12529)
@@ -125,5 +125,5 @@
 		 * we must clone the nodes on this partition so that the loads (Numericalflux)
 		 * will have access to their properties (dofs,...)*/
-		if(my_elements[(int)e1] && !isnan(e2) && !my_elements[(int)e2]){ 
+		if(my_elements[(int)e1] && !xIsNan<IssmDouble>(e2) && !my_elements[(int)e2]){ 
 
 			/*1: Get vertices ids*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 12529)
@@ -43,5 +43,5 @@
 		/*keep only this partition's nodes:*/
 		if((iomodel->my_vertices[i]==1)){
-			if (isnan(iomodel->Data(ThermalSpctemperatureEnum)[i])){ //No penalty applied on spc nodes!
+			if (xIsNan<IssmDouble>(iomodel->Data(ThermalSpctemperatureEnum)[i])){ //No penalty applied on spc nodes!
 				loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum));
 			}
Index: /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp	(revision 12529)
@@ -183,5 +183,5 @@
 			if(verbose>5) _printLine_("      processing hVertices");
 			for (i=0;i< nbv;i++){
-				if (!isnan(bamgopts->hVertices[i])){
+				if (!xIsNan<IssmDouble>(bamgopts->hVertices[i])){
 					vertices[i].m=Metric((double)bamgopts->hVertices[i]);
 				}
Index: /issm/trunk-jpl/src/c/objects/Constraints/SpcDynamic.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Constraints/SpcDynamic.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Constraints/SpcDynamic.cpp	(revision 12529)
@@ -124,5 +124,5 @@
 IssmDouble SpcDynamic::GetValue(){
 	_assert_(this->isset);
-	_assert_(!isnan(value));
+	_assert_(!xIsNan<IssmDouble>(value));
 	return value;
 }
Index: /issm/trunk-jpl/src/c/objects/Constraints/SpcStatic.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Constraints/SpcStatic.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Constraints/SpcStatic.cpp	(revision 12529)
@@ -124,5 +124,5 @@
 /*FUNCTION SpcStatic::GetValue {{{*/
 IssmDouble SpcStatic::GetValue(){
-	_assert_(!isnan(value));
+	_assert_(!xIsNan<IssmDouble>(value));
 	return value;
 }
Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 12529)
@@ -66,7 +66,7 @@
 
 	/*Build neighbors list*/
-	if (isnan(iomodel->Data(MeshUpperelementsEnum)[index])) penta_elements_ids[1]=this->id; //upper penta is the same penta
+	if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index])) penta_elements_ids[1]=this->id; //upper penta is the same penta
 	else                                    penta_elements_ids[1]=(int)(iomodel->Data(MeshUpperelementsEnum)[index]);
-	if (isnan(iomodel->Data(MeshLowerelementsEnum)[index])) penta_elements_ids[0]=this->id; //lower penta is the same penta
+	if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index])) penta_elements_ids[0]=this->id; //lower penta is the same penta
 	else                                    penta_elements_ids[0]=(int)(iomodel->Data(MeshLowerelementsEnum)[index]);
 	this->InitHookNeighbors(penta_elements_ids);
@@ -1742,5 +1742,5 @@
 	for(i=0;i<numdof2d;i++){
 		newthickness[i]=solution[doflist[i]];
-		if(isnan(newthickness[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(newthickness[i])) _error2_("NaN found in solution vector");
 		/*Constrain thickness to be at least 1m*/
 		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
@@ -1812,5 +1812,5 @@
 	for(int i=0;i<numdof;i++){
 		values[i]=solution[doflist[i]];
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -1842,5 +1842,5 @@
 		values[i]         =solution[doflist[i]];
 		values[i+numdof2d]=values[i];
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -4328,5 +4328,5 @@
 
 		/*Check solution*/
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 		//if(values[i]<0)      _printLine_("temperature < 0°K found in solution vector");
 		//if(values[i]>275)    _printLine_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)");
@@ -4398,5 +4398,5 @@
 
 		/*Check solution*/
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -4943,5 +4943,5 @@
 		/*Add gradje_g_gaussian vector to gradje_g: */
 		for(i=0;i<NUMVERTICES;i++){
-			_assert_(!isnan(grade_g[i]));
+			_assert_(!xIsNan<IssmDouble>(grade_g[i]));
 			grade_g[i]+=grade_g_gaussian[i];
 		}
@@ -5163,8 +5163,8 @@
 
 		/*Check solution*/
-		if(isnan(lambdax[i])) _error2_("NaN found in solution vector");
-		if(isnan(lambday[i])) _error2_("NaN found in solution vector");
-		if(isnan(lambdaz[i])) _error2_("NaN found in solution vector");
-		if(isnan(lambdap[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambdax[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambday[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambdaz[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambdap[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -5202,6 +5202,6 @@
 
 		/*Check solution*/
-		if(isnan(lambdax[i]))       _error2_("NaN found in solution vector");
-		if(isnan(lambday[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambdax[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambday[i]))       _error2_("NaN found in solution vector");
 	}
 
@@ -8128,6 +8128,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8226,6 +8226,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8314,8 +8314,8 @@
 
 		/*Check solution*/
-		if(isnan(vx[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vy[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vzstokes[i])) _error2_("NaN found in solution vector");
-		if(isnan(pressure[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzstokes[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(pressure[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8393,6 +8393,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8482,8 +8482,8 @@
 
 		/*Check solution*/
-		if(isnan(vx[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vy[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vzstokes[i])) _error2_("NaN found in solution vector");
-		if(isnan(pressure[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzstokes[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(pressure[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8558,6 +8558,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8627,5 +8627,5 @@
 
 		/*Check solution*/
-		if(isnan(vz[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vz[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -8725,8 +8725,8 @@
 
 		/*Check solution*/
-		if(isnan(vx[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vy[i]))       _error2_("NaN found in solution vector");
-		if(isnan(vz[i]))       _error2_("NaN found in solution vector");
-		if(isnan(pressure[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vz[i]))       _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(pressure[i])) _error2_("NaN found in solution vector");
 	}
 
Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp	(revision 12529)
@@ -1523,5 +1523,5 @@
 	for(int i=0;i<numdof;i++){
 		values[i]=solution[doflist[i]];
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -1556,5 +1556,5 @@
 	for(i=0;i<numdof;i++){
 		newthickness[i]=solution[doflist[i]];
-		if(isnan(newthickness[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(newthickness[i])) _error2_("NaN found in solution vector");
 		/*Constrain thickness to be at least 1m*/
 		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
@@ -3325,6 +3325,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -3385,6 +3385,6 @@
 
 		/*Check solution*/
-		if(isnan(vx[i])) _error2_("NaN found in solution vector");
-		if(isnan(vy[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vx[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -3731,5 +3731,5 @@
 		/*Add gradje_g_gaussian vector to gradje_g: */
 		for(i=0;i<NUMVERTICES;i++){
-			_assert_(!isnan(grade_g[i]));
+			_assert_(!xIsNan<IssmDouble>(grade_g[i]));
 			grade_g[i]+=grade_g_gaussian[i];
 		}
@@ -3792,5 +3792,5 @@
 		for (i=0;i<NUMVERTICES;i++){
 			grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
-			_assert_(!isnan(grade_g[i]));
+			_assert_(!xIsNan<IssmDouble>(grade_g[i]));
 		}
 	}
@@ -4947,6 +4947,6 @@
 
 		/*Check solution*/
-		if(isnan(lambdax[i])) _error2_("NaN found in solution vector");
-		if(isnan(lambday[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambdax[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambday[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -4978,5 +4978,5 @@
 	for(i=0;i<numdof;i++){
 		lambda[i]=values[i];
-		if(isnan(lambda[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(lambda[i])) _error2_("NaN found in solution vector");
 	}
 
@@ -5314,5 +5314,5 @@
 	for(i=0;i<numdof;i++){
 		values[i]=solution[doflist[i]];
-		if(isnan(values[i])) _error2_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(values[i])) _error2_("NaN found in solution vector");
 		if (values[i]<pow((IssmDouble)10,(IssmDouble)-10))values[i]=pow((IssmDouble)10,(IssmDouble)-10); //correcting the water column to positive values
  
Index: /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Inputs/BoolInput.cpp	(revision 12529)
@@ -169,6 +169,6 @@
 void BoolInput::Constrain(IssmPDouble cm_min, IssmPDouble cm_max){
 
-	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
-	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+	if(!xIsNan<IssmDouble>(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!xIsNan<IssmDouble>(cm_max)) if (this->value>cm_max)this->value=cm_max;
 
 }
Index: /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.cpp	(revision 12529)
@@ -225,6 +225,6 @@
 void DoubleInput::Constrain(IssmPDouble cm_min, IssmPDouble cm_max){
 
-	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
-	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+	if(!xIsNan<IssmDouble>(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!xIsNan<IssmDouble>(cm_max)) if (this->value>cm_max)this->value=cm_max;
 
 }
Index: /issm/trunk-jpl/src/c/objects/Inputs/IntInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Inputs/IntInput.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Inputs/IntInput.cpp	(revision 12529)
@@ -174,6 +174,6 @@
 void IntInput::Constrain(IssmPDouble cm_min, IssmPDouble cm_max){
 
-	if(!isnan(cm_min)) if (this->value<cm_min)this->value=(int)cm_min;
-	if(!isnan(cm_max)) if (this->value>cm_max)this->value=(int)cm_max;
+	if(!xIsNan<IssmDouble>(cm_min)) if (this->value<cm_min)this->value=(int)cm_min;
+	if(!xIsNan<IssmDouble>(cm_max)) if (this->value>cm_max)this->value=(int)cm_max;
 
 }
Index: /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.cpp	(revision 12529)
@@ -471,6 +471,6 @@
 	const int numnodes=6;
 		
-	if(!isnan(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
-	if(!isnan(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
+	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
+	if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
 
 }
Index: /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.cpp	(revision 12529)
@@ -334,6 +334,6 @@
 	const int numnodes=3;
 		
-	if(!isnan(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
-	if(!isnan(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
+	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
+	if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
 
 }
Index: /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp	(revision 12529)
@@ -113,5 +113,5 @@
 
 	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
-	_assert_(!isnan(alpha2));
+	_assert_(!xIsNan<IssmDouble>(alpha2));
 
 	/*Assign output pointers:*/
@@ -177,5 +177,5 @@
 
 	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
-	_assert_(!isnan(alpha2));
+	_assert_(!xIsNan<IssmDouble>(alpha2));
 
 	/*Assign output pointers:*/
@@ -243,5 +243,5 @@
 	if(vmag==0 && (s-1)<0) _error2_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");
 
-	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!isnan(alpha_complement));
+	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));
 
 	/*Assign output pointers:*/
@@ -309,5 +309,5 @@
 	if(vmag==0 && (s-1)<0) _error2_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");
 
-	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!isnan(alpha_complement));
+	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));
 
 	/*Assign output pointers:*/
Index: /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.cpp	(revision 12529)
@@ -60,5 +60,5 @@
 
 	/*First, see wether this is an internal or boundary edge (if e2=NaN)*/
-	if (isnan((IssmDouble)iomodel->Data(MeshEdgesEnum)[4*i+3])){ //edges are [node1 node2 elem1 elem2]
+	if (xIsNan<IssmDouble>((IssmDouble)iomodel->Data(MeshEdgesEnum)[4*i+3])){ //edges are [node1 node2 elem1 elem2]
 		/* Boundary edge, only one element */
 		e1=(int)iomodel->Data(MeshEdgesEnum)[4*i+2];
@@ -751,5 +751,5 @@
 		vyaverage_input->GetInputValue(&vy,gauss);
 		spcthickness_input->GetInputValue(&thickness,gauss);
-		if(isnan(thickness)) _error2_("Cannot weakly apply constraint because NaN was provided");
+		if(xIsNan<IssmDouble>(thickness)) _error2_("Cannot weakly apply constraint because NaN was provided");
 
 		UdotN=vx*normal[0]+vy*normal[1];
Index: /issm/trunk-jpl/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Node.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Node.cpp	(revision 12529)
@@ -524,5 +524,5 @@
 			if(this->indexing.s_set[i]){
 				values[count]=this->indexing.svalues[i];
-				_assert_(!isnan(values[count]));
+				_assert_(!xIsNan<IssmDouble>(values[count]));
 				count++;
 			}
Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 12529)
@@ -343,5 +343,5 @@
 	for (int i=0;i<this->nrows;i++){
 		for(int j=0;j<this->ncols;j++){
-			if (isnan(this->values[i*this->ncols+j])) _error2_("NaN found in Element Matrix");
+			if (xIsNan<IssmDouble>(this->values[i*this->ncols+j])) _error2_("NaN found in Element Matrix");
 			if (fabs(this->values[i*this->ncols+j])>1.e+50) _error2_("Element Matrix values exceeds 1.e+50");
 		}
Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp	(revision 12529)
@@ -210,5 +210,5 @@
 #ifdef _ISSM_DEBUG_ 
 	for (int i=0;i<this->nrows;i++){
-		if (isnan(this->values[i])) _error2_("NaN found in Element Vector");
+		if (xIsNan<IssmDouble>(this->values[i])) _error2_("NaN found in Element Vector");
 		if (fabs( this->values[i])>1.e+50) _error2_("Element Vector values exceeds 1.e+50");
 	}
Index: /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 12529)
@@ -44,5 +44,5 @@
 	iter=0;
 	fxmin = (*f)(xmin,optargs);
-	if (isnan(fxmin)) _error2_("Function evaluation returned NaN");
+	if (xIsNan<IssmDouble>(fxmin)) _error2_("Function evaluation returned NaN");
 	cout<<setprecision(5);
 	if(VerboseControl()) _pprintLine_("");
@@ -51,9 +51,9 @@
 	if(VerboseControl()) _pprintLine_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary");
 	fxmax = (*f)(xmax,optargs);
-	if (isnan(fxmax)) _error2_("Function evaluation returned NaN");
+	if (xIsNan<IssmDouble>(fxmax)) _error2_("Function evaluation returned NaN");
 	if(VerboseControl()) _pprintLine_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary");
 
 	/*test if jump option activated and xmin==0*/
-	if (!isnan(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
+	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
 		*psearch_scalar=xmax;
 		*pJ=fxmax;
@@ -75,5 +75,5 @@
 	/*2: call the function to be evaluated*/
 	fxbest = (*f)(x,optargs);
-	if(isnan(fxbest)) _error2_("Function evaluation returned NaN");
+	if(xIsNan<IssmDouble>(fxbest)) _error2_("Function evaluation returned NaN");
 	iter=iter+1;
 
@@ -90,5 +90,5 @@
 	if(VerboseControl())
 	 _pprintLine_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial");
-	if (!isnan(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
+	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
 		if(VerboseControl()) _pprintLine_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump);
 		loop=false;
@@ -159,5 +159,5 @@
 		//evaluate function on x
 		fx = (*f)(x,optargs);
-		if(isnan(fx)) _error2_("Function evaluation returned NaN");
+		if(xIsNan<IssmDouble>(fx)) _error2_("Function evaluation returned NaN");
 		iter=iter+1;
 
@@ -197,5 +197,5 @@
 			loop=false;
 		}
-		else if (!isnan(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
+		else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
 			if(VerboseControl()) _pprintLine_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump);
 			loop=false;
Index: /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp	(revision 12529)
@@ -41,5 +41,5 @@
 	//get the value of the function at the first boundary
 	fx1= (*f)(x1,optargs);
-	if (isnan(fx1)) _error2_("Function evaluation returned NaN");
+	if (xIsNan<IssmDouble>(fx1)) _error2_("Function evaluation returned NaN");
 	cout<<setprecision(5);
 	if(VerboseControl()) _pprintLine_("");
@@ -57,5 +57,5 @@
 		iter++;
 		fx2 = (*f)(x2,optargs);
-		if (isnan(fx2)) _error2_("Function evaluation returned NaN");
+		if (xIsNan<IssmDouble>(fx2)) _error2_("Function evaluation returned NaN");
 		if(VerboseControl())
 		 _pprintLine_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x2<<"  "<<setw(12)<<fx2<<"  "<<(fabs(x2-x1)>fabs(fx2-fx1)?fabs(fx2-fx1):fabs(x2-x1)));
Index: /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp	(revision 12529)
@@ -12,5 +12,5 @@
 
 	for(i=0;i<6;i++){
-		if(isnan(xzvectors[i])){
+		if(xIsNan<IssmDouble>(xzvectors[i])){
 			/*At least one NaN found: default to Id*/
 			T[0*3+0] = 1.0;	T[0*3+1] = 0.0;	T[0*3+2] = 0.0;
Index: /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/shared/Numerics/isnan.cpp	(revision 12529)
@@ -2,5 +2,5 @@
 #ifdef _INTEL_WIN_
 
-int isnan(IssmPDouble x){
+int xIsNan<IssmDouble>(IssmPDouble x){
 	return (x!=x)?1:0;
 }
Index: /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/solutions/controlconvergence.cpp	(revision 12529)
@@ -23,5 +23,5 @@
 
 	/*Has convergence been reached?*/
-	if (!isnan(tol_cm) && J<tol_cm){
+	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
 		converged=true;
 		if(VerboseConvergence()) _pprintString_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
Index: /issm/trunk-jpl/src/c/solutions/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/convergence.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/solutions/convergence.cpp	(revision 12529)
@@ -70,5 +70,5 @@
 	nF=pf->Norm(NORM_TWO);
 	res=nKUoldF/nF;
-	if (isnan(res)){
+	if (xIsNan<IssmDouble>(res)){
 		_pprintLine_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f");
 		_error2_("mechanical equilibrium convergence criterion is NaN!");
@@ -90,5 +90,5 @@
 
 	/*Relative criterion (optional)*/
-	if (!isnan(eps_rel) || (VerboseConvergence())){
+	if (!xIsNan<IssmDouble>(eps_rel) || (VerboseConvergence())){
 
 		//compute norm(du)/norm(u)
@@ -96,5 +96,5 @@
 		ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO);
 
-		if (isnan(ndu) || isnan(nu)) _error2_("convergence criterion is NaN!");
+		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error2_("convergence criterion is NaN!");
 
 		//clean up
@@ -102,5 +102,5 @@
 
 		//print
-		if (!isnan(eps_rel)){
+		if (!xIsNan<IssmDouble>(eps_rel)){
 			if((ndu/nu)<eps_rel){
 				if(VerboseConvergence()) _pprintLine_("   Convergence criterion: norm(du)/norm(u)" << "50s" << ndu/nu*100 << " < " << eps_rel*100 << " %");
@@ -116,10 +116,10 @@
 
 	/*Absolute criterion (Optional) = max(du)*/
-	if (!isnan(eps_abs) || (VerboseConvergence())){
+	if (!xIsNan<IssmDouble>(eps_abs) || (VerboseConvergence())){
 
 		//compute max(du)
 		duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
 		ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF);
-		if (isnan(ndu) || isnan(nu)) _error2_("convergence criterion is NaN!");
+		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error2_("convergence criterion is NaN!");
 
 		//clean up
@@ -127,5 +127,5 @@
 
 		//print
-		if (!isnan(eps_abs)){
+		if (!xIsNan<IssmDouble>(eps_abs)){
 			if ((nduinf*yts)<eps_abs){
 				if(VerboseConvergence()) _pprintLine_("   Convergence criterion: max(du)" << "50s" << nduinf*yts << " < " << eps_abs << " m/yr");
Index: /issm/trunk-jpl/src/c/solutions/gradient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/gradient_core.cpp	(revision 12528)
+++ /issm/trunk-jpl/src/c/solutions/gradient_core.cpp	(revision 12529)
@@ -38,5 +38,5 @@
 	norm_inf=new_gradient->Norm(NORM_INF);
 	if(norm_inf<=0)    _error2_("||∂J/∂α||∞ = 0    gradient norm is zero");
-	if(isnan(norm_inf))_error2_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
+	if(xIsNan<IssmDouble>(norm_inf))_error2_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
 
 	/*plug back into inputs: */