Changeset 27232 for issm/trunk/src/c/classes/Loads/Friction.cpp

Timestamp:

08/25/22 16:50:29 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 27230

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c/classes/Loads/Friction.cpp (modified) (6 diffs)

issm/trunk/src/c/classes/Loads/Friction.cpp

@@ -321 +321 @@
 
         /*This routine calculates the basal friction coefficient
-          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
+          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p
+          alpha2= min(drag^2 * Neff ^r * | vel | ^(s-1), drag_coulomb^2 * Neff*/
 
         /*diverse: */
         IssmDouble  r,s;
         IssmDouble  drag_p, drag_q;
-        IssmDouble  thickness,base,floatation_thickness,sealevel;
         IssmDouble  drag_coefficient,drag_coefficient_coulomb;
         IssmDouble  alpha2,alpha2_coulomb;
@@ -333 +333 @@
         element->GetInputValue(&drag_p,gauss,FrictionPEnum);
         element->GetInputValue(&drag_q,gauss,FrictionQEnum);
-        element->GetInputValue(&thickness, gauss,ThicknessEnum);
-        element->GetInputValue(&base, gauss,BaseEnum);
-        element->GetInputValue(&sealevel, gauss,SealevelEnum);
         element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
         element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
@@ -347 +344 @@
 
         /*Get effective pressure*/
-        IssmDouble Neff = EffectivePressure(gauss);
+        bool ispwStochastic;
+        IssmDouble Neff;
+        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
+        if(ispwStochastic){
+                /*Retrieve stochastic water pressure and compute ice pressure*/
+                IssmDouble p_ice,p_water,Neff_limit;
+                element->GetInputValue(&p_water,gauss,FrictionCoulombWaterPressureEnum);
+                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
+                p_ice = IcePressure(gauss);
+                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
+        }       
+        else{
+                /*Compute effective pressure directly*/
+                Neff = EffectivePressure(gauss);
+        }
+        
+        /*Get velocity magnitude*/
         IssmDouble vmag = VelMag(gauss);
 
@@ -360 +373 @@
         }
 
-        floatation_thickness=0;
-        if(base<0) floatation_thickness=-(rho_water/rho_ice)*base;
         if(vmag==0.){
                 alpha2_coulomb=0.;
         }
         else{
-                alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
+                //alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
+                alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*Neff/vmag;
         }
 
@@ -749 +761 @@
          *               C |u_b|^(m-1)
          * alpha2= __________________________
-         *          (1+(C/(Cmax N))^1/m |u_b| )^m
+         *          (1+(C/(Cmax Neff))^1/m |u_b| )^m
          *
          * */
@@ -764 +776 @@
         C = coeff*coeff;
 
-        /*Get effective pressure and velocity magnitude*/
-        IssmDouble N  = EffectivePressure(gauss);
+        /*Get effective pressure*/
+        bool ispwStochastic;
+        IssmDouble Neff;
+        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
+        if(ispwStochastic){
+                /*Retrieve stochastic water pressure and compute ice pressure*/
+                IssmDouble p_ice,p_water,Neff_limit;
+                element->GetInputValue(&p_water,gauss,FrictionSchoofWaterPressureEnum);
+                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
+                p_ice = IcePressure(gauss);
+                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
+        }       
+        else{
+                /*Compute effective pressure directly*/
+                Neff = EffectivePressure(gauss);
+        }
+
+        /*Get velocity magnitude*/
         IssmDouble ub = VelMag(gauss);
 
         /*Compute alpha^2*/
-        if((ub<1e-10) ||(N==0.0)){
+        if((ub<1e-10) ||(Neff==0.0)){
                 alpha2 = 0.;
         }
         else{
-                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
+                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*Neff),1./m)*ub,m);
         }