Changeset 21350

Timestamp:

11/07/16 17:22:11 (8 years ago)

Author:

adhikari

Message:

CHG: constants appearing in SLR rotational feedback double checked

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Tria.cpp (modified) (2 diffs)
• FemModel.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3930 +3930 @@
                 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 +3952 @@
                 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; 
         }
         

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -2539 +2539 @@
         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 +2554 @@
                 /*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