Changeset 28115 for issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

Timestamp:

02/29/24 15:22:20 (13 months ago)

Author:

poinelli

Message:

NEW: Add IntrusionMelt module to include melt under grounding zones, edits after Robel, Wilson, Seroussi TC (2022)

File:

• issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (modified) (11 diffs)

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

@@ -156 +156 @@
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
+        
         if(isgroundingline)     iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
         /*Initialize ThicknessResidual input*/
@@ -272 +273 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.min_thickness",MasstransportMinThicknessEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.penalty_factor",MasstransportPenaltyFactorEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.intrusion_distance",GroundinglineIntrusionDistanceEnum));
 
         iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.masstransport.requested_outputs");
@@ -622 +624 @@
         bool        mainlyfloating;
         IssmDouble  fraction1,fraction2;
-        IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,fmb_pert;
+        IssmDouble  Jdet,dt,intrusiondist;
+        IssmDouble  ms,mb,gmb,fmb,thickness,fmb_pert,gldistance;
         IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau;
         IssmDouble  dvx[2],dvy[2];
@@ -652 +654 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         element->FindParam(&stabilization,MasstransportStabilizationEnum);
+        element->FindParam(&intrusiondist,GroundinglineIntrusionDistanceEnum);
+
         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
@@ -660 +664 @@
         Input* vxaverage_input  = element->GetInput(VxAverageEnum);                                                                             _assert_(vxaverage_input);
         Input* vyaverage_input  = element->GetInput(VyAverageEnum);                                                                             _assert_(vyaverage_input);
-
+        //Input* gldistance_input = element->GetInput(DistanceToGroundinglineEnum);              _assert_(gldistance_input); 
         h=element->CharacteristicLength();
 
@@ -667 +671 @@
         if(melt_style==SubelementMelt2Enum){
                 element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
-           gauss = element->NewGauss(point1,fraction1,fraction2,3);
+            gauss = element->NewGauss(point1,fraction1,fraction2,3);
+        }
+        else if(melt_style==IntrusionMeltEnum){
+            element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+       gauss = element->NewGauss(point1,fraction1,fraction2,3);
         }
         else{
@@ -704 +712 @@
                         else mb=gmb;
                 }
+        else if(melt_style==IntrusionMeltEnum){
+                        Input* gldistance_input = element->GetInput(DistanceToGroundinglineEnum);              _assert_(gldistance_input); 
+                gldistance_input->GetInputValue(&gldistance,gauss);
+                        if (intrusiondist==0)
+                                if(gllevelset>0.) mb=gmb;
+                                else mb=fmb;
+                else if(gldistance>intrusiondist) 
+                                mb=gmb;
+                        else if(gldistance<=intrusiondist && gldistance>0) 
+                                mb=fmb*(1-gldistance/intrusiondist); 
+                        else
+                                mb=fmb;
+        }
                 else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
 
@@ -750 +771 @@
         bool mainlyfloating;
         IssmDouble  fraction1,fraction2,gllevelset;
-        IssmDouble  Jdet,dt;
-        IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
+        IssmDouble  Jdet,dt,intrusiondist;
+        IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.,gldistance;
         IssmDouble* xyz_list = NULL;
 
@@ -765 +786 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
+        element->FindParam(&intrusiondist,GroundinglineIntrusionDistanceEnum);
+
         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
@@ -778 +801 @@
       gauss = element->NewGauss(point1,fraction1,fraction2,3);
    }
+   else if(melt_style==IntrusionMeltEnum){
+            element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+       gauss = element->NewGauss(point1,fraction1,fraction2,3);
+        }
    else{
       gauss = element->NewGauss(3);
@@ -798 +825 @@
          else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
       }
-      else if(melt_style==SubelementMelt2Enum){
+        else if(melt_style==SubelementMelt2Enum){
          if(gllevelset>0.) mb=gmb;
          else mb=fmb;
       }
-      else if(melt_style==NoMeltOnPartiallyFloatingEnum){
+        else if(melt_style==NoMeltOnPartiallyFloatingEnum){
          if (phi<0.00000001) mb=fmb;
          else mb=gmb;
       }
-      else if(melt_style==FullMeltOnPartiallyFloatingEnum){
+        else if(melt_style==FullMeltOnPartiallyFloatingEnum){
          if (phi<0.99999999) mb=fmb;
          else mb=gmb;
-      }
-      else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
+        }
+                else if(melt_style==IntrusionMeltEnum){
+                        Input* gldistance_input = element->GetInput(DistanceToGroundinglineEnum);              _assert_(gldistance_input); 
+                gldistance_input->GetInputValue(&gldistance,gauss);
+                        if (intrusiondist==0)
+                                if(gllevelset>0.) mb=gmb;
+                                else mb=fmb;
+                else if(gldistance>intrusiondist) 
+                                mb=gmb;
+                        else if(gldistance<=intrusiondist && gldistance>0) 
+                                mb=fmb*(1-gldistance/intrusiondist); 
+                        else
+                                mb=fmb;
+        }
+        else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];