Changeset 24103

Timestamp:

07/22/19 14:15:42 (6 years ago)

Author:

schlegel

Message:

CHG: update some warnings for the Nye functions

Location:

issm/trunk-jpl/src/c/shared/Elements

Files:

• NyeCO2.cpp (modified) (1 diff)
• NyeH2O.cpp (modified) (2 diffs)

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

@@ -2 +2 @@
  * \brief figure out B of CO2 ice for a certain temperature
  *              INPUT function B=NyeCO2(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- ). 
+ *      where rigidigty (in s^(1/n)Pa) is the flow law paramter in the flow law sigma=B*e(1/n) (Nye, p2000). 
  */
 
-#include <math.h>
-#include "../Numerics/types.h"
+#include "../io/io.h" #include <math.h> #include "../Numerics/types.h"
 
 IssmDouble NyeCO2(IssmDouble temperature){
 
         /*Coefficients*/
-        const IssmPDouble Rg      = 8.3144598;     /* J mol^-1 K^-1 */
-        const IssmPDouble A_const = pow(10.,10.8); /*s^-1 MPa       */
-        const IssmPDouble Q       = 63000.;        /*J mol^-1       */
-        const IssmPDouble n       = 7.;            /*Glen's exponent*/
+        const IssmPDouble Rg      = 8.3144598;     /* J mol^-1 K^-1 */ const
+          IssmPDouble A_const = pow(10.,10.8); /*s^-1 MPa       */ const IssmPDouble
+          Q       = 63000.;        /*J mol^-1       */ const IssmPDouble n       =
+          7.;            /*Glen's exponent*/
 
         /*Arrhenius Law*/
-        IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa   */
+        IssmDouble A = A_const *exp(-Q/(temperature*Rg));  /*s^-1 MPa   */
         IssmDouble B = 1e6*pow(A,-1/n);                    /*s^(1/n) Pa */
 
+        /*Beyond-melting-point cases*/
+        if((temperature>200.)&&(temperature<220.)) printf("CO2 ICE - POSSIBLE
+                                MELTING. Some temperature values are between 200K and 220K.\n");
+        else if(temperature>=220.) _printf0_("CO2 ICE - GUARANTEED MELTING. Some
+                                temperature values are beyond 220K.\n");
+
         /*Return output*/
-        return B;
+        return B; 
 }

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

@@ -2 +2 @@
  * 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- ). 
-
- VALUES CALCULATED USING THE FOLLOWING MATLAB CODE
-
- Rg = 8.3144598; % J mol^-1 K^-1
- A_const = 9e4; % s^-1 MPa
- Q = 60000; % J mol^-1
- n = 3; % Glen's exponent
- T = -102.5:5:2.5; T = T + 273; % K
- A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
- B = A.^(-1/n)*1e6; % s^(1/n) Pa
+ *      where rigidigty (in s^(1/n)Pa) is the flow law paramter in the flow law sigma=B*e(1/n) (Nye, p2000). 
  */
 
+#include "../io/io.h"
 #include <math.h>
 #include "../Numerics/types.h"
@@ -27 +18 @@
 
         /*Arrhenius Law*/
-        IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa   */
+        IssmDouble A = A_const *exp(-Q/(temperature*Rg));  /*s^-1 MPa   */
         IssmDouble B = 1e6*pow(A,-1/n);                    /*s^(1/n) Pa */
+
+        /*Beyond-melting-point case*/
+        if(temperature>=273.15) _printf0_("H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.\n");
 
         /*Return output*/