Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26248)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26249)
@@ -164,5 +164,4 @@
 
 	if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
-		//Otherwise return the elastic value
 		IssmDouble ka=la0 + 2.0/3.0*mu0; //Bulk modulus
 		if (rheo==2){//EBM
@@ -203,5 +202,5 @@
 
 			mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega);
-			la=ka-2.0/3.0*mu0;
+			la=ka-2.0/3.0*mu;
 		}
 		else{//Maxwell
@@ -210,5 +209,8 @@
 		}
 	}
-
+	else{//Otherwise return the elastic value
+	la=la0;
+	mu=mu0;
+	}
 	*pla=la;
 	*pmu=mu;
@@ -444,5 +446,5 @@
 			   dy[1*6+1]= -2.0/x-fgr/rg0;
 			   }
-			   */ /*}}}*/
+	*/ /*}}}*/
 
 	if(issolid==1){
@@ -453,12 +455,13 @@
 
 	for (id=0;id<ny;id++){
+		dydx[id]=0.0;
 		for (iy=0;iy<ny;iy++){
 			nindex=layer_index*nstep*36+n*36;
-			dydx[id]=dydx[id]+yi_prefactor[nindex+id*6+iy]*y[iy];
+			dydx[id]+=yi_prefactor[nindex+id*6+iy]*y[iy];
 		}
 	}
 	return;
 }/*}}}*/
-template <typename doubletype> void        propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+template <typename doubletype> void        propagate_yi_euler(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
 	//yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
 	//euler method
@@ -466,5 +469,5 @@
 	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
 
-	IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
+	doubletype* dydx=xNewZeroInit<doubletype>(6);
 	IssmDouble dr = (xmax -xmin)/nstep;
 	IssmDouble x=xmin;
@@ -476,5 +479,5 @@
 		x = x + dr;
 	}
-	xDelete<IssmDouble>(dydx);
+	xDelete<doubletype>(dydx);
 }/*}}}*/
 template <typename doubletype> void        propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
@@ -681,5 +684,5 @@
 		if (matlitho->issolid[i]){
 			one = -1.0;
-			if (i>0 && !matlitho->issolid[i-1]) one = 1.0;
+			if (i>0) if (!matlitho->issolid[i-1]) one = 1.0;
 			for (int j=0;j<6;j++){
 				yi[(j+6*ici)+ nyi*(j+6*ici+3)] = one;
@@ -891,7 +894,7 @@
 		if (loveratio<=underflow_tol || xIsNan(loveratio) || xIsInf(loveratio)){//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(VerboseSolution()){
-				printf("\n   Degree: %i, surface/Depth Love number ratio small: %e\n",deg, loveratio);
-				printf("    Changing the interface where the integration starts\n");
-				printf("    New start interface: r=%gm\n\n", matlitho->radius[starting_layer+1]);
+				_printf_("\n   Degree: " << deg <<", surface/Depth Love number ratio small: "<<loveratio<<"\n");
+				_printf_("    Changing the interface where the integration starts\n");
+				_printf_("    New start interface: r="<< matlitho->radius[starting_layer+1] <<"m\n\n");
 			}
 			nyi-=6;
@@ -1005,5 +1008,5 @@
 	bool        allow_layer_deletion,love_kernels, istemporal;
 	bool        verbosemod = (int)VerboseModule();
-	IssmDouble  g0,r0,mu0, underflow_tol;
+	IssmDouble  mu0, underflow_tol;
 	IssmDouble *frequencies = NULL;
 	bool        save_results;
@@ -1040,6 +1043,4 @@
 	femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
 	femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
-	femmodel->parameters->FindParam(&g0,LoveG0Enum);
-	femmodel->parameters->FindParam(&r0,LoveR0Enum);
 	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
 	femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
@@ -1126,13 +1127,47 @@
 	femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
 
+	/*Only when love_kernels is on*/
+	if (love_kernels==1) {
+		femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
+	}
+
+	/*Free resources:*/
+	xDelete<IssmDouble>(frequencies);
+	xDelete<doubletype>(LoveKf);
+	xDelete<doubletype>(LoveHf);
+	xDelete<doubletype>(LoveLf);
+	xDelete<doubletype>(LoveKernels);
+
+	/* Legacy for fortran core, to be removed after complete validation */
+
+	//IssmDouble g0,r0;
+	//femmodel->parameters->FindParam(&g0,LoveG0Enum);
+	//femmodel->parameters->FindParam(&r0,LoveR0Enum);
+	//IssmDouble* LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	//IssmDouble* LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	//IssmDouble* LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	//IssmDouble* LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	//IssmDouble* LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	//IssmDouble* LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+
+	/*Initialize love kernels (real and imaginary parts): */
+	//IssmDouble* LoveKernelsr= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
+	//IssmDouble* LoveKernelsi= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
 
 	/*call the main module: */
 	//if (false){
-	//FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
+	//FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHr,LoveLr,LoveLi,LoveKernelsr,LoveKernelsi,  //output
 	//		nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
 	//		matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
 	//		matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->rheologymodel, matlitho->issolid //matlitho inputs
 	//		);
+
+
+	/*Only when love_kernels is on*/
+	//if (love_kernels==1) {
+		//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernelsr,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
+		//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
 	//}
+
 
 	/*Add love matrices to results:*/
@@ -1144,17 +1179,15 @@
 	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
 
-	/*Only when love_kernels is on*/
-	if (love_kernels==1) {
-		femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
-		//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
-	}
-
-	/*Free resources:*/
-	xDelete<IssmDouble>(frequencies);
-	xDelete<doubletype>(LoveKf);
-	xDelete<doubletype>(LoveHf);
-	xDelete<doubletype>(LoveLf);
+
+
+	//xDelete<IssmDouble>(LoveKr);
+	//xDelete<IssmDouble>(LoveHr);
+	//xDelete<IssmDouble>(LoveLr);
+	//xDelete<IssmDouble>(LoveKernelsr);
+	//xDelete<IssmDouble>(LoveKi);
+	//xDelete<IssmDouble>(LoveHi);
 	//xDelete<IssmDouble>(LoveLi);
-	xDelete<doubletype>(LoveKernels);
+	//xDelete<IssmDouble>(LoveKernelsi);
+
 } /*}}}*/
 
@@ -1173,4 +1206,6 @@
 template void        propagate_yi_RK2<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
 template void        propagate_yi_RK2<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_euler<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_euler<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
 template void        Innersphere_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
 template void        Innersphere_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);