Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26231)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26232)
@@ -6133,11 +6133,13 @@
 	IssmDouble* tide_love_h  = NULL;
 	IssmDouble* tide_love_k  = NULL;
-	IssmDouble* load_love_k  = NULL;
-	IssmDouble  tide_love_k2secular;
+	IssmDouble* tide_love_l  = NULL;
 	IssmDouble  moi_e, moi_p, omega;
-	IssmDouble Grotm1[3];
-	IssmDouble Grotm2[3];
-	IssmDouble Grotm3[3];
-	IssmDouble pre;
+	IssmDouble  Grotm1[3],GUrotm1[3],GNrotm1[3],GErotm1[3];
+	IssmDouble  Grotm2[3],GUrotm2[3],GNrotm2[3],GErotm2[3];
+	IssmDouble  Grotm3[3],GUrotm3[3],GNrotm3[3],GErotm3[3];
+	IssmDouble  Y21cos     , Y21sin     , Y20;
+	IssmDouble dY21cos_dlat,dY21sin_dlat,dY20_dlat;
+	IssmDouble dY21cos_dlon,dY21sin_dlon;
+	IssmDouble LoveRotRSL,LoveRotU,LoveRothoriz;
 	/*}}}*/
 	/*Recover parameters:{{{ */
@@ -6153,8 +6155,7 @@
 
 	if(computerotation){
-		parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
 		parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
 		parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
-		parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
+		parameters->FindParam(&tide_love_l,NULL,NULL,TidalLoveLEnum);
 		parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
 		parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
@@ -6270,16 +6271,61 @@
 		g=4.0/3.0*M_PI*rho_earth*NewtonG*planetradius;
 
+		//Amplitude of the rotational feedback
+		LoveRotRSL=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
+		LoveRotU=(tide_love_h[2]/g)*pow(omega*planetradius,2.0);
+		LoveRothoriz=(tide_love_l[2]/g)*pow(omega*planetradius,2.0);
+
 		for(int i=0;i<3;i++){
 			lati=latitude[i];
 			longi=longitude[i];
 
-			pre=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
-			Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
-			Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
-			Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati));
+			//Spherical harmonic functions of degree 2
+			Y21cos= -0.5*sin(2.*lati)*cos(longi);
+			Y21sin= -0.5*sin(2.*lati)*sin(longi);
+			Y20   = -(1.0/6.0 - 0.5*cos(2.0*lati));
+
+			Grotm1[i]=  LoveRotRSL*Y21cos;
+			Grotm2[i]=  LoveRotRSL*Y21sin;
+			Grotm3[i]=  LoveRotRSL*Y20;
+
+			if (computeelastic){
+				GUrotm1[i]=  LoveRotU*Y21cos;
+				GUrotm2[i]=  LoveRotU*Y21sin;
+				GUrotm3[i]=  LoveRotU*Y20;
+				if (horiz){
+				//bed_N = Love_l * d(Y)/dlat ;
+				dY21cos_dlat=-cos(2.*lati)*cos(longi);
+				dY21sin_dlat=-cos(2.*lati)*sin(longi);
+				dY20_dlat=-sin(2.0*lati);
+				GNrotm1[i]=  LoveRothoriz*dY21cos_dlat;
+				GNrotm2[i]=  LoveRothoriz*dY21sin_dlat;
+				GNrotm3[i]=  LoveRothoriz*dY20_dlat;
+
+				//bed_E = Love_l * 1/cos(lat) * d(Y)/dlon ;
+				dY21cos_dlon=-Y21sin;
+				dY21sin_dlon=Y21cos;
+				//dY20_dlon=0.;
+				GNrotm1[i]=  LoveRothoriz*dY21cos_dlon/cos(lati);
+				GNrotm2[i]=  LoveRothoriz*dY21sin_dlon/cos(lati);
+				GNrotm3[i]=  0.0;
+				}
+			}
 		}
 		this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
 		this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
 		this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
+		if (computeelastic){
+			this->AddInput(SealevelGUrotm1Enum,&GUrotm1[0],P1Enum);
+			this->AddInput(SealevelGUrotm2Enum,&GUrotm2[0],P1Enum);
+			this->AddInput(SealevelGUrotm3Enum,&GUrotm3[0],P1Enum);
+			if(horiz){
+				this->AddInput(SealevelGNrotm1Enum,&GNrotm1[0],P1Enum);
+				this->AddInput(SealevelGNrotm2Enum,&GNrotm2[0],P1Enum);
+				this->AddInput(SealevelGNrotm3Enum,&GNrotm3[0],P1Enum);
+				this->AddInput(SealevelGErotm1Enum,&GErotm1[0],P1Enum);
+				this->AddInput(SealevelGErotm2Enum,&GErotm2[0],P1Enum);
+				this->AddInput(SealevelGErotm3Enum,&GErotm3[0],P1Enum);
+			}
+		}
 	}
 	/*}}}*/
@@ -6995,4 +7041,13 @@
 	IssmDouble Grotm2[3];
 	IssmDouble Grotm3[3];
+	IssmDouble GUrotm1[3];
+	IssmDouble GUrotm2[3];
+	IssmDouble GUrotm3[3];
+	IssmDouble GNrotm1[3];
+	IssmDouble GNrotm2[3];
+	IssmDouble GNrotm3[3];
+	IssmDouble GErotm1[3];
+	IssmDouble GErotm2[3];
+	IssmDouble GErotm3[3];
 	bool rotation= false;
 	bool elastic=false;
@@ -7053,4 +7108,31 @@
 			}
 		}
+
+		if (elastic){
+			Element::GetInputListOnVertices(&GUrotm1[0],SealevelGUrotm1Enum);
+			Element::GetInputListOnVertices(&GUrotm2[0],SealevelGUrotm2Enum);
+			Element::GetInputListOnVertices(&GUrotm3[0],SealevelGUrotm3Enum);
+		
+			for(int i=0;i<NUMVERTICES;i++){
+				if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+					SealevelU[i]+=GUrotm1[i]*rotationvector[0]+GUrotm2[i]*rotationvector[1]+GUrotm3[i]*rotationvector[2];
+				}
+			}
+			if (horiz){
+				Element::GetInputListOnVertices(&GNrotm1[0],SealevelGNrotm1Enum);
+				Element::GetInputListOnVertices(&GNrotm2[0],SealevelGNrotm2Enum);
+				Element::GetInputListOnVertices(&GNrotm3[0],SealevelGNrotm3Enum);
+				Element::GetInputListOnVertices(&GErotm1[0],SealevelGErotm1Enum);
+				Element::GetInputListOnVertices(&GErotm2[0],SealevelGErotm2Enum);
+				Element::GetInputListOnVertices(&GErotm3[0],SealevelGErotm3Enum);
+		
+				for(int i=0;i<NUMVERTICES;i++){
+					if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+						SealevelN[i]+=GNrotm1[i]*rotationvector[0]+GNrotm2[i]*rotationvector[1]+GNrotm3[i]*rotationvector[2];
+						SealevelE[i]+=GErotm1[i]*rotationvector[0]+GErotm2[i]*rotationvector[1]+GErotm3[i]*rotationvector[2];
+					}
+				}
+			}
+		}
 	}
 
@@ -7062,6 +7144,6 @@
 	this->AddInput(BedGRDEnum,SealevelU,P1Enum);
 	if(horiz){
+		this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
 		this->AddInput(BedEastGRDEnum,SealevelE,P1Enum);
-		this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
 	}