Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 28091)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 28092)
@@ -352,5 +352,5 @@
 	//}
 
-	if (PW_test==0){ //elastic or fluid response: Love(t) = Love(s), we can do an early return
+	if (PW_test==0.0){ //elastic or fluid response: Love(t) = Love(s), we can do an early return
 		Lovet[t*(sh_nmax+1)+d]=Lovef[indf];
 		return;
@@ -384,5 +384,5 @@
 	doubletype dalpha=1.0/(nalpha-1); // alpha table resolution given 0 <= alpha <= 1
 	ialpha= static_cast<int>(DownCastVarToDouble(alpha/dalpha));
-	lincoef=alpha/dalpha-ialpha;//linear fraction in [0;1] for alpha interpolation
+	lincoef=alpha/dalpha-reCast<doubletype>(ialpha);//linear fraction in [0;1] for alpha interpolation
 	iz1=nz;
 	for (int i=0;i<nz;i++){
@@ -439,7 +439,7 @@
 		
 		//cubic spline functions
-		h00=2*pow(t,3)-3*pow(t,2)+1;
-		h10=pow(t,3)-2*pow(t,2)+t;
-		h01=-2*pow(t,3)+3*pow(t,2);
+		h00=2.0*pow(t,3)-3.0*pow(t,2)+1.0;
+		h10=pow(t,3)-2.0*pow(t,2)+t;
+		h01=-2.0*pow(t,3)+3.0*pow(t,2);
 		h11=pow(t,3)-pow(t,2);
 
@@ -589,5 +589,5 @@
 	int rheo=matlitho->rheologymodel[layer_index];
 
-	if(vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
+	if(vi!=0.0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
 		doubletype ka=la0 + 2.0/3.0*mu0; //Bulk modulus
 		if(rheo==2){//EBM
@@ -717,9 +717,9 @@
 				xmin=matlitho->radius[layer_index]/ra;
 				xmax=(matlitho->radius[layer_index+1])/ra;
-				dr = (xmax -xmin)/nstep;
+				dr = (xmax -xmin)/reCast<doubletype>(nstep);
 				x=xmin;
 
 				//fixme
-				g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
+				g=GetGravity<doubletype>((xmin+xmax)/2.0*ra,layer_index,femmodel,matlitho,vars);
 
 				for (int n=0;n<nstep;n++){
@@ -759,10 +759,10 @@
 			xmin=matlitho->radius[layer_index]/ra;
 			xmax=(matlitho->radius[layer_index+1])/ra;
-			dr = (xmax -xmin)/nstep;
+			dr = (xmax -xmin)/reCast<doubletype>(nstep);
 			x=xmin;
 
 
 				//fixme
-				g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
+				g=GetGravity<doubletype>((xmin+xmax)/2.0*ra,layer_index,femmodel,matlitho,vars);
 
 			for (int n=0;n<nstep;n++){
@@ -796,8 +796,8 @@
 			xmin=matlitho->radius[layer_index]/ra;
 			xmax=(matlitho->radius[layer_index+1])/ra;
-			dr = (xmax -xmin)/nstep;
+			dr = (xmax -xmin)/reCast<doubletype>(nstep);
 			x=xmin;
 				//fixme
-				g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
+				g=GetGravity<doubletype>((xmin+xmax)/2.0*ra,layer_index,femmodel,matlitho,vars);
 			for (int n=0;n<nstep;n++){
 				g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
@@ -923,5 +923,5 @@
 
 	doubletype* dydx=xNewZeroInit<doubletype>(6);
-	doubletype dr = (xmax -xmin)/nstep;
+	doubletype dr = (xmax -xmin)/reCast<doubletype>(nstep);
 	doubletype x=xmin;
 	for(int i = 0;i<nstep;i++){
@@ -949,5 +949,5 @@
 	doubletype y3[6]={0};
 
-	doubletype dr = (xmax -xmin)/nstep;
+	doubletype dr = (xmax -xmin)/reCast<doubletype>(nstep);
 	doubletype x=xmin;
 
@@ -978,5 +978,5 @@
 	doubletype y3[6]={0};
 
-	doubletype dr = (xmax -xmin)/nstep;
+	doubletype dr = (xmax -xmin)/reCast<doubletype>(nstep);
 	doubletype x=xmin;
 	for(int i = 0;i<nstep/2-1;i++){
@@ -1066,9 +1066,9 @@
 
 	yi[2+nyi*0]=(deg+3.0)/(deg*(deg+1.0))*r/ra;
-	yi[2+nyi*1]=1.0/(deg*r*ra);
+	yi[2+nyi*1]=1.0/(reCast<doubletype>(deg)*r*ra);
 	yi[2+nyi*2]=0.0;
 
-	yi[3+nyi*0]=2.0*mu*(deg+2.0)/((deg+1.0)*mu0);
-	yi[3+nyi*1]=2.0*mu*(deg-1.0)/(deg*r2*mu0);
+	yi[3+nyi*0]=2.0*mu*(reCast<doubletype>(deg)+2.0)/((deg+1.0)*mu0);
+	yi[3+nyi*1]=2.0*mu*(deg-1.0)/(reCast<doubletype>(deg)*r2*mu0);
 	yi[3+nyi*2]=0.0;
 
@@ -1079,5 +1079,5 @@
 	yi[5+nyi*0]=-cst*r/g0;
 	yi[5+nyi*1]=-cst/(r*g0);
-	yi[5+nyi*2]=deg/(r*g0);
+	yi[5+nyi*2]=reCast<doubletype>(deg)/(r*g0);
 
 
@@ -1326,5 +1326,5 @@
 		//-- forcings at the Inner Core Boundary
 		case 1:	//'ICB --Volumetric Potential'
-			yi_righthandside[6*icb+5]=(deg)/(rc*g0);
+			yi_righthandside[6*icb+5]=reCast<doubletype>(deg)/(rc*g0);
 			yi_righthandside[6*icb+4]=1.0/(ra*g0);
 			break;
@@ -1487,5 +1487,5 @@
 		}
 
-		if (omega==0){ // if running elastic love_numbers, record at which degree we must delete layers, this way we synch layer deletion between cpus next time we calculate love numbers
+		if (omega==0.0){ // if running elastic love_numbers, record at which degree we must delete layers, this way we synch layer deletion between cpus next time we calculate love numbers
 			//We need to delete a layer and try again if the ratio between deepest love number to surface love number is too low (risk of underflow) or garbage
 			if (loveratio<=underflow_tol || xIsNan(loveratio) || xIsInf(loveratio)) {