Index: issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21295)
+++ issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21344)
@@ -2459,5 +2459,5 @@
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
 	ns = elements->Size();
-	
+
 	/*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
 	for(int i=0;i<ns;i++){
@@ -2465,4 +2465,5 @@
 		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());
@@ -2475,16 +2476,88 @@
 	for(int i=0;i<nsmax;i++){
 		if(i<ns){
-
 			if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
-
 			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,eartharea);
+			element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
 		}
 		if(i%100==0)pSgo->Assemble();
 	}
 	if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-
+	
 	/*Free ressources:*/
 	xDelete<IssmDouble>(Sg_old);
+}
+/*}}}*/
+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, 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]; 
+
+	for(int i=0;i<vertices->Size();i++){
+		//Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
+		Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
+		
+		lati=latitude[i]/180*PI;	longi=longitude[i]/180*PI;
+
+		/*only first order terms are considered now: */ 
+		value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radius[i],2.0)*
+						(-m3/6.0 - m3*cos(2.0*lati)/2.0 - sin(lati)*cos(lati)*(m1*cos(longi)+m2*sin(longi))); 
+	
+		pSgo_rot->SetValue(vertex->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);
+	
 }
 /*}}}*/
