Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 13080)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 13081)
@@ -370,5 +370,5 @@
 		if(stabilization==2){
 			/*Streamline upwinding*/
-			vel=sqrt(pow(vx,2.)+pow(vy,2.))+1.e-8;
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
 			K[0][0]=h/(2*vel)*vx*vx;
 			K[1][0]=h/(2*vel)*vy*vx;
@@ -2328,5 +2328,4 @@
 IssmDouble Tria::SurfaceArea(void){
 
-	int    i;
 	IssmDouble S;
 	IssmDouble normal[3];
@@ -2339,5 +2338,5 @@
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 
-	for (i=0;i<3;i++){
+	for(int i=0;i<3;i++){
 		v13[i]=xyz_list[0][i]-xyz_list[2][i];
 		v23[i]=xyz_list[1][i]-xyz_list[2][i];
@@ -2348,5 +2347,5 @@
 	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
 
-	S = 0.5 * sqrt(pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2));
+	S = 0.5 * sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
 
 	/*Return: */
@@ -2357,10 +2356,9 @@
 void Tria::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){
 
-	int i;
 	IssmDouble v13[3],v23[3];
 	IssmDouble normal[3];
 	IssmDouble normal_norm;
 
-	for (i=0;i<3;i++){
+	for(int i=0;i<3;i++){
 		v13[i]=xyz_list[0][i]-xyz_list[2][i];
 		v23[i]=xyz_list[1][i]-xyz_list[2][i];
@@ -2371,9 +2369,9 @@
 	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
 
-	normal_norm=sqrt( pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2) );
-
-	*(surface_normal)=normal[0]/normal_norm;
-	*(surface_normal+1)=normal[1]/normal_norm;
-	*(surface_normal+2)=normal[2]/normal_norm;
+	normal_norm=sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
+
+	*(surface_normal+0) = normal[0]/normal_norm;
+	*(surface_normal+1) = normal[1]/normal_norm;
+	*(surface_normal+2) = normal[2]/normal_norm;
 }
 /*}}}*/
@@ -2608,5 +2606,5 @@
 	normal[1]=sin(atan2(x1-x2,y2-y1));
 
-	length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
+	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
 
 	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
@@ -2974,5 +2972,5 @@
 		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
 		surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+		slope_magnitude=sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
 		if(slope_magnitude>MAXSLOPE) alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
 		else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
@@ -3304,5 +3302,5 @@
 	/*Get Vz and compute vel*/
 	GetInputListOnVertices(&vz[0],VzEnum,0);
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
@@ -3364,5 +3362,5 @@
 	/*Now Compute vel*/
 	GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
@@ -3983,5 +3981,5 @@
 		 *      S                obs            obs
 		 */
-		misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
+		misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2));
 
 		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
@@ -4046,7 +4044,7 @@
 		 *                            obs
 		 */
-		velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
-		obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
-		misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
+		velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
+		obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
+		misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
 
 		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum);
@@ -4112,7 +4110,7 @@
 		 *                              obs                       obs
 		 */
-		misfit=0.5*pow(meanvel,2.)*(
-					pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
-					pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
+		misfit=0.5*pow(meanvel,2)*(
+					pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2) +
+					pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2) );
 
 		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum);
@@ -4175,5 +4173,5 @@
 		 *
 		 */
-		misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
+		misfit=0.5*( pow(vx-vxobs,2) + pow(vy-vyobs,2) );
 
 		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
@@ -4238,7 +4236,7 @@
 		 *              obs                        obs                      
 		 */
-		scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
-		scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
-		misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
+		scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+		scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
+		misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
 		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum);
 
@@ -4288,5 +4286,5 @@
 
 		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
-		Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+		Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
 	}
 
@@ -4437,5 +4435,5 @@
 
 		/*compute ThicknessAbsMisfit*/
-		Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
+		Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight;
 	}
 
@@ -4625,6 +4623,6 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
-						scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
+						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
 						dux=scalex*(vxobs-vx);
 						duy=scaley*(vyobs-vy);
@@ -4646,7 +4644,7 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
-						obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
-						scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
+						velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
+						obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
+						scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
 						dux=scale*vx;
 						duy=scale*vy;
@@ -4666,5 +4664,5 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
+						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
 						dux=scale*(vxobs-vx);
 						duy=scale*(vyobs-vy);
@@ -4684,6 +4682,6 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
 						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
 						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
@@ -4808,6 +4806,6 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
-						scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
+						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
 						dux=scalex*(vxobs-vx);
 						duy=scaley*(vyobs-vy);
@@ -4829,7 +4827,7 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
-						obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
-						scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
+						velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
+						obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
+						scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
 						dux=scale*vx;
 						duy=scale*vy;
@@ -4849,5 +4847,5 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
+						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
 						dux=scale*(vxobs-vx);
 						duy=scale*(vyobs-vy);
@@ -4867,6 +4865,6 @@
 					 */
 					for (i=0;i<NUMVERTICES;i++){
-						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
 						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
 						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
@@ -4936,5 +4934,5 @@
 
 		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
-		Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+		Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
 	}
 
@@ -5193,5 +5191,5 @@
 
 	/* compute VelocityFactor */
-	VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
+	VelocityFactor= n_man*CR*CR*rho_water*g/mu_water;
 	
 	gauss=new GaussTria();
@@ -5205,9 +5203,8 @@
 
 		/* Water velocity x and y components */
-	//	vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-	//	vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
-	
-		vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-		vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+	//	vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+	//	vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+		vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+		vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
 	}
 
@@ -5304,5 +5301,5 @@
 
 		/*Artificial diffusivity*/
-		vel=sqrt(pow(vx,2.)+pow(vy,2.));
+		vel=sqrt(vx*vx+vy*vy);
 		K[0][0]=diffusivity*h/(2*vel)*vx*vx;
 		K[1][0]=diffusivity*h/(2*vel)*vy*vx;
@@ -5414,9 +5411,7 @@
 
 	/*Intermediaries*/
-	const int numdof = NDOF1*NUMVERTICES;
-
-	int       i;
-	int*      doflist=NULL;
-	IssmDouble    values[numdof];
+	const int   numdof         = NDOF1 *NUMVERTICES;
+	int        *doflist        = NULL;
+	IssmDouble  values[numdof];
 
 	/*Get dof list: */
@@ -5424,9 +5419,8 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++){
+	for(int i=0;i<numdof;i++){
 		values[i]=solution[doflist[i]];
 		if(xIsNan<IssmDouble>(values[i])) _error_("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
- 
+		if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
 	}
 
@@ -5679,5 +5673,5 @@
 		if(stabilization==1){
 			/*Streamline upwinding*/
-			vel=sqrt(pow(vx,2.)+pow(vy,2.));
+			vel=sqrt(vx*vx+vy*vy);
 			K[0][0]=h/(2*vel)*vx*vx;
 			K[1][0]=h/(2*vel)*vy*vx;