Changeset 24101

Timestamp:

07/22/19 14:05:11 (6 years ago)

Author:

schlegel

Message:

CHG: minor updates to nye functions

Location:

issm/trunk-jpl/src/m/materials

Files:

• nye.m (modified) (1 diff)
• nye.py (modified) (2 diffs)

issm/trunk-jpl/src/m/materials/nye.m ¶

@@ -3 +3 @@
 %
 %   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

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

@@ -4 +4 @@
         """
    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 +11 @@
         """
 
-        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