source: issm/oecreview/Archive/24684-25833/ISSM-24793-24794.diff@ 25834

../trunk-jpl/src/m/materials/nye.py

@@ -1 +1 @@
 import numpy as np
+from warnings import warn
 
+
 def nye(temperature, ice_type):
-        """
-   NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
-        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
+    """
+    NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
+        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
-        """
+           rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
+        """
 
-        # 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)
+    # 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.')
+    # Beyond-melting-point cases
+    if (ice_type == 1):
+        for i in range(len(T)):
+            if (200 < T[i] < 220):
+                warn('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.')
+            break
+        if ((T >= 220).any()):
+            warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.')
+    elif (ice_type == 2) and ((T > 273.15).any()):
+        warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.')
 
-        Rg = 8.3144598              # J mol^-1 K^-1
+    Rg = 8.3144598              # J mol^-1 K^-1
 
-        if ice_type == 1:           # CO2 ice
-            A_const = 10**(13.0)    # s^-1 MPa
-            Q = 66900.              # J mol^-1
-            n = 8.                  # 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 = 1.0e13    # s^-1 MPa
+        Q = 66900.              # J mol^-1
+        n = 8.                  # Glen's exponent
+    elif ice_type == 2:         # H2O ice
+        A_const = 9 * 1.0e4       # s^-1 MPa
+        Q = 60000.              #  J mol^-1
+        n = 3.                  # Glen's exponent
+    else:
+        raise RuntimeError('Ice type not supported')
 
-        # Arrhenius Law
-        A=A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa
-        rigidity=A**(-1/n)*10**6 # s^(1/n) Pa
+    # Arrhenius Law
+    A = A_const * np.exp(-1 * Q / (T * Rg))  # s^-1 MPa
+    rigidity = A**(-1 / n) * 1.0e6  # s^(1/n) Pa
 
-        # Return output
-        return rigidity
+    # Return output
+    return rigidity