source: issm/branches/trunk-larour-NatGeoScience2016/src/m/materials/paterson.py@ 21759

Last change on this file since 21759 was 21759, checked in by Eric.Larour, 8 years ago

CHG: merged branch back to trunk-jpl 21754.
File size: 3.5 KB
Line 
1import numpy as np
2
3def paterson(temperature):
4 """
5 PATERSON - figure out the rigidity of ice for a given temperature
6
7 rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).
8 temperature is in Kelvin degrees
9
10 Usage:
11 rigidity=paterson(temperature)
12 """
13
14 if np.any(temperature<0.):
15 raise RuntimeError("input temperature should be in Kelvin (positive)")
16
17 if np.ndim(temperature)==2:
18 #T = temperature.reshape(-1,)-273.15
19 T = temperature.flatten()-273.15
20 elif isinstance(temperature,float) or isinstance(temperature,int):
21 T = np.array([temperature])-273.15
22 else:
23 T = temperature-273.15
24
25 #The routine below is equivalent to:
26
27 # n=3; T=temperature-273;
28 # %From paterson,
29 # Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50];
30 # A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4*
31 # 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3)
32 # %Convert into rigidity B
33 # B=A.^(-1/n)*10^3; %s^(1/3)Pa
34 # %Now, do a cubic fit between Temp and B:
35 # fittedmodel=fit(Temp,B,'cubicspline');
36 # rigidity=fittedmodel(temperature);
37
38 rigidity=np.zeros_like(T)
39 pos1=np.nonzero(T<=-45)[0]
40 if len(pos1):
41 rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101)
42 pos2=np.nonzero(np.logical_and(-45<=T,T<-40))[0]
43 if len(pos2):
44 rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554)
45 pos3=np.nonzero(np.logical_and(-40<=T,T<-35))[0]
46 if len(pos3):
47 rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214)
48 pos4=np.nonzero(np.logical_and(-35<=T,T<-30))[0]
49 if len(pos4):
50 rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331)
51 pos5=np.nonzero(np.logical_and(-30<=T,T<-25))[0]
52 if len(pos5):
53 rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191)
54 pos6=np.nonzero(np.logical_and(-25<=T,T<-20))[0]
55 if len(pos6):
56 rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181)
57 pos7=np.nonzero(np.logical_and(-20<=T,T<-15))[0]
58 if len(pos7):
59 rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111)
60 pos8=np.nonzero(np.logical_and(-15<=T,T<-10))[0]
61 if len(pos8):
62 rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523)
63 pos9=np.nonzero(np.logical_and(-10<=T,T<-5))[0]
64 if len(pos9):
65 rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714)
66 pos10=np.nonzero(np.logical_and(-5<=T,T<-2))[0]
67 if len(pos10):
68 rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348)
69 pos11=np.nonzero(-2<=T)[0]
70 if len(pos11):
71 rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860)
72
73 #Now make sure that rigidity is positive
74 pos=np.nonzero(rigidity<0)[0]
75 if len(pos):
76 rigidity[pos]=1.e6
77
78 return rigidity
79
Note: See TracBrowser for help on using the repository browser.