Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.cpp
===================================================================
--- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.cpp	(revision 21394)
+++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.cpp	(revision 21395)
@@ -2402,4 +2402,85 @@
 }
 /*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+
+	/*serialized vectors:*/
+	IssmDouble* Sg_old=NULL;
+	IssmDouble  eartharea=0;
+	IssmDouble  eartharea_cpu=0;
+	IssmDouble	tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
+	IssmDouble	load_love_k2 = -0.30922675; //degree 2 load Love number 
+	IssmDouble	m1, m2, m3; 
+	IssmDouble	lati, longi, radi, value; 
+
+	/*Serialize vectors from previous iteration:*/
+	Sg_old=pSg_old->ToMPISerial();
+
+	/*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
+	for(int i=0;i<elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+		eartharea_cpu += element->GetAreaSpherical();
+	}
+	ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+	IssmDouble moi_list[3]={0,0,0}; 
+	IssmDouble moi_list_cpu[3]={0,0,0}; 
+	for(int i=0;i<elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+		element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
+		moi_list_cpu[0] += moi_list[0]; 
+		moi_list_cpu[1] += moi_list[1]; 
+		moi_list_cpu[2] += moi_list[2]; 
+	}
+	ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	// 	
+	ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	// 	
+	ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	
+	/*pull out some useful parameters: */
+	parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
+	parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
+	parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum);
+	parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum);
+	parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum);
+	parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
+
+	/*compute perturbation terms for angular velocity vector: */
+	m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; 
+	m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; 
+	m3 = -(1+load_love_k2)/moi_p * moi_list[2];	// term associated with fluid number (3-order-of-magnitude smaller) is negelected  
+
+	/* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
+	 * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) 
+	 * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) 
+	 * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] 
+	 */
+	for(int i=0;i<vertices->Size();i++){
+		int sid;
+		//Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
+		Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
+		sid=vertex->Sid();
+
+		lati=latitude[sid]/180*PI;	longi=longitude[sid]/180*PI; radi=radius[sid];
+
+		/*only first order terms are considered now: */ 
+		value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)*
+						(-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); 
+	
+		pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
+	}
+
+	/*Assemble mesh velocity*/
+	pSgo_rot->Assemble();
+	
+	/*Free ressources:*/
+	xDelete<IssmDouble>(Sg_old);
+	
+}
+/*}}}*/
 void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
 
Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.h
===================================================================
--- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.h	(revision 21394)
+++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.h	(revision 21395)
@@ -116,4 +116,5 @@
 		void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
 		void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
+		void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
 		void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
 		IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);