Changeset 24073 for issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp

Timestamp:

07/09/19 08:28:35 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: improved poorly coded law

File:

• issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp

@@ -1 +1 @@
 /* \file NyeH2O.cpp
- * \    brief figure out B of H2O ice for a certain temperature
+ * brief figure out B of H2O ice for a certain temperature
  *              INPUT function B=NyeH2O(temperature)
  *      where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Nye, p -COMPLETE- ). 
@@ -20 +20 @@
 IssmDouble NyeH2O(IssmDouble temperature){
 
-        IssmDouble B,T;
+        /*Coefficients*/
+        const IssmPDouble Rg      = 8.3144598; /* J mol^-1 K^-1 */
+        const IssmPDouble A_const = 9.e4;      /*s^-1 MPa       */
+        const IssmPDouble Q       = 60000.;    /*J mol^-1       */
+        const IssmPDouble n       = 3.;        /*Glen's exponent*/
 
-        /*Switch to celsius from Kelvin: */
-        T=temperature-273.15;
+        /*Arrhenius Law*/
+        IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa   */
+        IssmDouble B = 1e6*pow(A,-1/n);                    /*s^(1/n) Pa */
 
-        if(T<=-100){
-                B=29901899089.3257;
-        }
-        else if((T>=-100.0) && (T<=-95.0)){
-                B=20004935214.6693;
-        }
-        else if((T>=-95.0) && (T<=-90.0)){
-                B=13685054490.7813;
-        }
-        else if((T>=-90.0) && (T<=-85.0)){
-                B=9555312498.71763;
-        }
-        else if((T>=-85.0) && (T<=-80.0)){
-                B=6798800864.99471;
-        }
-        else if((T>=-80.0) && (T<=-75.0)){
-                B=4922440695.83365;
-        }
-        else if((T>=-75.0) && (T<=-70.0)){
-                B=3621794191.06010;
-        }
-        else if((T>=-70.0) && (T<=-65.0)){
-                B=2704902125.98878;
-        }
-        else if((T>=-65.0) && (T<=-60.0)){
-                B=2048338724.47737;
-        }
-        else if((T>=-60.0) && (T<=-55.0)){
-                B=1571285974.28101;
-        }
-        else if((T>=-55.0) && (T<=-50.0)){
-                B=1219918332.04930;
-        }
-        else if((T>=-50.0) && (T<=-45.0)){
-                B=957813559.547873;
-        }
-        else if((T>=-45.0) && (T<=-40.0)){
-                B=759956402.731125;
-        }
-        else if((T>=-40.0) && (T<=-35.0)){
-                B=608924663.484756;
-        }
-        else if((T>=-35.0) && (T<=-30.0)){
-                B=492424321.822968;
-        }
-        else if((T>=-30.0) && (T<=-25.0)){
-                B=401672414.952195;
-        }
-        else if((T>=-25.0) && (T<=-20.0)){
-                B=330320977.865296;
-        }
-        else if((T>=-20.0) && (T<=-15.0)){
-                B=273731390.267536;
-        }
-        else if((T>=-15.0) && (T<=-10.0)){
-                B=228478811.983490;
-        }
-        else if((T>=-10.0) && (T<=-5.0)){
-                B=192009688.386368;
-        }
-        else if((T>=-5.0) && (T<=0.0)){
-                B=162402355.392225;
-        }
-        else{
-                B=138197905.740112;
-        }
-
+        /*Return output*/
         return B;
 }