Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21349)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21350)
@@ -3930,7 +3930,12 @@
 		S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
 		
+		/* Perturbation terms for moment of inertia (moi_list): 
+		 * computed analytically (see Wu & Peltier, eqs 10 & 32) 
+		 * also consistent with my GMD formulation!
+		 * ALL in geographic coordinates 
+		 * */
 		dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea; 
 		dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea; 
-		dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(cos(late),2))/eartharea; 
+		dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; 
 	}
 	else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
@@ -3947,5 +3952,5 @@
 		dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea; 
 		dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea; 
-		dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(cos(late),2))/eartharea; 
+		dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; 
 	}
 	
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21349)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21350)
@@ -2539,6 +2539,11 @@
 	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]; 
-
+	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++){
 		//Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
@@ -2549,5 +2554,5 @@
 		/*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))); 
+						(-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*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