Index: /issm/trunk-jpl/src/m/materials/nye.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/nye.m	(revision 24076)
+++ /issm/trunk-jpl/src/m/materials/nye.m	(revision 24077)
@@ -8,7 +8,8 @@
 %   Usage:
 %      rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
-
-	if (any(temperature < 130) || any(temperature > 273))
-		error('Input temperature out of bounds (130<T<273)');
+	if ((ice_type == 2) && (any(temperature > 195)))
+		error('Input temperature for CO2 ice out of bounds (T>195)');
+	elseif ((ice_type == 2) && (any(temperature > 273.15)))
+		error('Input temperature for H2O ice out of bounds (T>273.15)');
 	end
 
@@ -31,5 +32,5 @@
 	%Arhenius law
 	A = A_const*exp(-Q./(temperature*Rg)); % s^-1 MPa
-	B = A.^(-1/n)*1e6; % s^(1/n) Pa
+	rigidity = A.^(-1/n)*1e6; % s^(1/n) Pa
 
 end
Index: /issm/trunk-jpl/src/m/materials/nye.py
===================================================================
--- /issm/trunk-jpl/src/m/materials/nye.py	(revision 24077)
+++ /issm/trunk-jpl/src/m/materials/nye.py	(revision 24077)
@@ -0,0 +1,42 @@
+import numpy as np
+
+def nye(temperature, ice_type):
+	"""
+   NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
+	rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) (Nye, p2000).
+	temperature is in Kelvin degrees
+
+   Usage:
+	   rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
+	"""
+
+	if (ice_type == 1) and (np.any(temperature > 195)):
+		raise RuntimeError("Input temperature for CO2 ice out of bounds (T>195)")
+	elif (ice_type == 2) and (np.any(temperature > 273)):
+		raise RuntimeError("Input temperature for H2O ice out of bounds (T>273)")
+
+	Rg = 8.3144598 # J mol^-1 K^-1
+
+	if np.ndim(temperature)==2:
+		T = temperature.flatten()
+	elif isinstance(temperature,float) or isinstance(temperature,int):
+		T = np.array([temperature])
+	else:
+		T = temperature
+		
+	B=np.zeros_like(T)
+
+	if ice_type == 1: # CO2 ice
+		A_const = 10**(10.8) # s^-1 MPa
+		Q = 63000 # J mol^-1
+		n = 7 # Glen's exponent
+
+	if ice_type == 2: # H2O ice
+		A_const = 9*10**4 # s^-1 MPa
+		Q = 60000 #  J mol^-1
+		n = 3 # Glen's exponent
+
+	A = A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa
+	rigidity = A**(-1/n)*10**6 # s^(1/n) Pa
+
+	return rigidity