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

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 3.9 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/m/materials/nye.py
2===================================================================
3--- ../trunk-jpl/src/m/materials/nye.py (revision 24793)
4+++ ../trunk-jpl/src/m/materials/nye.py (revision 24794)
5@@ -1,51 +1,53 @@
6 import numpy as np
7+from warnings import warn
8
9+
10 def nye(temperature, ice_type):
11- """
12- NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
13- rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law sigma=B*e(1/n) (Nye, p2000).
14- temperature is in Kelvin degrees
15+ """
16+ NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
17+ rigidity (in s^(1/n)Pa) is the flow law parameter in the flow law sigma=B*e(1/n) (Nye, p2000).
18+ temperature is in Kelvin degrees
19
20 Usage:
21- rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
22- """
23+ rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
24+ """
25
26- # Declaring temperature and rigidity arrays
27- if np.ndim(temperature)==2:
28- T=temperature.flatten()
29- elif isinstance(temperature,float) or isinstance(temperature,int):
30- T=np.array([temperature])
31- else:
32- T=temperature
33- rigidity=np.zeros_like(T)
34+ # Declaring temperature and rigidity arrays
35+ if np.ndim(temperature) == 2:
36+ T = temperature.flatten()
37+ elif isinstance(temperature, float) or isinstance(temperature, int):
38+ T = np.array([temperature])
39+ else:
40+ T = temperature
41+ rigidity = np.zeros_like(T)
42
43- # Beyond-melting-point cases
44- if (ice_type==1):
45- for i in range(len(T)):
46- if (200<T[i]<220):
47- warnings.warn('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.')
48- break
49- if ((T>=220).any()):
50- warnings.warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.')
51- elif (ice_type==2) and ((T>273.15).any()):
52- warnings.warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.')
53+ # Beyond-melting-point cases
54+ if (ice_type == 1):
55+ for i in range(len(T)):
56+ if (200 < T[i] < 220):
57+ warn('CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.')
58+ break
59+ if ((T >= 220).any()):
60+ warn('CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.')
61+ elif (ice_type == 2) and ((T > 273.15).any()):
62+ warn('H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.')
63
64- Rg = 8.3144598 # J mol^-1 K^-1
65+ Rg = 8.3144598 # J mol^-1 K^-1
66
67- if ice_type == 1: # CO2 ice
68- A_const = 10**(13.0) # s^-1 MPa
69- Q = 66900. # J mol^-1
70- n = 8. # Glen's exponent
71- elif ice_type == 2: # H2O ice
72- A_const = 9*10**4 # s^-1 MPa
73- Q = 60000. # J mol^-1
74- n = 3. # Glen's exponent
75- else:
76- raise RuntimeError('Ice type not supported')
77+ if ice_type == 1: # CO2 ice
78+ A_const = 1.0e13 # s^-1 MPa
79+ Q = 66900. # J mol^-1
80+ n = 8. # Glen's exponent
81+ elif ice_type == 2: # H2O ice
82+ A_const = 9 * 1.0e4 # s^-1 MPa
83+ Q = 60000. # J mol^-1
84+ n = 3. # Glen's exponent
85+ else:
86+ raise RuntimeError('Ice type not supported')
87
88- # Arrhenius Law
89- A=A_const*np.exp(-1*Q/(T*Rg)) # s^-1 MPa
90- rigidity=A**(-1/n)*10**6 # s^(1/n) Pa
91+ # Arrhenius Law
92+ A = A_const * np.exp(-1 * Q / (T * Rg)) # s^-1 MPa
93+ rigidity = A**(-1 / n) * 1.0e6 # s^(1/n) Pa
94
95- # Return output
96- return rigidity
97+ # Return output
98+ return rigidity
Note: See TracBrowser for help on using the repository browser.