source: issm/oecreview/Archive/24307-24683/ISSM-24550-24551.diff

../trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -378 +378 @@
          */
 
         /*Intermediaries: */
-        IssmDouble  T,Tpmp,deltaT,deltaTref,pressure;
-        IssmDouble  alpha2,time,gamma;
+        IssmDouble  T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient;
+        IssmDouble  alpha2,time,gamma,ref,alp_new,alphascaled;
         const IssmDouble yts = 365*24*3600.;
 
         /*Get viscous part*/
@@ -387 +387 @@
 
         /*Get delta Refs*/
         element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
+        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
+        /*New*/
+        /*element->GetInputValue(&deltaTrefsfc,gauss,FrictionSurfaceTemperatureEnum);
+         *    element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum);
+         *       */
 
         /*Compute delta T*/
         element->GetInputValue(&T,gauss,TemperatureEnum);
@@ -394 +399 @@
         Tpmp = element->TMeltingPoint(pressure);
         deltaT = T-Tpmp;
 
+
         /*Compute gamma*/
         element->parameters->FindParam(&time,TimeEnum);
         element->parameters->FindParam(&gamma,FrictionGammaEnum);
-        //if(time<25e3*yts){
-        //      gamma = 10.;
-        //}
-        //else{
-        //      gamma = 5.;
-        //}
-        //gamma = 5.;
 
-        /*Compute scaling parameter*/
-        alpha2 = alpha2 * exp((deltaTref - deltaT)/(2*gamma));
+        ref = exp(deltaTref/gamma);
+        alp_new = ref/exp(deltaT/gamma);
 
+        alphascaled = sqrt(alp_new)*drag_coefficient;
+        if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient);
+
+        alp_new=alp_new*alpha2;
+
         /*Assign output pointers:*/
-        *palpha2=alpha2;
+        *palpha2=alp_new;
 }/*}}}*/
 void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/