Index: /issm/trunk-jpl/src/m/materials/nye.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/nye.m	(revision 24100)
+++ /issm/trunk-jpl/src/m/materials/nye.m	(revision 24101)
@@ -3,34 +3,41 @@
 %
 %   Compute 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
+%   rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law
+%   sigma=B*e(1/n) (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) && (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)');
+
+	% Beyond-melting-point cases
+	warning OFF BACKTRACE
+	if (ice_type==1)
+		if (any(temperature>200&temperature<220))
+			warning('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.\nLook at indexes: %s', mat2str(find(temperature>200 & temperature<220))');
+		end
+		if (any(temperature>=220))
+			warning('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.\nLook at indexes: %s', mat2str(find(temperature>=220))');
+		end
+	elseif ((ice_type==2)&&(any(temperature>273.15)))
+		warning('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.\nLook at indexes: %s', mat2str(find(temperature>273.15))');
 	end
 
-	Rg = 8.3144598; % J mol^-1 K^-1
+	% Coefficients
+	Rg=8.3144598;       % J mol^-1 K^-1
 
-	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==1)     % CO2 ice
+		A_const     = 10^(10.8);    % s^-1 MPa
+		Q           = 63000;        % J mol^-1
+		n           = 7;            % Glen's exponent
 	elseif(ice_type==2) % H2O ice
-		A_const = 9e4;   % s^-1 MPa
-		Q       = 60000; % J mol^-1
-		n       = 3;     % Glen's exponent
-
+		A_const     = 9e4;          % s^-1 MPa
+		Q           = 60000;        % J mol^-1
+		n           = 3;            % Glen's exponent
 	else
-		error('ice type not supported');
+		error('Ice type not supported');
 	end
 
-	%Arhenius law
-	A = A_const*exp(-Q./(temperature*Rg)); % s^-1 MPa
-	rigidity = A.^(-1/n)*1e6; % s^(1/n) Pa
+	% Arrhenius Law
+	A=A_const*exp(-Q./(temperature*Rg));  % s^-1 MPa
+	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 24100)
+++ /issm/trunk-jpl/src/m/materials/nye.py	(revision 24101)
@@ -4,5 +4,5 @@
 	"""
    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).
+	rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law sigma=B*e(1/n) (Nye, p2000).
 	temperature is in Kelvin degrees
 
@@ -11,32 +11,41 @@
 	"""
 
-	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)")
+        # Declaring temperature and rigidity arrays
+        if np.ndim(temperature)==2:
+            T=temperature.flatten()
+        elif isinstance(temperature,float) or isinstance(temperature,int):
+            T=np.array([temperature])
+        else:
+            T=temperature
+        rigidity=np.zeros_like(T)
+
+        # Beyond-melting-point cases
+        if (ice_type==1):
+            for i in range(len(T)):
+                if (200<T[i]<220):
+                    warnings.warn('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.')
+                break
+            if ((T>=220).any()):
+                warnings.warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.')
+        elif (ice_type==2) and ((T>273.15).any()):
+            warnings.warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.')
 
 	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
+	elif ice_type == 2: # H2O ice
+	    A_const = 9*10**4 # s^-1 MPa
+	    Q = 60000. #  J mol^-1
+	    n = 3. # Glen's exponent
+        else:
+            raise RuntimeError('Ice type not supported')
 
-	if ice_type == 1: # CO2 ice
-		A_const = 10**(10.8) # s^-1 MPa
-		Q = 63000. # J mol^-1
-		n = 7. # Glen's exponent
+        # Arrhenius Law
+        A=A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa
+        rigidity=A**(-1/n)*10**6 # s^(1/n) Pa
 
-	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 
+        # Return output
+        return rigidity