Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26235)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26236)
@@ -9,20 +9,976 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
-void love_core(FemModel* femmodel){
-
-	Vector<IssmDouble> *wg    = NULL;
-	Vector<IssmDouble> *dwdtg = NULL;
-	IssmDouble          *x    = NULL;
-	IssmDouble          *y    = NULL;
+#include "petscblaslapack.h"
+
+
+void love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/
+	// initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system
+	bool        verbosemod = (int)VerboseModule();
+
+	int numlayers = matlitho->numlayers;
+	IssmDouble* r=matlitho->radius;
+	IssmDouble* EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
+	IssmDouble r1,r2,ro, GG;
+	femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+
+	for (int i=0;i<numlayers;i++){
+		  r2 = r[i+1];
+		  ro = matlitho->density[i];
+		if (i==0){
+			EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0;
+		}else{
+			r1=r[i];
+			EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;;
+		}
+	}
+
+	IssmDouble g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0);
+	femmodel->parameters->SetParam(g0,LoveG0Enum);
+	femmodel->parameters->SetParam(EarthMass,numlayers,LoveEarthMassEnum);
+	xDelete<IssmDouble>(EarthMass);
+
+	int nyi=6*(numlayers+1);
+	int starting_layer=0;
+	femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
+}
+
+void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+
+	//returns lame parameters (material rigity) lambda and mu for the right frequency and layer
+ 	IssmDouble vi=matlitho->viscosity[layer_index];
+ 	IssmDouble mu=matlitho->lame_mu[layer_index];
+ 	IssmDouble la=matlitho->lame_lambda[layer_index];
+ 	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
+				   //Otherwise return the elastic value
+		IssmDouble ka=la + 2.0/3.0*mu; //Bulk modulus
+		if (rheo==2){//EBM
+			IssmDouble alpha=matlitho->ebm_alpha[layer_index];
+			IssmDouble delta=matlitho->ebm_delta[layer_index];
+			IssmDouble taul=matlitho->ebm_taul[layer_index];
+			IssmDouble tauh=matlitho->ebm_tauh[layer_index];
+			IssmDouble hf1,hf2,U1,U2;
+			IssmDouble* a=xNew<IssmDouble>(2);
+			IssmDouble* b=xNew<IssmDouble>(1);
+
+			a[0]=1.0;a[1]=1.0+alpha/2.0;
+			b[0]=2.0+alpha/2.0;
+			//hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
+			//hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
+			//hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
+			//hf2=hypergeometric_pFq_({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
+			U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf1-pow(taul,2.0+alpha)*hf2);
+
+			a[0]=1.0;a[1]=.5+alpha/2.0;
+			b[0]=1.5+alpha/2.0;
+			//hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
+			//hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
+			//hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
+			//hf2=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
+			U2=(pow(tauh,1.0+alpha)*hf1-pow(taul,1.0+alpha)*hf2)/(1.0+alpha);
+
+			mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega);
+		   	la=ka-2.0/3.0*mu;
+
+			xDelete<IssmDouble>(a);
+			xDelete<IssmDouble>(b);
+		} 
+		else if (rheo==1){//Burgers
+			IssmDouble vi2=matlitho->burgers_viscosity[layer_index];
+			IssmDouble mu2=matlitho->burgers_mu[layer_index];
+
+		  	mu=mu*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu/vi)+mu/vi2*omega);
+		   	la=ka-2.0/3.0*mu;
+	  	}
+		else{//Maxwell
+			la = (la + mu*ka/vi/omega)/(1.0 + mu/vi/omega);
+			mu = mu/(1.0+mu/vi/omega);
+		}
+	 }
+
+	 *pla=la;
+	 *pmu=mu;
+
+} /*}}}*/
+
+
+IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+	//computes gravity at radius r2
+	IssmDouble* EarthMass;
+	IssmDouble g, GG;
+
+	femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
+	femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+	IssmDouble ro=matlitho->density[layer_index];
+	IssmDouble M=0;
+	IssmDouble r1=0;
+	if (layer_index==0){
+		M=4.0/3.0*PI*ro*pow(r2,3.0);
+	}
+	else{ 
+		r1=matlitho->radius[layer_index];
+		M=EarthMass[layer_index-1]+4.0/3.0*PI*ro*(pow(r2,3.0)-pow(r1,3.0));
+	}
+	return	g= GG*M/pow(r2,2.0);
+}/*}}}*/
+
+void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+	//precalculates partial derivative factors for function yi_derivatives
+ IssmDouble ra=matlitho->radius[matlitho->numlayers];
+ IssmDouble  g0,r0,mu0, GG;
+ int nstep,nindex, starting_layer;
+ femmodel->parameters->FindParam(&g0,LoveG0Enum);
+ femmodel->parameters->FindParam(&r0,LoveR0Enum);
+ femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+ femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+ femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
+ femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
+
+ IssmDouble frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
+ IssmDouble xmin,xmax,x,dr;
+ IssmDouble g,ro, issolid;
+
+	if (pomega) { //frequency and degree dependent terms /*{{{*/
+		IssmDouble la,mu;
+		IssmDouble* f=xNewZeroInit<IssmDouble>(12);
+		int deg=*pdeg;
+		IssmDouble omega=*pomega;	
+		fn=deg*(deg+1.0);
+
+		for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
+
+				ro=matlitho->density[layer_index];
+				issolid=matlitho->issolid[layer_index];
+				if(issolid==1){
+		 		GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
+	
+				/*_______Expressions*/
+				flm=(la+2.0*mu);
+				rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
+				rm0=mu/mu0;
+				frh=ro*ra/mu0;
+
+				f[0]=(-2.0*la/flm);
+				f[1]=mu0/flm;
+				f[2]=(la*fn/flm);
+				f[3]=rm0*rlm;
+				f[4]=ro*pow(omega,2.0)*ra/mu0;
+				f[5]=(-4.0*mu/flm);
+				f[6]=fn*frh;
+				f[7]=-(2.0*rm0*rlm)*fn;
+				f[8]=1.0/rm0;
+				f[9]=-2.0*rm0*rlm;
+				f[10]=-la/flm;
+				f[11]=2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm;
+
+				xmin=matlitho->radius[layer_index]/ra;
+				xmax=(matlitho->radius[layer_index+1])/ra;
+				dr = (xmax -xmin)/nstep;
+			     	x=xmin;
+				for (int n=0;n<nstep;n++){
+					g=GetGravity(x*ra,layer_index,femmodel,matlitho);
+
+					nindex=layer_index*nstep*36+n*36;
+					yi_prefactor[nindex+ 0*6+0]= f[0]/x;                      // in dy[0*6+0]
+					yi_prefactor[nindex+ 0*6+1]= f[1];                        // in dy[0*6+1]
+					yi_prefactor[nindex+ 0*6+2]= f[2]/x;                      // in dy[0*6+2]
+					yi_prefactor[nindex+ 1*6+0]= 4.0*(-frh*g+f[3]/x)/x + f[4];// in dy[1*6+0]
+					yi_prefactor[nindex+ 1*6+1]= f[5]/x;                      // in dy[1*6+1]
+					yi_prefactor[nindex+ 1*6+2]= (f[6]*g+f[7]/x)/x;           // in dy[1*6+2]
+					yi_prefactor[nindex+ 2*6+3]= f[8];                        // in dy[2*6+3]
+					yi_prefactor[nindex+ 3*6+0]= (frh*g+f[9]/x)/x;            // in dy[3*6+0]
+					yi_prefactor[nindex+ 3*6+1]= f[10]/x;                     // in dy[3*6+1]
+					yi_prefactor[nindex+ 3*6+2]= f[11]/(x*x) + f[4];          // in dy[3*6+2]
+					x=x+dr;
+				}
+			}
+		}
+
+		xDelete<IssmDouble>(f);
+	/*}}}*/
+	} else if (pdeg) { // degree dependent terms /*{{{*/
+	 	int deg=*pdeg;
+	 	fn=(deg*(deg+1.0));
+
+		for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
+			ro=matlitho->density[layer_index];
+			issolid=matlitho->issolid[layer_index];
+	
+			/*_______Expressions*/
+			fgr=4.0*PI*GG*ro*ra;
+	 
+			xmin=matlitho->radius[layer_index]/ra;
+			xmax=(matlitho->radius[layer_index+1])/ra;
+			dr = (xmax -xmin)/nstep;
+		     	x=xmin;
+			for (int n=0;n<nstep;n++){
+				 nindex=layer_index*nstep*36+n*36;
+				 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
+	
+				 if(issolid==1){
+			 	 yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
+				 yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
+				 yi_prefactor[nindex+ 5*6+4]= fn/(x*x);		     // in dy[5*6+4]
+				 } else {
+				 yi_prefactor[nindex+ 1*6+0]= (-4.0*(fgr/g)+fn/x)/x; // in dy[1*6+0] liquid layer
+				 }
+				x=x+dr;
+			}
+		}
+	/*}}}*/
+	} else { // static terms /*{{{*/
+		for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
+			ro=matlitho->density[layer_index];
+			issolid=matlitho->issolid[layer_index];
+	
+			/*_______Expressions*/
+			frhg0=ro*g0*ra/mu0;
+			fgr=4.0*PI*GG*ro*ra;
+
+			xmin=matlitho->radius[layer_index]/ra;
+			xmax=(matlitho->radius[layer_index+1])/ra;
+			dr = (xmax -xmin)/nstep;
+		     	x=xmin;
+	  	 	for (int n=0;n<nstep;n++){
+				 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
+				 nindex=layer_index*nstep*36+n*36;
+				 if(issolid==1){
+					 yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
+					 yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
+					 yi_prefactor[nindex+ 2*6+2]= 1.0/x;        // in dy[2*6+2]
+					 yi_prefactor[nindex+ 3*6+3]= -3.0/x;       // in dy[3*6+3]
+					 yi_prefactor[nindex+ 3*6+4]= -frhg0/x;     // in dy[3*6+4]
+					 yi_prefactor[nindex+ 4*6+0]= fgr/g0;       // in dy[4*6+0]
+					 yi_prefactor[nindex+ 4*6+5]= 1.0;          // in dy[4*6+5]
+					 yi_prefactor[nindex+ 5*6+5]= -2.0/x;       // in dy[5*6+5]
+				 } else {
+					 yi_prefactor[nindex+ 0*6+0]= fgr/g;        // in dy[0*6+0] liquid layer
+					 yi_prefactor[nindex+ 0*6+1]= 1.0;          // in dy[0*6+1] liquid layer
+					 yi_prefactor[nindex+ 1*6+1]= -2.0/x-fgr/g; // in dy[1*6+1] liquid layer
+				 }
+				x=x+dr;
+			}
+		}
+	/*}}}*/
+	}
+}/*}}}*/
+
+IssmDouble* yi_derivatives(IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+	//computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
+
+	IssmDouble issolid=matlitho->issolid[layer_index];
+	int iy,id,ny, nindex, nstep;
+ 	IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* dy=xNewZeroInit<IssmDouble>(36);
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
+
+	/*{{{*/ /* For reference:
+	 flm=(la+2.0*mu);
+	 rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
+	 rm0=mu/mu0;
+	 rg0=g/g0;
+	 frh=ro*g*ra/mu0;
+	 fgr=4.0*PI*GG*ro*ra/g0;
+	 fn=(deg*(deg+1.0));
+	 
+
+	 if(issolid==1){
+	 ny = 6;
+
+	 dy[0*6+0]= (-2.0*la/flm)/x;
+	 dy[0*6+1]= mu0/flm;
+	 dy[0*6+2]= (la*fn/flm)/x;
+	 dy[0*6+3]= 0.0;
+	 dy[0*6+4]= 0.0;
+	 dy[0*6+5]= 0.0;
+
+	 dy[1*6+0]=  4.0*(-frh+rm0*rlm/x)/x + ro*pow(omega,2.0)*ra/mu0;
+	 dy[1*6+1]=(-4.0*mu/flm)/x;
+	 dy[1*6+2]= fn*(frh-2.0*rm0*rlm/x)/x;
+	 dy[1*6+3]= fn/x;
+	 dy[1*6+4]= 0.0;
+	 dy[1*6+5]= -frh/rg0;
+
+	 dy[2*6+0]= -1.0/x;
+	 dy[2*6+1]= 0.0;
+	 dy[2*6+2]= 1.0/x;
+	 dy[2*6+3]= 1/rm0;
+	 dy[2*6+4]= 0.0;
+	 dy[2*6+5]= 0.0;
+
+	 dy[3*6+0]= (frh-2.0*rm0*rlm/x)/x;
+	 dy[3*6+1]= ( -la/flm)/x;
+	 dy[3*6+2]= (2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm)/(x*x) + ro*pow(omega,2.0)*ra/mu0;
+	 dy[3*6+3]= -3.0/x;
+	 dy[3*6+4]= -(frh/rg0)/x;
+	 dy[3*6+5]= 0.0;
+
+	 dy[4*6+0]= fgr;
+	 dy[4*6+1]= 0.0;
+	 dy[4*6+2]= 0.0;
+	 dy[4*6+3]= 0.0;
+	 dy[4*6+4]= 0.0;
+	 dy[4*6+5]= 1.0;
+
+	 dy[5*6+0]= 0.0;
+	 dy[5*6+1]= 0.0;
+	 dy[5*6+2]= -(fgr*fn)/x;
+	 dy[5*6+3]= 0.0;
+	 dy[5*6+4]= fn/(x*x);
+	 dy[5*6+5]= -2.0/x;
+	  
+	 } else {
+	 ny = 2;
+
+	 dy[0*6+0]= fgr/rg0;
+	 dy[0*6+1]= 1.0;
+	 dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x;
+	 dy[1*6+1]= -2.0/x-fgr/rg0;
+	}
+	*/ /*}}}*/
+
+	 if(issolid==1){
+		 ny = 6;
+	 } else {
+		 ny = 2;
+	 }
+	 
+	 for (id=0;id<ny;id++){
+		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];
+	 	}
+	 }
+	xDelete<IssmDouble>(dy);
+	return dydx;
+}/*}}}*/
+
+void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* 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
+	int nstep;
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
+
+ 	IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
+	IssmDouble dr = (xmax -xmin)/nstep;
+     	IssmDouble x=xmin;
+    for(int i = 0;i<nstep;i++){
+	dydx=yi_derivatives(y,layer_index, i,yi_prefactor,femmodel,matlitho);
+  	for (int j=0;j<6;j++){
+		y[j]+=dydx[j]*dr;
+	}
+  	x = x + dr;
+    }
+	xDelete<IssmDouble>(dydx);
+}/*}}}*/
+
+void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* 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?
+	//Implements Runge-Kutta 2nd order (midpoint method)
+	int nstep;
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
+
+ 	IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
+
+	IssmDouble dr = (xmax -xmin)/nstep;
+     	IssmDouble x=xmin;
+    for(int i = 0;i<nstep/2;i++){
+	k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
+	k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);		
+	
+  	for (int j=0;j<6;j++){
+		  y[j]+=k2[j]*2.0*dr;
+	}
+  	x = x + 2.0*dr;
+    }
+	xDelete<IssmDouble>(k1);
+	xDelete<IssmDouble>(k2);
+	xDelete<IssmDouble>(k3);
+	xDelete<IssmDouble>(k4);
+	xDelete<IssmDouble>(y1);
+	xDelete<IssmDouble>(y2);
+	xDelete<IssmDouble>(y3);
+}/*}}}*/
+
+void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* 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?
+	//Implements Runge-Kutta 4th order
+	int nstep;
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
+
+ 	IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
+ 	IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
+
+	IssmDouble dr = (xmax -xmin)/nstep;
+     	IssmDouble x=xmin;
+    for(int i = 0;i<nstep/2-1;i++){
+	k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
+	k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
+	k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
+	k4=yi_derivatives(y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho);		
+	
+  	for (int j=0;j<6;j++){
+		y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;		
+	}
+  	x = x + 2.0*dr;
+    }
+
+	//Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr
+	int i=nstep/2;
+	k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
+	k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
+	k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+	for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
+	k4=yi_derivatives(y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);		
+	
+  	for (int j=0;j<6;j++){
+		y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;		
+	}
+
+  	x = x + 2.0*dr;
+
+	xDelete<IssmDouble>(k1);
+	xDelete<IssmDouble>(k2);
+	xDelete<IssmDouble>(k3);
+	xDelete<IssmDouble>(k4);
+	xDelete<IssmDouble>(y1);
+	xDelete<IssmDouble>(y2);
+	xDelete<IssmDouble>(y3);
+}/*}}}*/
+
+void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+	//fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5]
+
+	IssmDouble r = matlitho->radius[layer_index];
+	IssmDouble ra=matlitho->radius[matlitho->numlayers];
+	IssmDouble  g0,r0,mu0, GG;
+	int nyi;
+ 	femmodel->parameters->FindParam(&g0,LoveG0Enum);
+ 	femmodel->parameters->FindParam(&r0,LoveR0Enum);
+ 	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+ 	femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+	IssmDouble ro=matlitho->density[layer_index];
+	IssmDouble issolid=matlitho->issolid[layer_index];
+	IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho);
+	IssmDouble la,mu;
+
+	if (layer_index==0){
+		// radius[0] cannot be 0 for numerical reasons, but below our first interface at radius[0] would in reality be the same material as in the first layer
+		GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
+	} else {
+		GetEarthRheology(&la, &mu,layer_index-1,omega,femmodel,matlitho);   
+	}    
+
+	IssmDouble cst = 4.0*PI*GG*ro;
+	IssmDouble r2=pow(r,2.0);
+
+	yi[0+nyi*0]=1.0*r/ra;
+	yi[0+nyi*1]=1.0/(r*ra);
+	yi[0+nyi*2]=0.0;
+
+	yi[1+nyi*0]=(2.0*mu*(deg-1.0-3.0/deg) + cst/3.0*ro*r2)/mu0;
+	yi[1+nyi*1]=(2.0*mu*(deg-1.0)/r2 + cst/3.0*ro)/mu0;
+	yi[1+nyi*2]=-ro/mu0;
+
+	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*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*2]=0.0;
+
+	yi[4+nyi*0]=0.0;
+	yi[4+nyi*1]=0.0;
+	yi[4+nyi*2]=1.0/(g0*ra);
+
+	yi[5+nyi*0]=-cst*r/g0;
+	yi[5+nyi*1]=-cst/(r*g0);
+	yi[5+nyi*2]=deg/(r*g0);
+}/*}}}*/
+
+void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho) { /*{{{*/
+
+	IssmDouble  g0,r0,mu0,x,ro1, GG;
+	int nyi,starting_layer, nstep;
+	femmodel->parameters->FindParam(&g0,LoveG0Enum);
+	femmodel->parameters->FindParam(&r0,LoveR0Enum);
+	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+	femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+	femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
+
+	  IssmDouble xmin,xmax,one,ro,g;
+	  IssmDouble ra=matlitho->radius[matlitho->numlayers];
+
+	for (int i=0;i<6*(matlitho->numlayers+1);i++){
+		for(int j=0;j<6*(matlitho->numlayers+1);j++){
+			yi[i+6*(matlitho->numlayers+1)*j]=0.0;
+		}
+	}
+
+	int ny,is,ii,jj;
+	IssmDouble ystart[6];
+	for (int k=0;k<6;k++) ystart[k]=0.0;		
+  
+	int ici = 0;   // Index of current interface 
+	for (int i = starting_layer; i<matlitho->numlayers;i++){ 
+
+		xmin=matlitho->radius[i]/ra;
+		xmax=(matlitho->radius[i+1])/ra;
+
+		  
+		if (matlitho->issolid[i]){
+		   		ny = 6;
+		   		is = 0;
+		   		one= 1.0;
+		} else {	
+		   		ny = 2;
+		   		is = 4;
+		   		one= -1.0;
+		}
+	
+
+		for (int j = 0;j<ny;j++){
+			for (int k=0;k<6;k++){ystart[k]=0.0;}
+			ystart[j]= 1.0;
+			  		
+			// Numerical Integration 
+			//propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
+			propagate_yi_RK2(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
+			//propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
+			// Boundary Condition matrix - propagation part 
+			ii = 6*(ici+1)+is;
+			jj = 6*(ici+1)+j+is-3;
+			for (int kk=0;kk<ny;kk++){
+		  		yi[(ii+kk)+nyi*jj] = ystart[kk]*one;
+			}
+		}
+		  
+		  
+		// Boundary Condition matrix - solid regions
+		if (matlitho->issolid[i]){
+			one = -1.0;
+			if (i>0 & !matlitho->issolid[i-1]) one = 1.0;
+			for (int j=0;j<6;j++){
+				yi[(j+6*ici)+ nyi*(j+6*ici+3)] = one;
+			}
+		} else { // Boundary Condition matrix - liquid regions
+			ro1=matlitho->density[i];
+			g=GetGravity(matlitho->radius[i], i, femmodel,matlitho);
+			ii = 6*ici;
+			yi[ii+nyi*(ii+3)] = -1.0;
+			yi[ii+nyi*(ii+4+3)] = -g0/g;
+			yi[(ii+1)+nyi*(ii+3)]=-ro1*g*ra/mu0;
+			yi[(ii+2)+nyi*(ii+1+3)]=-1.0;
+			yi[(ii+5)+nyi*(ii+3)]= 4.0*PI*GG*ro1*ra/g0;
+			yi[(ii+4)+nyi*(ii+4+3)]=-1.0;
+			yi[(ii+5)+nyi*(ii+5+3)]=-1.0;
+			g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho);
+			ii = 6*(ici+1);
+			yi[ii+nyi*(ii-1)]=-1.0;
+			yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB
+			yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB
+			  // yi(13,..) y1 CMB
+			yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0;
+			yi[(ii+2)+nyi*(ii)]=-1.0;
+			yi[(ii+5)+nyi*(ii-1)]= 4.0*PI*GG*ro1*ra/g0;
+		}	
+		ici = ici+1;
+	}
+
+
+  	//-- Internal sphere: integration starts here rather than r=0 for numerical reasons
+	Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho);
+
+  	//-- Surface conditions
+	yi[(nyi-6)+nyi*(nyi-3)]=-1.0;
+	yi[(nyi-4)+nyi*(nyi-2)]=-1.0;
+	yi[(nyi-2)+nyi*(nyi-1)]=-1.0;
+	yi[(nyi-1)+nyi*(nyi-1)]=deg+1.0;
+
+  	//-- Degree 1 special case
+	if(deg==1){
+		for (int i=0;i<nyi;i++){
+			yi[(nyi-1)+nyi*i]=0.0;
+		}
+		yi[(nyi-1)+nyi*(nyi-1)]=1.0;
+	}
+}/*}}}*/
+
+void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+
+	IssmDouble  g0,r0,mu0,ra,rb,rc;
+	int nyi, forcing_type,icb,cmb;
+	IssmDouble* EarthMass;
+	femmodel->parameters->FindParam(&g0,LoveG0Enum);
+	femmodel->parameters->FindParam(&r0,LoveR0Enum);
+	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+	femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
+	femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum);
+	femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
+	femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
+	  
+   	// In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
+	ra=matlitho->radius[matlitho->numlayers];	
+	rb=0;
+	rc=0;
+	if (forcing_type<=4){
+		rc=matlitho->radius[icb];
+	} 
+	else if (forcing_type<=8){
+		rb=matlitho->radius[cmb];
+	}
+
+	IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0));
+	xDelete<IssmDouble>(EarthMass);
+
+	for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0;
+
+	switch (forcing_type) {
+
+	//-- forcings at the Inner Core Boundary
+   	case 1:	//'ICB --Volumetric Potential'
+		yi_righthandside[6*icb+5]=(deg)/(rc*g0);
+		yi_righthandside[6*icb+4]=1.0/(ra*g0);
+		break;
+	case 2: //'ICB --Pressure'
+		yi_righthandside[6*icb+1]=-ro_mean/mu0;
+		break;
+	case 3://'ICB --Loading'
+		yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc;
+		yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0);
+		break;
+	case 4://'ICB --Tangential Traction'
+		yi_righthandside[6*icb+3]= ro_mean/mu0;
+		break;
+
+	//--forcings at the Core Mantle Boundary
+	case 5://'CMB --Volumetric Potential'
+		yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb;
+		yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
+		break;
+	case 6://'CMB --Pressure'
+		yi_righthandside[6*cmb+1]=-ro_mean/mu0;
+		break;
+	case 7://'CMB --Loading'
+		yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb;
+		yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
+		break;
+	case 8://'CMB --Tangential Traction'
+		yi_righthandside[6*cmb+3]=-ro_mean/mu0;
+		break;
+
+	//--forcings at the surface
+	case 9://'SURF--Volumetric Potential'
+		if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0);
+		break;
+	case 10://'SURF--Pressure'
+		yi_righthandside[nyi-5]=-ro_mean/mu0;
+		break;
+	case 11://'SURF--Loading'
+		yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0);
+		if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0);
+		break;
+	case 12://'SURF--Tangential Traction'
+		yi_righthandside[nyi-3]= ro_mean/mu0;
+		break;
+	default:
+		_error_("love core error: forcing_type not supported yet");
+	}
+}/*}}}*/
+
+
+void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+
+	IssmDouble  g0,r0,mu0,loveratio,underflow_tol;
+	IssmDouble* frequencies;
+	int nyi,starting_layer, dummy;
+	bool allow_layer_deletion;
+	femmodel->parameters->FindParam(&g0,LoveG0Enum);
+	femmodel->parameters->FindParam(&r0,LoveR0Enum);
+	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+	femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
+	femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
+	femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum);
+	femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
+ 	IssmDouble ra=matlitho->radius[matlitho->numlayers];
+	bool exit=false;
+	int lda,ldb;
+
+	for(;!exit;){ //cycles of: attempt to solve the yi system, then delete a layer if necessary
+		lda=nyi;
+		ldb=nyi;
+		IssmDouble*  yilocal=xNew<IssmDouble>(nyi*nyi); // we will need to redeclare these inside here as nyi changes
+		IssmDouble*  rhslocal=xNew<IssmDouble>(nyi);
+
+		//we need to do a local copy of yi,rhs to send them to LAPACK with the appropriate size and to keep the original matrices in case layers are deleted and we need to try again
+		for (int i=0;i<nyi;i++){ 
+			rhslocal[i]=rhs[i];
+			for (int j=0;j<nyi;j++){
+				yilocal[i+j*nyi]=yi[i+j*nyi];
+			}
+		}
+		//-- Resolution
+		int* ipiv=xNewZeroInit<int>(nyi); //pivot index vector
+		int info = 0;// error checker
+		int nrhs=1; // number of right hand size columns
+
+		dgesv_(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info);
+		xDelete<int>(ipiv);
+
+		if(VerboseSolution() && info!=0){ 
+			_printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
+			printf("%s %s %s %s\n", "i","j","yi[i+nyi*j]","rhs[i]");
+			for (int i=0;i<nyi;i++){
+				for (int j=0;j<nyi;j++){
+					printf("%i %i %e %e\n", i,j,yi[i+nyi*j],rhs[i]);
+				}
+			}
+		}
+
+		*loveh = rhslocal[nyi-3]*ra*g0;
+		*lovel = rhslocal[nyi-2]*ra*g0;
+		*lovek = rhslocal[nyi-1]*ra*g0;
+
+		IssmDouble loveh1 = rhslocal[3];
+		IssmDouble lovel1 = rhslocal[5];
+		IssmDouble lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra);
+
+		IssmDouble loveh1s = rhslocal[nyi-3];
+		IssmDouble lovel1s = rhslocal[nyi-2];
+		IssmDouble lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra);
+
+		loveratio = abs(loveh1/loveh1s); //ratio of center to surface love numbers, determines if we should remove layers
+		if (abs(lovel1/lovel1s) < loveratio) loveratio = abs(lovel1/lovel1s); 
+		if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s);
+
+		if (!allow_layer_deletion || nyi<=12 || omega!=2.0*PI*frequencies[0]){ //We are not allowed to delete layers, or there is only one layer left. We also don't want to delete layers in the middle of a loop on frequencies, as that can lead to a jump that would compromise the inverse laplace transform.
+			goto save_results;
+		}
+
+		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]);
+			}
+			nyi-=6;
+			starting_layer+=1;
+			femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
+			femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
+
+
+			for (int i=0;i<nyi;i++){//shift everything down by 1 layer
+				rhs[i]=rhs[i+6];
+				for (int j=0;j<nyi;j++){
+					yi[j+i*nyi]=yi[j+6+(i+6)*(nyi+6)];
+				}
+			}
+
+			Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different
+		} else { //we are ready to save the outputs and break the main loop
+
+		save_results:
+			for (int i=0;i<nyi;i++){
+				rhs[i]=rhslocal[i];
+				for (int j=0;j<nyi;j++){
+					yi[j+i*nyi]=yilocal[j+i*nyi];
+				}
+			}
+
+			//make sure we can't output numbers from deleted layers
+			for (int i=nyi;i<(matlitho->numlayers+1)*6;i++){ 
+				rhs[i]=0.0;
+				for (int j=0;j<(matlitho->numlayers+1)*6;j++){
+					yi[j+i*(matlitho->numlayers+1)*6]=0.0;
+				}
+			}
+			for (int i=0;i<nyi;i++){
+				for (int j=nyi;j<(matlitho->numlayers+1)*6;j++){
+					yi[j+i*(matlitho->numlayers+1)*6]=0.0;
+				}
+			}
+
+			exit = true;
+	  	}
+		xDelete<IssmDouble>(yilocal);
+		xDelete<IssmDouble>(rhslocal);
+  	}
+	xDelete<IssmDouble>(frequencies);	
+}/*}}}*/
+
+IssmDouble factorial(int n){ /*{{{*/
+	IssmDouble prod=1;
+	for (int i=2;i<n+1;i++) prod*=i;
+	return prod;
+}/*}}}*/
+
+IssmDouble n_C_r(int n, int r){ /*{{{*/ 
+	//n choose r
+        int primes[169] = 
+  {2,    3,    5,    7,   11,   13,   17,   19,   23,   29,  
+  31,   37,   41,   43,   47,   53,   59,   61,   67,   71,  
+  73,   79,   83,   89,   97,  101,  103,  107,  109,  113,  
+ 127,  131,  137,  139,  149,  151,  157,  163,  167,  173,  
+ 179,  181,  191,  193,  197,  199,  211,  223,  227,  229,  
+ 233,  239,  241,  251,  257,  263,  269,  271,  277,  281,  
+ 283,  293,  307,  311,  313,  317,  331,  337,  347,  349,  
+ 353,  359,  367,  373,  379,  383,  389,  397,  401,  409,  
+ 419,  421,  431,  433,  439,  443,  449,  457,  461,  463,  
+ 467,  479,  487,  491,  499,  503,  509,  521,  523,  541,  
+ 547,  557,  563,  569,  571,  577,  587,  593,  599,  601,  
+ 607,  613,  617,  619,  631,  641,  643,  647,  653,  659,  
+ 661,  673,  677,  683,  691,  701,  709,  719,  727,  733,  
+ 739,  743,  751,  757,  761,  769,  773,  787,  797,  809,  
+ 811,  821,  823,  827,  829,  839,  853,  857,  859,  863,  
+ 877,  881,  883,  887,  907,  911,  919,  929,  937,  941,  
+ 947,  953,  967,  971,  977,  983,  991,  997, 1009};
+	int num, den;
+        num = 1;
+        den = 1;
+
+	for (int i=0;i<r;i++){
+		num = num*(n-i);
+		den = den*(i+1);
+		if (i>0) {
+			// Divide out common prime factors
+			for (int k=0;k<169;k++){ //169 is the length of the prime vector here
+				if ( i % primes[k] == 0) { // modulo
+			                num = num/primes[k];
+			                den = den/primes[k];
+		        	}
+		        }
+		}
+	}
+
+	IssmDouble res;        
+	return res = num/den;
+}/*}}}*/
+
+
+IssmDouble* postwidder_coef(int NTit){ /*{{{*/
+	//Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform:
+	//The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k))
+	//The method is based on equations (2), (6), (7) in: 
+	//Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004)
+	//Note that the coefficients xi lack the factor s=k*log(2)/t. 
+	//That is because we are computing the heaviside response of the system rather than its impulse response, 
+	//and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M.
+
+	IssmDouble* xi=xNewZeroInit<IssmDouble>(2*NTit*NTit);
+	int indxi;
+	for (int M=1;M<NTit+1;M++){
+		for (int k=1;k<2*M+1;k++){
+		indxi=(M-1)*(2*NTit)+k-1;
+		    for (int j=floor((k+1)/2);j<min(k,M)+1;j++){
+			    xi[indxi]+=pow(j,M+1.0)/factorial(M)*n_C_r(M,j)*n_C_r(2*j,j)*n_C_r(j,k-j);
+		    }
+		    xi[indxi]*=pow(-1.0,k+M)/k;
+		}
+	}
+	return xi;
+}/*}}}*/
+
+void postwidder_transform(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel){ /*{{{*/
+	//Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef
+
+	int nfreq, indxi, indf;
+	femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
+	int nt=nfreq/2/NTit;
+
+	indf=d*nfreq+t*2*NTit;
+	IssmDouble LoveM[NTit];
+
+	for (int M=1;M<NTit+1;M++){
+		LoveM[M-1]=0.0;
+		for (int k=1;k<2*M+1;k++){
+			indxi=(M-1)*(2*NTit)+k-1;
+			LoveM[M-1]+=xi[indxi]*Lovef[indf+k-1];
+		}
+		
+		//Make sure we are not getting into numerical instability
+		//Diverging once: ok, we'll give that the benefit of the doubt, it can be an inflexion point in the convergence series
+		//Diverging twice in a row: we are definitely propagating numerical error: revert to the last stable value and exit
+		if (M>3){ 
+			if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) &&
+			     abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){
+				Lovet[d*nt+t]=LoveM[M-3];
+				return;
+			}
+		}
+	}
+	Lovet[d*nt+t]=LoveM[NTit-1];
+}/*}}}*/
+
+void love_freq_to_temporal(IssmDouble* LoveHt,IssmDouble* LoveLt,IssmDouble* LoveKt,IssmDouble* LoveHf,IssmDouble* LoveLf,IssmDouble* LoveKf, FemModel* femmodel){ /*{{{*/
+	//Transforms all frequency-dependent love numbers into time-dependent love numbers
+	int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit;
+	IssmDouble meanh,meanl,meank,dh,dl,dk;
+
+	femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
+	femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
+	femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
+	femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
+	
+	int nt=nfreq/2/NTit;
+	IssmDouble* xi=postwidder_coef(NTit);
+
+	for (int d=sh_nmin;d<sh_nmax+1;d++){
+		for (int t=0;t<nt;t++){
+			postwidder_transform(LoveHt,LoveHf,d,t,NTit,xi,femmodel);
+			postwidder_transform(LoveLt,LoveLf,d,t,NTit,xi,femmodel);
+			postwidder_transform(LoveKt,LoveKf,d,t,NTit,xi,femmodel);
+		}
+	}
+	xDelete<IssmDouble>(xi);
+}/*}}}*/
+
+
+
+void love_core(FemModel* femmodel){ /*{{{*/
+
+	/*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
+	Matlitho* matlitho=NULL;
+	for (Object* & object: femmodel->materials->objects){
+		Material* material=xDynamicCast<Material*>(object);
+		if(material->ObjectEnum()==MatlithoEnum){
+			matlitho=xDynamicCast<Matlitho*>(material);
+			break;
+		}
+	}
+	_assert_(matlitho);
 
 	/*love parameters: */
+	int         nfreq,nyi,starting_layer,nstep,forcing_type,dummy;
+	int         sh_nmin,sh_nmax,kernel_index,deleted_layer_offset;
+	bool        allow_layer_deletion,love_kernels, istemporal;
+	bool        verbosemod = (int)VerboseModule();
+	IssmDouble lovek, loveh, lovel, loveratio;
+	IssmDouble  g0,r0,mu0, underflow_tol, omega;
 	IssmDouble *frequencies = NULL;
-	int         nfreq,dummy;
-	int         sh_nmin,sh_nmax;
-	IssmDouble  g0,r0,mu0;
-	bool        allow_layer_deletion;
-	bool        love_kernels;
-	int         forcing_type;
-	bool        verbosemod = (int)VerboseModule();
 
 	/*parameters: */
@@ -35,4 +991,5 @@
 
 	/*recover love number parameters: */
+
 	femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
 	femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
@@ -45,23 +1002,11 @@
 	femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
 	femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
-
-	/*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
-	Matlitho* matlitho=NULL;
-	for (Object* & object: femmodel->materials->objects){
-		Material* material=xDynamicCast<Material*>(object);
-		if(material->ObjectEnum()==MatlithoEnum){
-			matlitho=xDynamicCast<Matlitho*>(material);
-			break;
-		}
-	}
-	_assert_(matlitho);
+	femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
+	femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
 
 	/*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
-	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));
+	IssmDouble*  LoveKf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	IssmDouble*  LoveHf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+	IssmDouble*  LoveLf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
 
 	/*Initialize love kernels (real and imaginary parts): */
@@ -69,33 +1014,109 @@
 	IssmDouble*  LoveKernelsImag = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
 
+	love_init(femmodel,matlitho);
+	femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+	femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
+	IssmDouble* yi_prefactor=xNewZeroInit<IssmDouble>(6*6*nstep*matlitho->numlayers);
+	IssmDouble* yi_righthandside=xNewZeroInit<IssmDouble>(nyi);
+	IssmDouble* yi=xNewZeroInit<IssmDouble>(nyi*nyi);
+
+	//precalculate yi coefficients that do not depend on degree or frequency
+	fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho); 
+	for(int deg=sh_nmin;deg<sh_nmax+1;deg++){
+
+		//precalculate yi coefficients that depend on degree but not frequency
+		fill_yi_prefactor(yi_prefactor, &deg, NULL,femmodel, matlitho); 
+
+		for (int fr=0;fr<nfreq;fr++){
+			omega=2.0*PI*frequencies[fr]; //angular frequency
+
+			//precalculate yi coefficients that depend on degree and frequency
+			fill_yi_prefactor(yi_prefactor, &deg,&omega,femmodel, matlitho); 
+
+			yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho); 
+
+			build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho);
+
+			solve_yi_system(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho);
+
+			LoveHf[deg*nfreq+fr]=loveh;
+			LoveKf[deg*nfreq+fr]=lovek-1.0;
+			LoveLf[deg*nfreq+fr]=lovel;
+
+			femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); //synchronize nyi in case we deleted a layer
+			deleted_layer_offset=(matlitho->numlayers+1)*6-nyi;// =6 per deleted layer
+			kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 +
+					    fr*(matlitho->numlayers+1)*6 +
+					        deleted_layer_offset;
+			for (int i=0;i<nyi;i++){
+				LoveKernelsReal[kernel_index+i]=yi_righthandside[i];
+			}
+		}
+	}
+
+
+	xDelete<IssmDouble>(yi);
+	xDelete<IssmDouble>(yi_righthandside);
+
+
+	//Temporal love numbers
+	if (istemporal){
+		int NTit;
+		femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
+		int nt = nfreq/2/NTit;
+
+		IssmDouble* LoveHt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
+		IssmDouble* LoveLt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
+		IssmDouble* LoveKt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
+
+		love_freq_to_temporal(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
+
+		femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
+		femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
+		femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0));
+
+		xDelete<IssmDouble>(LoveHt);
+		xDelete<IssmDouble>(LoveLt);
+		xDelete<IssmDouble>(LoveKt);
+	}
+
+	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
+	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));
+	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
+
+
 	/*call the main module: */
-	FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //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->isburgers, matlitho->issolid //matlitho inputs
-			);
+	//if (false){
+	//FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //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
+	//		);
+	//}
 
 	/*Add love matrices to results:*/
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
-	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
-	/*Only when love_kernels is set unity*/
+	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
+	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
+	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
+	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
+	//femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
+	//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<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(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 ressources:*/
+		//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<IssmDouble>(LoveKr);
-	xDelete<IssmDouble>(LoveHr);
-	xDelete<IssmDouble>(LoveLr);
-	xDelete<IssmDouble>(LoveKi);
-	xDelete<IssmDouble>(LoveHi);
-	xDelete<IssmDouble>(LoveLi);
+	xDelete<IssmDouble>(LoveKf);
+	xDelete<IssmDouble>(LoveHf);
+	xDelete<IssmDouble>(LoveLf);
+	//xDelete<IssmDouble>(LoveKi);
+	//xDelete<IssmDouble>(LoveHi);
+	//xDelete<IssmDouble>(LoveLi);
 	xDelete<IssmDouble>(LoveKernelsReal);
 	xDelete<IssmDouble>(LoveKernelsImag);
-}
+} /*}}}*/
+