Changeset 21375

Timestamp:

11/14/16 13:12:25 (8 years ago)

Author:

Eric.Larour

Message:

CHG: better handling of fraction of grounded area that contributes SLR.

File:

• issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp (modified) (4 diffs)

issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp

@@ -1066 +1066 @@
 
         }
-        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum){
+        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
                 /*Check that not all nodes are grounded or floating*/
                 if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
@@ -3622 +3622 @@
         IssmDouble late,longe,re;
         IssmDouble lati,longi,ri;
+        bool notfullygrounded=false;
 
         /*elastic green function:*/
@@ -3642 +3643 @@
                 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
                 return;
+        }
+
+        /*early return if we are fully floating: */
+        if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0){
+                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
+                return;
+        }
+        
+        /*If we are here, we are on ice that is fully grounded or half-way to floating: */
+        if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
+                notfullygrounded=true; //used later on.
         }
 
@@ -3703 +3715 @@
         /*}}}*/
 
-        /*Compute area of element:*/
+        /*Compute area of element. Scale it by grounded fraction if not fully grounded: */
         area=GetAreaSpherical();
-
+        if(notfullygrounded){
+                IssmDouble phi=0;
+                IssmDouble xyz_list[NUMVERTICES][3];
+                ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+                phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
+                area*=phi;
+        }
+        
         /*Compute ice thickness change: */
         Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); 
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
-        deltathickness_input->GetInputAverage(&I);
-
-        /*Are we on an ice shelf? if so, does not contribute to sea level rise!: */
-        if (this->IsFloating()) I=0;
+
+        /*If we are fully grounded, take the average over the element: */
+        if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+        else{ 
+                IssmDouble total_weight=0;
+                bool mainlyfloating = true;
+                int         point1;
+                IssmDouble  fraction1,fraction2;
+                
+                /*Recover portion of element that is grounded*/
+                this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+                Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+
+                /* Start  looping on the number of gaussian points and average over these gaussian points: */
+                total_weight=0;
+                I=0;
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+                        IssmDouble Ig=0;
+                        gauss->GaussPoint(ig);
+                        deltathickness_input->GetInputValue(&Ig);
+                        I+=Ig*gauss->weight;
+                        total_weight+=gauss->weight;
+                }
+                I=I/total_weight;
+                delete gauss;
+        }
 
         /*Compute eustatic compoent:*/