source: issm/oecreview/Archive/21724-22754/ISSM-21793-21794.diff@ 22755

../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -469 +469 @@
 
         /*Intermediaries */
         IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+        IssmDouble  ms,mb,gmb,fmb,thickness;
+        IssmDouble  gllevelset,phi=1.;
         IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -484 +485 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
-        Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(groundedice_input);
+        Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
         Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
 
@@ -499 +500 @@
                 ms_input->GetInputValue(&ms,gauss);
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                groundedice_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&phi,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
                 if(phi>0.) mb=gmb;
                 else mb=fmb;
@@ -519 +520 @@
         if(!element->IsIceInElement()) return NULL;
 
         /*Intermediaries */
-        IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+        int         migration_style,point1;
+        bool        mainlyfloating;
+        IssmDouble  Jdet,fraction1,fraction2;
+        IssmDouble  dt,gllevelset;
+        IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
         IssmDouble* xyz_list = NULL;
+        Gauss*      gauss     = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = element->GetNumberOfNodes();
@@ -536 +541 @@
         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);          _assert_(ms_input);
-        Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(groundedice_input);
+        Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
         Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
 
+        /*Recover portion of element that is grounded*/
+        phi=element->GetGroundedPortion(xyz_list);
+        if(migration_style==SubelementMigration2Enum){
+                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
+        }
+        else{
+                gauss = element->NewGauss(3);
+        }
+
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
 
@@ -550 +565 @@
                 ms_input->GetInputValue(&ms,gauss);
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                groundedice_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&gllevelset,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
-                if(phi>0) mb=gmb;
-                else mb=fmb;
+                if(migration_style==SubelementMigrationEnum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+                }
+                else if(migration_style==SubelementMigration2Enum){
+                        if(gllevelset>0.) mb=gmb;
+                        else mb=fmb;
+                }
+                else if(migration_style==SubelementMigration3Enum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=gmb;
+                }
+                else{
+                        if(gllevelset>0) mb=gmb;
+                        else mb=fmb;
+                }
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
         }