Index: /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp	(revision 24102)
+++ /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp	(revision 24103)
@@ -2,23 +2,28 @@
  * \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; 
 }
Index: /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp	(revision 24102)
+++ /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp	(revision 24103)
@@ -2,17 +2,8 @@
  * 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,6 +18,9 @@
 
 	/*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*/