Index: /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.cpp	(revision 12423)
+++ /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.cpp	(revision 12424)
@@ -42,6 +42,6 @@
 
 	/*Intermediaties*/
-	double *seg_coords  = NULL;
-	double *seg_weights = NULL;
+	IssmPDouble *seg_coords  = NULL;
+	IssmPDouble *seg_weights = NULL;
 	int     i,index3;
 
@@ -51,8 +51,8 @@
 
 	/*Allocate GaussTria fields*/
-	coords1=(double*)xmalloc(numgauss*sizeof(double));
-	coords2=(double*)xmalloc(numgauss*sizeof(double));
-	coords3=(double*)xmalloc(numgauss*sizeof(double));
-	weights=(double*)xmalloc(numgauss*sizeof(double));
+	coords1=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
+	coords2=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
+	coords3=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
+	weights=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
 
 	/*Reverse index1 and 2 if necessary*/
@@ -71,5 +71,5 @@
 	else if (index1==0 && index2==2){
 		for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
-		for(i=0;i<numgauss;i++) coords2[i]=0;
+		for(i=0;i<numgauss;i++) coords2[i]= 0 ;
 		for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
 		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
@@ -97,8 +97,8 @@
 /*FUNCTION GaussTria::~GaussTria(){{{*/
 GaussTria::~GaussTria(){
-	xfree((void**)&weights);
-	xfree((void**)&coords1);
-	xfree((void**)&coords2);
-	xfree((void**)&coords3);
+	xDelete<IssmDouble>(weights);
+	xDelete<IssmDouble>(coords1);
+	xDelete<IssmDouble>(coords2);
+	xDelete<IssmDouble>(coords3);
 }
 /*}}}*/
@@ -199,9 +199,9 @@
 /*}}}*/
 /*FUNCTION GaussTria::GaussFromCoords{{{*/
-void GaussTria::GaussFromCoords(double x,double y,double* xyz_list){
+void GaussTria::GaussFromCoords(IssmDouble x,IssmDouble y,IssmDouble* xyz_list){
 
 	/*Intermediaries*/
-	double    area = 0;
-	double    x1,y1,x2,y2,x3,y3;
+	IssmDouble    area = 0;
+	IssmDouble    x1,y1,x2,y2,x3,y3;
 
 	/*in debugging mode: check that the default constructor has been called*/
Index: /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.h	(revision 12423)
+++ /issm/trunk-jpl/src/c/objects/Gauss/GaussTria.h	(revision 12424)
@@ -15,14 +15,14 @@
 	private:
 		int numgauss;
-		double* weights;
-		double* coords1;
-		double* coords2;
-		double* coords3;
+		IssmDouble* weights;
+		IssmDouble* coords1;
+		IssmDouble* coords2;
+		IssmDouble* coords3;
 
 	public:
-		double weight;
-		double coord1;
-		double coord2;
-		double coord3;
+		IssmDouble weight;
+		IssmDouble coord1;
+		IssmDouble coord2;
+		IssmDouble coord3;
 		
 	public:
@@ -38,5 +38,5 @@
 		int  end(void);
 		void Echo(void);
-		void GaussFromCoords(double x1,double y1,double* xyz_list);
+		void GaussFromCoords(IssmDouble x1,IssmDouble y1,IssmDouble* xyz_list);
 		void GaussPoint(int ig);
 		void GaussVertex(int iv);
Index: /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp	(revision 12423)
+++ /issm/trunk-jpl/src/c/objects/Loads/Friction.cpp	(revision 12424)
@@ -55,6 +55,6 @@
 }
 /*}}}*/
-/*FUNCTION Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
-void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
+/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
+void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
 
 	/*This routine calculates the basal friction coefficient 
@@ -62,12 +62,12 @@
 
 	/*diverse: */
-	double  r,s;
-	double  drag_p, drag_q;
-	double  gravity,rho_ice,rho_water;
-	double  Neff;
-	double  thickness,bed;
-	double  vx,vy,vz,vmag;
-	double  drag_coefficient;
-	double  alpha2;
+	IssmDouble  r,s;
+	IssmDouble  drag_p, drag_q;
+	IssmDouble  gravity,rho_ice,rho_water;
+	IssmDouble  Neff;
+	IssmDouble  thickness,bed;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_coefficient;
+	IssmDouble  alpha2;
 
 	/*Recover parameters: */
@@ -119,6 +119,6 @@
 }
 /*}}}*/
-/*FUNCTION Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
-void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
+/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
+void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
 
 	/*This routine calculates the basal friction coefficient 
@@ -126,12 +126,12 @@
 
 	/*diverse: */
-	double  r,s;
-	double  drag_p, drag_q;
-	double  gravity,rho_ice,rho_water;
-	double  Neff;
-	double  thickness,bed;
-	double  vx,vy,vz,vmag;
-	double  drag_coefficient;
-	double  alpha2;
+	IssmDouble  r,s;
+	IssmDouble  drag_p, drag_q;
+	IssmDouble  gravity,rho_ice,rho_water;
+	IssmDouble  Neff;
+	IssmDouble  thickness,bed;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_coefficient;
+	IssmDouble  alpha2;
 
 	/*Recover parameters: */
@@ -183,6 +183,6 @@
 }
 /*}}}*/
-/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
-void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
+/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
+void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
 
 	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
@@ -192,12 +192,12 @@
 	/*diverse: */
 	int     i;
-	double  r,s;
-	double  vx,vy,vz,vmag;
-	double  drag_p,drag_q;
-	double  Neff;
-	double  drag_coefficient;
-	double  bed,thickness;
-	double  gravity,rho_ice,rho_water;
-	double  alpha_complement;
+	IssmDouble  r,s;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_p,drag_q;
+	IssmDouble  Neff;
+	IssmDouble  drag_coefficient;
+	IssmDouble  bed,thickness;
+	IssmDouble  gravity,rho_ice,rho_water;
+	IssmDouble  alpha_complement;
 
 	/*Recover parameters: */
@@ -249,6 +249,6 @@
 }
 /*}}}*/
-/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
-void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
+/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
+void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
 
 	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
@@ -258,12 +258,12 @@
 	/*diverse: */
 	int     i;
-	double  r,s;
-	double  vx,vy,vz,vmag;
-	double  drag_p,drag_q;
-	double  Neff;
-	double  drag_coefficient;
-	double  bed,thickness;
-	double  gravity,rho_ice,rho_water;
-	double  alpha_complement;
+	IssmDouble  r,s;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_p,drag_q;
+	IssmDouble  Neff;
+	IssmDouble  drag_coefficient;
+	IssmDouble  bed,thickness;
+	IssmDouble  gravity,rho_ice,rho_water;
+	IssmDouble  alpha_complement;
 
 	/*Recover parameters: */
@@ -316,5 +316,5 @@
 /*}}}*/
 /*FUNCTION Friction::GetInputValue{{{*/
-void Friction::GetInputValue(double* pvalue,GaussTria* gauss,int enum_type){
+void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){
 
 	Input* input=inputs->GetInput(enum_type);
@@ -325,5 +325,5 @@
 /*}}}*/
 /*FUNCTION Friction::GetInputValue{{{*/
-void Friction::GetInputValue(double* pvalue,GaussPenta* gauss,int enum_type){
+void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){
 
 	Input* input=inputs->GetInput(enum_type);
Index: /issm/trunk-jpl/src/c/objects/Loads/Friction.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Friction.h	(revision 12423)
+++ /issm/trunk-jpl/src/c/objects/Loads/Friction.h	(revision 12424)
@@ -27,10 +27,10 @@
 	
 		void  Echo(void);
-		void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
-		void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
-		void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
-		void  GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
-		void  GetInputValue(double* pvalue,GaussTria* gauss,int enum_type);
-		void  GetInputValue(double* pvalue,GaussPenta* gauss,int enum_type);
+		void  GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetAlphaComplement(IssmDouble* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetAlphaComplement(IssmDouble* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type);
+		void  GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type);
 
 };
Index: /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp	(revision 12423)
+++ /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp	(revision 12424)
@@ -1,7 +1,7 @@
 /*  Gauss point structures and prototypes  */
 
+#include "../../include/include.h"
 #include "./GaussPoints.h"
 #include "../Alloc/alloc.h"
-#include "../../include/include.h"
 #include "../../io/io.h"
 #include "../Exceptions/exceptions.h"
@@ -11,5 +11,5 @@
 /*General Gauss points*/
 /*FUNCTION GaussLegendreLinear {{{*/
-void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus){
+void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus){
 	/* Gauss-Legendre quadrature points.
 
@@ -25,33 +25,33 @@
 	/*Intermediaries*/
 	int i;
-	double *alpha,*beta;
+	IssmPDouble *alpha,*beta;
 
 	/*p= 1, npoint= 1*/
-	static double wgt1[]={2.000000000000000};
-	static double xi1[]={0.000000000000000};
+	static IssmPDouble wgt1[]={2.000000000000000};
+	static IssmPDouble xi1[]={0.000000000000000};
 
 	/*p= 3, npoint= 2*/
-	static double wgt2[]={1.000000000000000, 1.000000000000000};
-	static double xi2[]={-0.577350269189626, 0.577350269189626};
+	static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
+	static IssmPDouble xi2[]={-0.577350269189626, 0.577350269189626};
 
 	/*p= 5, npoint= 3*/
-	static double wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};
-	static double xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};
+	static IssmPDouble wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};
+	static IssmPDouble xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};
 
 	/*p= 7, npoint= 4*/
-	static double wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
-	static double xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};
-
-	static double* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 };
-	static double* xip [MAX_LINE_GAUS_PTS]={xi1  ,xi2  ,xi3  ,xi4  };
-
-	static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof(double),
-		sizeof(wgt2 )/sizeof(double),
-		sizeof(wgt3 )/sizeof(double),
-		sizeof(wgt4 )/sizeof(double)};
+	static IssmPDouble wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
+	static IssmPDouble xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};
+
+	static IssmPDouble* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 };
+	static IssmPDouble* xip [MAX_LINE_GAUS_PTS]={xi1  ,xi2  ,xi3  ,xi4  };
+
+	static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble),
+		sizeof(wgt2 )/sizeof(IssmPDouble),
+		sizeof(wgt3 )/sizeof(IssmPDouble),
+		sizeof(wgt4 )/sizeof(IssmPDouble)};
 
 	//	_printf_(true,"Gauss-Legendre recurrence coefficients ngaus=%d\n",ngaus);
-	*pxgaus = (double *) xmalloc(ngaus*sizeof(double));
-	*pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
+	*pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
+	*pxwgt  = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
 
 	/*  check to see if Gauss points need to be calculated  */
@@ -70,6 +70,6 @@
 
 		/*  calculate the Gauss points using recurrence relations  */
-		alpha=(double *) xmalloc(ngaus*sizeof(double));
-		beta =(double *) xmalloc(ngaus*sizeof(double));
+		alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
+		beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
 
 		/*  calculate the Legendre recurrence coefficients  */
@@ -89,5 +89,5 @@
 }/*}}}*/
 /*FUNCTION GaussLegendreTria{{{*/
-void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
+void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord ) {
 	/*Gauss quadrature points for the triangle.
 
@@ -98,70 +98,70 @@
 	/*Intermediaries*/
 	int i,j,ipt,nigaus;
-	double xi,eta;
-	double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
+	IssmPDouble xi,eta;
+	IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
 
 	/*Hardcoded Gauss points declaration*/
 	/*p= 1, npoint= 1{{{*/
-	static double wgt1[]={
+	static IssmPDouble wgt1[]={
 		1.732050807568877};
-	static double l11[]={
+	static IssmPDouble l11[]={
 		0.333333333333333};
-	static double l21[]={
+	static IssmPDouble l21[]={
 		0.333333333333333};
-	static double l31[]={
+	static IssmPDouble l31[]={
 		0.333333333333333};
 	/*}}}*/
 	/*p= 2, npoint= 3  {{{*/
-	static double wgt2[]={
+	static IssmPDouble wgt2[]={
 		0.577350269189625, 0.577350269189625, 0.577350269189625};
-	static double l12[]={
+	static IssmPDouble l12[]={
 		0.666666666666667, 0.166666666666667, 0.166666666666667};
-	static double l22[]={
+	static IssmPDouble l22[]={
 		0.166666666666667, 0.666666666666667, 0.166666666666667};
-	static double l32[]={
+	static IssmPDouble l32[]={
 		0.166666666666667, 0.166666666666667, 0.666666666666667};
 	/*}}}*/
 	/*p= 3, npoint= 4  {{{*/
-	static double wgt3[]={
+	static IssmPDouble wgt3[]={
 		-0.974278579257493, 0.902109795608790, 0.902109795608790, 
 		0.902109795608790};
-	static double l13[]={
+	static IssmPDouble l13[]={
 		0.333333333333333, 0.600000000000000, 0.200000000000000, 
 		0.200000000000000};
-	static double l23[]={
+	static IssmPDouble l23[]={
 		0.333333333333333, 0.200000000000000, 0.600000000000000, 
 		0.200000000000000};
-	static double l33[]={
+	static IssmPDouble l33[]={
 		0.333333333333333, 0.200000000000000, 0.200000000000000, 
 		0.600000000000000};
 	/*}}}*/
 	/*p= 4, npoint= 6  {{{*/
-	static double wgt4[]={
+	static IssmPDouble wgt4[]={
 		0.386908262797819, 0.386908262797819, 0.386908262797819, 
 		0.190442006391807, 0.190442006391807, 0.190442006391807};
-	static double l14[]={
+	static IssmPDouble l14[]={
 		0.108103018168070, 0.445948490915965, 0.445948490915965, 
 		0.816847572980459, 0.091576213509771, 0.091576213509771};
-	static double l24[]={
+	static IssmPDouble l24[]={
 		0.445948490915965, 0.108103018168070, 0.445948490915965, 
 		0.091576213509771, 0.816847572980459, 0.091576213509771};
-	static double l34[]={
+	static IssmPDouble l34[]={
 		0.445948490915965, 0.445948490915965, 0.108103018168070, 
 		0.091576213509771, 0.091576213509771, 0.816847572980459};
 	/*}}}*/
 	/*p= 5, npoint= 7  {{{*/
-	static double wgt5[]={
+	static IssmPDouble wgt5[]={
 		0.389711431702997, 0.229313399254729, 0.229313399254729, 
 		0.229313399254729, 0.218133059367230, 0.218133059367230, 
 		0.218133059367230};
-	static double l15[]={
+	static IssmPDouble l15[]={
 		0.333333333333333, 0.059715871789770, 0.470142064105115, 
 		0.470142064105115, 0.797426985353087, 0.101286507323456, 
 		0.101286507323456};
-	static double l25[]={
+	static IssmPDouble l25[]={
 		0.333333333333333, 0.470142064105115, 0.059715871789770, 
 		0.470142064105115, 0.101286507323456, 0.797426985353087, 
 		0.101286507323456};
-	static double l35[]={
+	static IssmPDouble l35[]={
 		0.333333333333333, 0.470142064105115, 0.470142064105115, 
 		0.059715871789770, 0.101286507323456, 0.101286507323456, 
@@ -169,20 +169,20 @@
 	/*}}}*/
 	/*p= 6, npoint=12  {{{*/
-	static double wgt6[]={
+	static IssmPDouble wgt6[]={
 		0.202279763184836, 0.202279763184836, 0.202279763184836, 
 		0.088065961139281, 0.088065961139281, 0.088065961139281, 
 		0.143502272432755, 0.143502272432755, 0.143502272432755, 
 		0.143502272432755, 0.143502272432755, 0.143502272432755};
-	static double l16[]={
+	static IssmPDouble l16[]={
 		0.501426509658179, 0.249286745170910, 0.249286745170910, 
 		0.873821971016996, 0.063089014491502, 0.063089014491502, 
 		0.053145049844817, 0.053145049844817, 0.310352451033784, 
 		0.636502499121399, 0.310352451033784, 0.636502499121399};
-	static double l26[]={
+	static IssmPDouble l26[]={
 		0.249286745170910, 0.501426509658179, 0.249286745170910, 
 		0.063089014491502, 0.873821971016996, 0.063089014491502, 
 		0.310352451033784, 0.636502499121399, 0.053145049844817, 
 		0.053145049844817, 0.636502499121399, 0.310352451033784};
-	static double l36[]={
+	static IssmPDouble l36[]={
 		0.249286745170910, 0.249286745170910, 0.501426509658179, 
 		0.063089014491502, 0.063089014491502, 0.873821971016996, 
@@ -191,5 +191,5 @@
 	/*}}}*/
 	/*p= 7, npoint=13  {{{*/
-	static double wgt7[]={
+	static IssmPDouble wgt7[]={
 		-0.259062916308362, 0.304174548458604, 0.304174548458604, 
 		0.304174548458604, 0.092400122517855, 0.092400122517855, 
@@ -197,5 +197,5 @@
 		0.133564951824643, 0.133564951824643, 0.133564951824643, 
 		0.133564951824643};
-	static double l17[]={
+	static IssmPDouble l17[]={
 		0.333333333333333, 0.479308067841920, 0.260345966079040, 
 		0.260345966079040, 0.869739794195568, 0.065130102902216, 
@@ -203,5 +203,5 @@
 		0.312865496004874, 0.638444188569810, 0.312865496004874, 
 		0.638444188569810};
-	static double l27[]={
+	static IssmPDouble l27[]={
 		0.333333333333333, 0.260345966079040, 0.479308067841920, 
 		0.260345966079040, 0.065130102902216, 0.869739794195568, 
@@ -209,5 +209,5 @@
 		0.048690315425316, 0.048690315425316, 0.638444188569810, 
 		0.312865496004874};
-	static double l37[]={
+	static IssmPDouble l37[]={
 		0.333333333333333, 0.260345966079040, 0.260345966079040, 
 		0.479308067841920, 0.065130102902216, 0.065130102902216, 
@@ -217,5 +217,5 @@
 	/*}}}*/
 	/*p= 8, npoint=16  {{{*/
-	static double wgt8[]={
+	static IssmPDouble wgt8[]={
 		0.249961964823104, 0.164703541925695, 0.164703541925695, 
 		0.164703541925695, 0.178777729989794, 0.178777729989794, 
@@ -224,5 +224,5 @@
 		0.047164287656184, 0.047164287656184, 0.047164287656184, 
 		0.047164287656184};
-	static double l18[]={
+	static IssmPDouble l18[]={
 		0.333333333333333, 0.081414823414554, 0.459292588292723, 
 		0.459292588292723, 0.658861384496480, 0.170569307751760, 
@@ -231,5 +231,5 @@
 		0.263112829634638, 0.728492392955404, 0.263112829634638, 
 		0.728492392955404};
-	static double l28[]={
+	static IssmPDouble l28[]={
 		0.333333333333333, 0.459292588292723, 0.081414823414554, 
 		0.459292588292723, 0.170569307751760, 0.658861384496480, 
@@ -238,5 +238,5 @@
 		0.008394777409958, 0.008394777409958, 0.728492392955404, 
 		0.263112829634638};
-	static double l38[]={
+	static IssmPDouble l38[]={
 		0.333333333333333, 0.459292588292723, 0.459292588292723, 
 		0.081414823414554, 0.170569307751760, 0.170569307751760, 
@@ -247,5 +247,5 @@
 	/*}}}*/
 	/*p= 9, npoint=19  {{{*/
-	static double wgt9[]={
+	static IssmPDouble wgt9[]={
 		0.168244134395468, 0.054273292833345, 0.054273292833345, 
 		0.054273292833345, 0.134801255248419, 0.134801255248419, 
@@ -255,5 +255,5 @@
 		0.074969289332873, 0.074969289332873, 0.074969289332873, 
 		0.074969289332873};
-	static double l19[]={
+	static IssmPDouble l19[]={
 		0.333333333333333, 0.020634961602525, 0.489682519198738, 
 		0.489682519198738, 0.125820817014127, 0.437089591492937, 
@@ -263,5 +263,5 @@
 		0.221962989160766, 0.741198598784498, 0.221962989160766, 
 		0.741198598784498};
-	static double l29[]={
+	static IssmPDouble l29[]={
 		0.333333333333333, 0.489682519198738, 0.020634961602525, 
 		0.489682519198738, 0.437089591492937, 0.125820817014127, 
@@ -271,5 +271,5 @@
 		0.036838412054736, 0.036838412054736, 0.741198598784498, 
 		0.221962989160766};
-	static double l39[]={
+	static IssmPDouble l39[]={
 		0.333333333333333, 0.489682519198738, 0.489682519198738, 
 		0.020634961602525, 0.437089591492937, 0.437089591492937, 
@@ -281,5 +281,5 @@
 	/*}}}*/
 	/*p=10, npoint=25  {{{*/
-	static double wgt10[]={
+	static IssmPDouble wgt10[]={
 		0.157301373584232, 0.063611224790829, 0.063611224790829, 
 		0.063611224790829, 0.078498377595183, 0.078498377595183, 
@@ -291,5 +291,5 @@
 		0.016318805873179, 0.016318805873179, 0.016318805873179, 
 		0.016318805873179};
-	static double l110[]={
+	static IssmPDouble l110[]={
 		0.333333333333333, 0.028844733232685, 0.485577633383657, 
 		0.485577633383657, 0.781036849029926, 0.109481575485037, 
@@ -301,5 +301,5 @@
 		0.066803251012200, 0.923655933587500, 0.066803251012200, 
 		0.923655933587500};
-	static double l210[]={
+	static IssmPDouble l210[]={
 		0.333333333333333, 0.485577633383657, 0.028844733232685, 
 		0.485577633383657, 0.109481575485037, 0.781036849029926, 
@@ -311,5 +311,5 @@
 		0.009540815400299, 0.009540815400299, 0.923655933587500, 
 		0.066803251012200};
-	static double l310[]={
+	static IssmPDouble l310[]={
 		0.333333333333333, 0.485577633383657, 0.485577633383657, 
 		0.028844733232685, 0.109481575485037, 0.109481575485037, 
@@ -323,5 +323,5 @@
 	/*}}}*/
 	/*p=11, npoint=27  {{{*/
-	static double wgt11[]={
+	static IssmPDouble wgt11[]={
 		0.001605622060698, 0.001605622060698, 0.001605622060698, 
 		0.133626914252765, 0.133626914252765, 0.133626914252765, 
@@ -333,5 +333,5 @@
 		0.035866718600836, 0.035866718600836, 0.035866718600836, 
 		0.035866718600836, 0.035866718600836, 0.035866718600836};
-	static double l111[]={
+	static IssmPDouble l111[]={
 		-0.069222096541517, 0.534611048270758, 0.534611048270758, 
 		0.202061394068290, 0.398969302965855, 0.398969302965855, 
@@ -343,5 +343,5 @@
 		0.021022016536166, 0.021022016536166, 0.171488980304042, 
 		0.807489003159792, 0.171488980304042, 0.807489003159792};
-	static double l211[]={
+	static IssmPDouble l211[]={
 		0.534611048270758,-0.069222096541517, 0.534611048270758, 
 		0.398969302965855, 0.202061394068290, 0.398969302965855, 
@@ -353,5 +353,5 @@
 		0.171488980304042, 0.807489003159792, 0.021022016536166, 
 		0.021022016536166, 0.807489003159792, 0.171488980304042};
-	static double l311[]={
+	static IssmPDouble l311[]={
 		0.534611048270758, 0.534611048270758,-0.069222096541517, 
 		0.398969302965855, 0.398969302965855, 0.202061394068290, 
@@ -365,5 +365,5 @@
 	/*}}}*/
 	/*p=12, npoint=33  {{{*/
-	static double wgt12[]={
+	static IssmPDouble wgt12[]={
 		0.044567514407799, 0.044567514407799, 0.044567514407799, 
 		0.075677707051848, 0.075677707051848, 0.075677707051848, 
@@ -377,5 +377,5 @@
 		0.029992592075802, 0.029992592075802, 0.029992592075802, 
 		0.029992592075802, 0.029992592075802, 0.029992592075802};
-	static double l112[]={
+	static IssmPDouble l112[]={
 		0.023565220452390, 0.488217389773805, 0.488217389773805, 
 		0.120551215411079, 0.439724392294460, 0.439724392294460, 
@@ -389,5 +389,5 @@
 		0.025734050548330, 0.025734050548330, 0.116251915907597, 
 		0.858014033544073, 0.116251915907597, 0.858014033544073};
-	static double l212[]={
+	static IssmPDouble l212[]={
 		0.488217389773805, 0.023565220452390, 0.488217389773805, 
 		0.439724392294460, 0.120551215411079, 0.439724392294460, 
@@ -401,5 +401,5 @@
 		0.116251915907597, 0.858014033544073, 0.025734050548330, 
 		0.025734050548330, 0.858014033544073, 0.116251915907597};
-	static double l312[]={
+	static IssmPDouble l312[]={
 		0.488217389773805, 0.488217389773805, 0.023565220452390, 
 		0.439724392294460, 0.439724392294460, 0.120551215411079, 
@@ -415,5 +415,5 @@
 	/*}}}*/
 	/*  p=13, npoint=37  {{{*/
-	static double wgt13[]={
+	static IssmPDouble wgt13[]={
 		0.090968907790622, 0.019537784619314, 0.019537784619314, 
 		0.019537784619314, 0.054427130356344, 0.054427130356344, 
@@ -429,5 +429,5 @@
 		0.026884523429480, 0.026884523429480, 0.026884523429480, 
 		0.026884523429480};
-	static double l113[]={
+	static IssmPDouble l113[]={
 		0.333333333333333, 0.009903630120591, 0.495048184939705, 
 		0.495048184939705, 0.062566729780852, 0.468716635109574, 
@@ -443,5 +443,5 @@
 		0.126357385491669, 0.851409537834241, 0.126357385491669, 
 		0.851409537834241};
-	static double l213[]={
+	static IssmPDouble l213[]={
 		0.333333333333333, 0.495048184939705, 0.009903630120591, 
 		0.495048184939705, 0.468716635109574, 0.062566729780852, 
@@ -457,5 +457,5 @@
 		0.022233076674090, 0.022233076674090, 0.851409537834241, 
 		0.126357385491669};
-	static double l313[]={
+	static IssmPDouble l313[]={
 		0.333333333333333, 0.495048184939705, 0.495048184939705, 
 		0.009903630120591, 0.468716635109574, 0.468716635109574, 
@@ -473,5 +473,5 @@
 	/*}}}*/
 	/*p=14, npoint=42{{{*/
-	static double wgt14[]={
+	static IssmPDouble wgt14[]={
 		0.037903474783419, 0.037903474783419, 0.037903474783419, 
 		0.056791094234956, 0.056791094234956, 0.056791094234956, 
@@ -488,5 +488,5 @@
 		0.008677970905831, 0.008677970905831, 0.008677970905831, 
 		0.008677970905831, 0.008677970905831, 0.008677970905831};
-	static double l114[]={
+	static IssmPDouble l114[]={
 		0.022072179275643, 0.488963910362179, 0.488963910362179, 
 		0.164710561319092, 0.417644719340454, 0.417644719340454, 
@@ -503,5 +503,5 @@
 		0.001268330932872, 0.001268330932872, 0.118974497696957, 
 		0.879757171370171, 0.118974497696957, 0.879757171370171};
-	static double l214[]={
+	static IssmPDouble l214[]={
 		0.488963910362179, 0.022072179275643, 0.488963910362179, 
 		0.417644719340454, 0.164710561319092, 0.417644719340454, 
@@ -518,5 +518,5 @@
 		0.118974497696957, 0.879757171370171, 0.001268330932872, 
 		0.001268330932872, 0.879757171370171, 0.118974497696957};
-	static double l314[]={
+	static IssmPDouble l314[]={
 		0.488963910362179, 0.488963910362179, 0.022072179275643, 
 		0.417644719340454, 0.417644719340454, 0.164710561319092, 
@@ -535,5 +535,5 @@
 	/*}}}*/
 	/*p=15, npoint=48{{{*/
-	static double wgt15[]={
+	static IssmPDouble wgt15[]={
 		0.003320126005206, 0.003320126005206, 0.003320126005206, 
 		0.076641563419124, 0.076641563419124, 0.076641563419124, 
@@ -552,5 +552,5 @@
 		0.013291658531346, 0.013291658531346, 0.013291658531346, 
 		0.013291658531346, 0.013291658531346, 0.013291658531346};
-	static double l115[]={
+	static IssmPDouble l115[]={
 		-0.013945833716486, 0.506972916858243, 0.506972916858243, 
 		0.137187291433955, 0.431406354283023, 0.431406354283023, 
@@ -569,5 +569,5 @@
 		0.012459809331199, 0.012459809331199, 0.103575616576386, 
 		0.883964574092416, 0.103575616576386, 0.883964574092416};
-	static double l215[]={
+	static IssmPDouble l215[]={
 		0.506972916858243,-0.013945833716486, 0.506972916858243, 
 		0.431406354283023, 0.137187291433955, 0.431406354283023, 
@@ -586,5 +586,5 @@
 		0.103575616576386, 0.883964574092416, 0.012459809331199, 
 		0.012459809331199, 0.883964574092416, 0.103575616576386};
-	static double l315[]={
+	static IssmPDouble l315[]={
 		0.506972916858243, 0.506972916858243,-0.013945833716486, 
 		0.431406354283023, 0.431406354283023, 0.137187291433955, 
@@ -605,5 +605,5 @@
 	/*}}}*/
 	/*p=16, npoint=52  {{{*/
-	static double wgt16[]={
+	static IssmPDouble wgt16[]={
 		0.081191089584902, 0.011095307165226, 0.011095307165226, 
 		0.011095307165226, 0.072244353151393, 0.072244353151393, 
@@ -624,5 +624,5 @@
 		0.011864642509229, 0.011864642509229, 0.011864642509229, 
 		0.011864642509229};
-	static double l116[]={
+	static IssmPDouble l116[]={
 		0.333333333333333, 0.005238916103123, 0.497380541948438, 
 		0.497380541948438, 0.173061122901295, 0.413469438549352, 
@@ -643,5 +643,5 @@
 		0.085283615682657, 0.900399064086661, 0.085283615682657, 
 		0.900399064086661};
-	static double l216[]={
+	static IssmPDouble l216[]={
 		0.333333333333333, 0.497380541948438, 0.005238916103123, 
 		0.497380541948438, 0.413469438549352, 0.173061122901295, 
@@ -662,5 +662,5 @@
 		0.014317320230681, 0.014317320230681, 0.900399064086661, 
 		0.085283615682657};
-	static double l316[]={
+	static IssmPDouble l316[]={
 		0.333333333333333, 0.497380541948438, 0.497380541948438, 
 		0.005238916103123, 0.413469438549352, 0.413469438549352, 
@@ -683,5 +683,5 @@
 	/*}}}*/
 	/*p=17, npoint=61{{{*/
-	static double wgt17[]={
+	static IssmPDouble wgt17[]={
 		0.057914928034477, 0.008822054327014, 0.008822054327014, 
 		0.008822054327014, 0.025410682752829, 0.025410682752829, 
@@ -705,5 +705,5 @@
 		0.011545213295771, 0.011545213295771, 0.011545213295771, 
 		0.011545213295771};
-	static double l117[]={
+	static IssmPDouble l117[]={
 		0.333333333333333, 0.005658918886452, 0.497170540556774, 
 		0.497170540556774, 0.035647354750751, 0.482176322624625, 
@@ -727,5 +727,5 @@
 		0.080711313679564, 0.904625504095608, 0.080711313679564, 
 		0.904625504095608};
-	static double l217[]={
+	static IssmPDouble l217[]={
 		0.333333333333333, 0.497170540556774, 0.005658918886452, 
 		0.497170540556774, 0.482176322624625, 0.035647354750751, 
@@ -749,5 +749,5 @@
 		0.014663182224828, 0.014663182224828, 0.904625504095608, 
 		0.080711313679564};
-	static double l317[]={
+	static IssmPDouble l317[]={
 		0.333333333333333, 0.497170540556774, 0.497170540556774, 
 		0.005658918886452, 0.482176322624625, 0.482176322624625, 
@@ -774,5 +774,5 @@
 	/*  p=18, npoint=70  {{{*/
 
-	static double wgt18[]={
+	static IssmPDouble wgt18[]={
 		0.053364381350150, 0.015713921277179, 0.015713921277179, 
 		0.015713921277179, 0.032495554156279, 0.032495554156279, 
@@ -799,5 +799,5 @@
 		0.000079999375178, 0.000079999375178, 0.000079999375178, 
 		0.000079999375178};
-	static double l118[]={
+	static IssmPDouble l118[]={
 		0.333333333333333, 0.013310382738157, 0.493344808630921, 
 		0.493344808630921, 0.061578811516086, 0.469210594241957, 
@@ -824,5 +824,5 @@
 		0.020874755282586, 1.014347260005363, 0.020874755282586, 
 		1.014347260005363};
-	static double l218[]={
+	static IssmPDouble l218[]={
 		0.333333333333333, 0.493344808630921, 0.013310382738157, 
 		0.493344808630921, 0.469210594241957, 0.061578811516086, 
@@ -849,5 +849,5 @@
 		-0.035222015287949,-0.035222015287949, 1.014347260005363, 
 		0.020874755282586};
-	static double l318[]={
+	static IssmPDouble l318[]={
 		0.333333333333333, 0.493344808630921, 0.493344808630921, 
 		0.013310382738157, 0.469210594241957, 0.469210594241957, 
@@ -877,5 +877,5 @@
 	/*p=19, npoint=73  {{{*/
 
-	static double wgt19[]={
+	static IssmPDouble wgt19[]={
 		0.056995437856306, 0.017893352515055, 0.017893352515055, 
 		0.017893352515055, 0.038775849701151, 0.038775849701151, 
@@ -903,5 +903,5 @@
 		0.006581669842530, 0.006581669842530, 0.006581669842530, 
 		0.006581669842530};
-	static double l119[]={
+	static IssmPDouble l119[]={
 		0.333333333333333, 0.020780025853987, 0.489609987073006, 
 		0.489609987073006, 0.090926214604215, 0.454536892697893, 
@@ -929,5 +929,5 @@
 		0.065494628082938, 0.924344252620784, 0.065494628082938, 
 		0.924344252620784};
-	static double l219[]={
+	static IssmPDouble l219[]={
 		0.333333333333333, 0.489609987073006, 0.020780025853987, 
 		0.489609987073006, 0.454536892697893, 0.090926214604215, 
@@ -955,5 +955,5 @@
 		0.010161119296278, 0.010161119296278, 0.924344252620784, 
 		0.065494628082938};
-	static double l319[]={
+	static IssmPDouble l319[]={
 		0.333333333333333, 0.489609987073006, 0.489609987073006, 
 		0.020780025853987, 0.454536892697893, 0.454536892697893, 
@@ -983,5 +983,5 @@
 	/*}}}*/
 	/*p=20, npoint=79 {{{*/
-	static double wgt20[]={
+	static IssmPDouble wgt20[]={
 		0.057256499746719, 0.001501721280705, 0.001501721280705, 
 		0.001501721280705, 0.020195803723819, 0.020195803723819, 
@@ -1011,5 +1011,5 @@
 		0.006190192638113, 0.006190192638113, 0.006190192638113, 
 		0.006190192638113};
-	static double l120[]={
+	static IssmPDouble l120[]={
 		0.333333333333333,-0.001900928704400, 0.500950464352200, 
 		0.500950464352200, 0.023574084130543, 0.488212957934729, 
@@ -1039,5 +1039,5 @@
 		0.059696109149007, 0.929756171556853, 0.059696109149007, 
 		0.929756171556853};
-	static double l220[]={
+	static IssmPDouble l220[]={
 		0.333333333333333, 0.500950464352200,-0.001900928704400, 
 		0.500950464352200, 0.488212957934729, 0.023574084130543, 
@@ -1067,5 +1067,5 @@
 		0.010547719294141, 0.010547719294141, 0.929756171556853, 
 		0.059696109149007};
-	static double l320[]={
+	static IssmPDouble l320[]={
 		0.333333333333333, 0.500950464352200, 0.500950464352200, 
 		-0.001900928704400, 0.488212957934729, 0.488212957934729, 
@@ -1097,20 +1097,20 @@
 	/*}}}*/
 
-	static double* wgtp[MAX_TRIA_SYM_ORD]={
+	static IssmPDouble* wgtp[MAX_TRIA_SYM_ORD]={
 		wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 ,
 		wgt6 ,wgt7 ,wgt8 ,wgt9 ,wgt10,
 		wgt11,wgt12,wgt13,wgt14,wgt15,
 		wgt16,wgt17,wgt18,wgt19,wgt20};
-	static double* l1p [MAX_TRIA_SYM_ORD]={
+	static IssmPDouble* l1p [MAX_TRIA_SYM_ORD]={
 		l11  ,l12  ,l13  ,l14  ,l15  ,
 		l16  ,l17  ,l18  ,l19  ,l110 ,
 		l111 ,l112 ,l113 ,l114 ,l115 ,
 		l116 ,l117 ,l118 ,l119 ,l120 };
-	static double* l2p [MAX_TRIA_SYM_ORD]={
+	static IssmPDouble* l2p [MAX_TRIA_SYM_ORD]={
 		l21  ,l22  ,l23  ,l24  ,l25  ,
 		l26  ,l27  ,l28  ,l29  ,l210 ,
 		l211 ,l212 ,l213 ,l214 ,l215 ,
 		l216 ,l217 ,l218 ,l219 ,l220 };
-	static double* l3p [MAX_TRIA_SYM_ORD]={
+	static IssmPDouble* l3p [MAX_TRIA_SYM_ORD]={
 		l31  ,l32  ,l33  ,l34  ,l35  ,
 		l36  ,l37  ,l38  ,l39  ,l310 ,
@@ -1118,24 +1118,24 @@
 		l316 ,l317 ,l318 ,l319 ,l320 };
 
-	static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(double),
-		sizeof(wgt2 )/sizeof(double),
-		sizeof(wgt3 )/sizeof(double),
-		sizeof(wgt4 )/sizeof(double),
-		sizeof(wgt5 )/sizeof(double),
-		sizeof(wgt6 )/sizeof(double),
-		sizeof(wgt7 )/sizeof(double),
-		sizeof(wgt8 )/sizeof(double),
-		sizeof(wgt9 )/sizeof(double),
-		sizeof(wgt10)/sizeof(double),
-		sizeof(wgt11)/sizeof(double),
-		sizeof(wgt12)/sizeof(double),
-		sizeof(wgt13)/sizeof(double),
-		sizeof(wgt14)/sizeof(double),
-		sizeof(wgt15)/sizeof(double),
-		sizeof(wgt16)/sizeof(double),
-		sizeof(wgt17)/sizeof(double),
-		sizeof(wgt18)/sizeof(double),
-		sizeof(wgt19)/sizeof(double),
-		sizeof(wgt20)/sizeof(double)};
+	static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
+		sizeof(wgt2 )/sizeof(IssmPDouble),
+		sizeof(wgt3 )/sizeof(IssmPDouble),
+		sizeof(wgt4 )/sizeof(IssmPDouble),
+		sizeof(wgt5 )/sizeof(IssmPDouble),
+		sizeof(wgt6 )/sizeof(IssmPDouble),
+		sizeof(wgt7 )/sizeof(IssmPDouble),
+		sizeof(wgt8 )/sizeof(IssmPDouble),
+		sizeof(wgt9 )/sizeof(IssmPDouble),
+		sizeof(wgt10)/sizeof(IssmPDouble),
+		sizeof(wgt11)/sizeof(IssmPDouble),
+		sizeof(wgt12)/sizeof(IssmPDouble),
+		sizeof(wgt13)/sizeof(IssmPDouble),
+		sizeof(wgt14)/sizeof(IssmPDouble),
+		sizeof(wgt15)/sizeof(IssmPDouble),
+		sizeof(wgt16)/sizeof(IssmPDouble),
+		sizeof(wgt17)/sizeof(IssmPDouble),
+		sizeof(wgt18)/sizeof(IssmPDouble),
+		sizeof(wgt19)/sizeof(IssmPDouble),
+		sizeof(wgt20)/sizeof(IssmPDouble)};
 
 	//	_printf_(true,"GaussLegendreTria: iord=%d\n",iord);
@@ -1150,8 +1150,8 @@
 		*pngaus=np[iord-1];
 
-		*pl1  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl2  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl3  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pwgt = (double *) xmalloc(*pngaus*sizeof(double));
+		*pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
 
 		for (i=0; i<*pngaus; i++) {
@@ -1168,8 +1168,8 @@
 		*pngaus=nigaus*nigaus;
 
-		*pl1  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl2  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl3  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pwgt = (double *) xmalloc(*pngaus*sizeof(double));
+		*pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
 
 		/*  get the gauss points in each direction  */
@@ -1207,5 +1207,5 @@
 }/*}}}*/
 /*FUNCTION GaussLegendreTetra{{{*/
-void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
+void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord ) {
 	/* Gauss quadrature points for the tetrahedron.
 
@@ -1220,52 +1220,52 @@
 	/*Intermediaries*/
 	int i,j,k,ipt,nigaus;
-	double xi,eta,zeta;
-	double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
+	IssmPDouble xi,eta,zeta;
+	IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
 
 	/*Hardcoded Gauss points definition*/
 	/*p= 1, npoint= 1  {{{*/
-	static double wgt1[]={
+	static IssmPDouble wgt1[]={
 		1.000000000000000};
-	static double l11[]={
+	static IssmPDouble l11[]={
 		0.250000000000000};
-	static double l21[]={
+	static IssmPDouble l21[]={
 		0.250000000000000};
-	static double l31[]={
+	static IssmPDouble l31[]={
 		0.250000000000000};
-	static double l41[]={
+	static IssmPDouble l41[]={
 		0.250000000000000};
 	/*}}}*/
 	/*p= 2, npoint= 4  {{{*/
 
-	static double wgt2[]={
+	static IssmPDouble wgt2[]={
 		0.250000000000000, 0.250000000000000, 0.250000000000000,
 		0.250000000000000};
-	static double l12[]={
+	static IssmPDouble l12[]={
 		0.585410196624969, 0.138196601125011, 0.138196601125011,
 		0.138196601125011};
-	static double l22[]={
+	static IssmPDouble l22[]={
 		0.138196601125011, 0.585410196624969, 0.138196601125011,
 		0.138196601125011};
-	static double l32[]={
+	static IssmPDouble l32[]={
 		0.138196601125011, 0.138196601125011, 0.585410196624969,
 		0.138196601125011};
-	static double l42[]={
+	static IssmPDouble l42[]={
 		0.138196601125011, 0.138196601125011, 0.138196601125011,
 		0.585410196624969};
 	/*}}}*/
 	/*p= 3, npoint= 5  {{{*/
-	static double wgt3[]={
+	static IssmPDouble wgt3[]={
 		-0.800000000000000, 0.450000000000000, 0.450000000000000, 
 		0.450000000000000, 0.450000000000000};
-	static double l13[]={
+	static IssmPDouble l13[]={
 		0.250000000000000, 0.500000000000000, 0.166666666666667, 
 		0.166666666666667, 0.166666666666667};
-	static double l23[]={
+	static IssmPDouble l23[]={
 		0.250000000000000, 0.166666666666667, 0.500000000000000, 
 		0.166666666666667, 0.166666666666667};
-	static double l33[]={
+	static IssmPDouble l33[]={
 		0.250000000000000, 0.166666666666667, 0.166666666666667, 
 		0.500000000000000, 0.166666666666667};
-	static double l43[]={
+	static IssmPDouble l43[]={
 		0.250000000000000, 0.166666666666667, 0.166666666666667, 
 		0.166666666666667, 0.500000000000000};
@@ -1273,25 +1273,25 @@
 	/*p= 4, npoint=11  {{{*/
 
-	static double wgt4[]={
+	static IssmPDouble wgt4[]={
 		-0.013155555555556, 0.007622222222222, 0.007622222222222, 
 		0.007622222222222, 0.007622222222222, 0.024888888888889,
 		0.024888888888889, 0.024888888888889, 0.024888888888889, 
 		0.024888888888889, 0.024888888888889};
-	static double l14[]={
+	static IssmPDouble l14[]={
 		0.250000000000000, 0.785714285714286, 0.071428571428571, 
 		0.071428571428571, 0.071428571428571, 0.399403576166799,
 		0.399403576166799, 0.399403576166799, 0.100596423833201, 
 		0.100596423833201, 0.100596423833201};
-	static double l24[]={
+	static IssmPDouble l24[]={
 		0.250000000000000, 0.071428571428571, 0.785714285714286, 
 		0.071428571428571, 0.071428571428571, 0.399403576166799,
 		0.100596423833201, 0.100596423833201, 0.399403576166799, 
 		0.399403576166799, 0.100596423833201};
-	static double l34[]={
+	static IssmPDouble l34[]={
 		0.250000000000000, 0.071428571428571, 0.071428571428571, 
 		0.785714285714286, 0.071428571428571, 0.100596423833201, 
 		0.399403576166799, 0.100596423833201, 0.399403576166799,
 		0.100596423833201, 0.399403576166799};
-	static double l44[]={
+	static IssmPDouble l44[]={
 		0.250000000000000, 0.071428571428571, 0.071428571428571, 
 		0.071428571428571, 0.785714285714286, 0.100596423833201, 
@@ -1301,5 +1301,5 @@
 	/*p= 5, npoint=15  {{{*/
 
-	static double wgt5[]={
+	static IssmPDouble wgt5[]={
 		0.030283678097089, 0.006026785714286, 0.006026785714286, 
 		0.006026785714286, 0.006026785714286, 0.011645249086029, 
@@ -1307,5 +1307,5 @@
 		0.010949141561386, 0.010949141561386, 0.010949141561386,
 		0.010949141561386, 0.010949141561386, 0.010949141561386};
-	static double l15[]={
+	static IssmPDouble l15[]={
 		0.250000000000000, 0.000000000000000, 0.333333333333333, 
 		0.333333333333333, 0.333333333333333, 0.727272727272727, 
@@ -1313,5 +1313,5 @@
 		0.066550153573664, 0.066550153573664, 0.066550153573664, 
 		0.433449846426336, 0.433449846426336, 0.433449846426336};
-	static double l25[]={
+	static IssmPDouble l25[]={
 		0.250000000000000, 0.333333333333333, 0.000000000000000, 
 		0.333333333333333, 0.333333333333333, 0.090909090909091, 
@@ -1319,5 +1319,5 @@
 		0.066550153573664, 0.433449846426336, 0.433449846426336, 
 		0.066550153573664, 0.066550153573664, 0.433449846426336};
-	static double l35[]={
+	static IssmPDouble l35[]={
 		0.250000000000000, 0.333333333333333, 0.333333333333333, 
 		0.000000000000000, 0.333333333333333, 0.090909090909091, 
@@ -1325,5 +1325,5 @@
 		0.433449846426336, 0.066550153573664, 0.433449846426336, 
 		0.066550153573664, 0.433449846426336, 0.066550153573664};
-	static double l45[]={
+	static IssmPDouble l45[]={
 		0.250000000000000, 0.333333333333333, 0.333333333333333, 
 		0.333333333333333, 0.000000000000000, 0.090909090909091, 
@@ -1334,5 +1334,5 @@
 	/*p= 6, npoint=24  {{{*/
 
-	static double wgt6[]={
+	static IssmPDouble wgt6[]={
 		0.006653791709695, 0.006653791709695, 0.006653791709695, 
 		0.006653791709695, 0.001679535175887, 0.001679535175887, 
@@ -1343,5 +1343,5 @@
 		0.008035714285714, 0.008035714285714, 0.008035714285714, 
 		0.008035714285714, 0.008035714285714, 0.008035714285714};
-	static double l16[]={
+	static IssmPDouble l16[]={
 		0.356191386222545, 0.214602871259152, 0.214602871259152, 
 		0.214602871259152, 0.877978124396166, 0.040673958534611, 
@@ -1353,5 +1353,5 @@
 		0.269672331458316, 0.603005664791649, 0.269672331458316, 
 		0.603005664791649, 0.269672331458316, 0.603005664791649};
-	static double l26[]={
+	static IssmPDouble l26[]={
 		0.214602871259152, 0.356191386222545, 0.214602871259152, 
 		0.214602871259152, 0.040673958534611, 0.877978124396166, 
@@ -1363,5 +1363,5 @@
 		0.063661001875018, 0.063661001875018, 0.063661001875018, 
 		0.063661001875018, 0.603005664791649, 0.269672331458316};
-	static double l36[]={
+	static IssmPDouble l36[]={
 		0.214602871259152, 0.214602871259152, 0.356191386222545, 
 		0.214602871259152, 0.040673958534611, 0.040673958534611, 
@@ -1373,5 +1373,5 @@
 		0.063661001875018, 0.063661001875018, 0.603005664791649, 
 		0.269672331458316, 0.063661001875018, 0.063661001875018};
-	static double l46[]={
+	static IssmPDouble l46[]={
 		0.214602871259152, 0.214602871259152, 0.214602871259152, 
 		0.356191386222545, 0.040673958534611, 0.040673958534611, 
@@ -1385,16 +1385,16 @@
 	/*}}}*/
 
-	static double* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};
-	static double* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };
-	static double* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 };
-	static double* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };
-	static double* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };
-
-	static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(double),
-		sizeof(wgt2 )/sizeof(double),
-		sizeof(wgt3 )/sizeof(double),
-		sizeof(wgt4 )/sizeof(double),
-		sizeof(wgt5 )/sizeof(double),
-		sizeof(wgt6 )/sizeof(double)};
+	static IssmPDouble* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};
+	static IssmPDouble* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };
+	static IssmPDouble* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 };
+	static IssmPDouble* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };
+	static IssmPDouble* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };
+
+	static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
+		sizeof(wgt2 )/sizeof(IssmPDouble),
+		sizeof(wgt3 )/sizeof(IssmPDouble),
+		sizeof(wgt4 )/sizeof(IssmPDouble),
+		sizeof(wgt5 )/sizeof(IssmPDouble),
+		sizeof(wgt6 )/sizeof(IssmPDouble)};
 
 	//	_printf_(true,"GaussLegendreTetra: iord=%d\n",iord);
@@ -1410,9 +1410,9 @@
 		*pngaus=np[iord-1];
 
-		*pl1  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl2  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl3  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl4  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pwgt = (double *) xmalloc(*pngaus*sizeof(double));
+		*pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl4  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
 
 		for (i=0; i<*pngaus; i++) {
@@ -1430,9 +1430,9 @@
 		*pngaus=nigaus*nigaus*nigaus;
 
-		*pl1  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl2  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl3  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pl4  = (double *) xmalloc(*pngaus*sizeof(double));
-		*pwgt = (double *) xmalloc(*pngaus*sizeof(double));
+		*pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pl4  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
+		*pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
 
 		/*  get the gauss points in each direction  */
@@ -1472,5 +1472,5 @@
 }/*}}}*/
 /*FUNCTION GaussLobatto{{{*/
-void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
+void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus ) {
 	/*Gauss-Lobatto quadrature points.
 
@@ -1488,40 +1488,40 @@
 
 	int i;
-	double *alpha,*beta;
-	double left=-1.,right= 1.;
-	double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
+	IssmPDouble *alpha,*beta;
+	IssmPDouble left=-1.,right= 1.;
+	IssmPDouble p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
 
 	/*p= 1, npoint= 1 (Gauss-Legendre)*/
-	static double wgt1[]={2.000000000000000};
-	static double xi1[]={0.000000000000000};
+	static IssmPDouble wgt1[]={2.000000000000000};
+	static IssmPDouble xi1[]={0.000000000000000};
 
 	/*p= 1, npoint= 2*/
-	static double wgt2[]={1.000000000000000, 1.000000000000000};
-	static double xi2[]={-1.000000000000000, 1.000000000000000};
+	static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
+	static IssmPDouble xi2[]={-1.000000000000000, 1.000000000000000};
 
 	/*p= 3, npoint= 3*/
-	static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
-	static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
+	static IssmPDouble wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
+	static IssmPDouble xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
 
 	/*p= 5, npoint= 4*/
-	static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
-	static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
+	static IssmPDouble wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
+	static IssmPDouble xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
 
 	/*p= 7, npoint= 5*/
-	static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
-	static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
-
-	static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
-	static double* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
-
-	static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(double),
-		sizeof(wgt2 )/sizeof(double),
-		sizeof(wgt3 )/sizeof(double),
-		sizeof(wgt4 )/sizeof(double),
-		sizeof(wgt5 )/sizeof(double)};
+	static IssmPDouble wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
+	static IssmPDouble xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
+
+	static IssmPDouble* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
+	static IssmPDouble* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
+
+	static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble),
+		sizeof(wgt2 )/sizeof(IssmPDouble),
+		sizeof(wgt3 )/sizeof(IssmPDouble),
+		sizeof(wgt4 )/sizeof(IssmPDouble),
+		sizeof(wgt5 )/sizeof(IssmPDouble)};
 
 	//	_printf_(true,"Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
-	*pxgaus = (double *) xmalloc(ngaus*sizeof(double));
-	*pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
+	*pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
+	*pxwgt  = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
 
 	/*  check to see if Gauss points need to be calculated  */
@@ -1539,6 +1539,6 @@
 
 		/*  calculate the Gauss points using recurrence relations  */
-		alpha=(double *) xmalloc(ngaus*sizeof(double));
-		beta =(double *) xmalloc(ngaus*sizeof(double));
+		alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
+		beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
 
 		/*  calculate the Legendre recurrence coefficients  */
@@ -1581,5 +1581,5 @@
 }/*}}}*/
 /*FUNCTION GaussRecur{{{*/
-void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
+void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta ) {
 	/*Gauss quadrature points from recursion coefficients.
 	 *
@@ -1589,6 +1589,6 @@
 	/*Intermediaries*/
 	int i,j,k,l,m,ii,mml,iter;
-	double p,g,r,s,c,f,b;
-	double* work;
+	IssmPDouble p,g,r,s,c,f,b;
+	IssmPDouble* work;
 
 	if (n==1){
@@ -1598,5 +1598,5 @@
 	}
 
-	work=(double*)xmalloc(n*sizeof(double));
+	work=(IssmPDouble*)xmalloc(n*sizeof(IssmPDouble));
 
 	zero[0]  =alpha[0];
@@ -1704,5 +1704,5 @@
 /*Element Gauss points TO BE REMOVED*/
 /*FUNCTION gaussQuad{{{*/
-void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) { 
+void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ) {
 	/*Gauss quadrature points for the quadrilaterial.*/
 
@@ -1712,5 +1712,5 @@
 }/*}}}*/
 /*FUNCTION gaussHexa{{{*/
-void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
+void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ) {
 	/*Gauss quadrature points for the hexahedron.*/
 
Index: /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.h	(revision 12423)
+++ /issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.h	(revision 12424)
@@ -7,16 +7,16 @@
 
 #define MAX_LINE_GAUS_PTS    4
-void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus ); 
+void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus );
 #define MAX_TRIA_SYM_ORD    20
-void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
+void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord );
 #define MAX_TETRA_SYM_ORD    6
-void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
+void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord );
 #define MAX_LINE_GLOB_PTS    5
-void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ); 
+void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus );
 #define MAX_GAUS_ITER   30
-void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );
+void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta );
 
-void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
-void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
+void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus );
+void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus );
 
 #endif