1 | #!/usr/bin/env python3 |
---|
2 | ''' |
---|
3 | Script to check water conservation in the IPSL coupled model |
---|
4 | |
---|
5 | Here, some common utilities to all scripts |
---|
6 | |
---|
7 | Warning, to install, configure, run, use any of included software or |
---|
8 | to read the associated documentation you'll need at least one (1) |
---|
9 | brain in a reasonably working order. Lack of this implement will |
---|
10 | void any warranties (either express or implied). Authors assumes |
---|
11 | no responsability for errors, omissions, data loss, or any other |
---|
12 | consequences caused directly or indirectly by the usage of his |
---|
13 | software by incorrectly or partially configured personal |
---|
14 | |
---|
15 | |
---|
16 | SVN information |
---|
17 | $Author$ |
---|
18 | $Date$ |
---|
19 | $Revision$ |
---|
20 | $Id: ATM_waterbudget.py 6277 2022-12-08 09:24:05Z omamce $ |
---|
21 | $HeadURL$ |
---|
22 | ''' |
---|
23 | |
---|
24 | import numpy as np |
---|
25 | import configparser, re |
---|
26 | |
---|
27 | VarInt = 10 |
---|
28 | Rho = 1000.0 |
---|
29 | Ra = 6366197.7236758135 #-- Earth Radius (m) |
---|
30 | Grav = 9.81 #-- Gravity (m^2/s |
---|
31 | ICE_rho_ice = 917.0 #-- Ice volumic mass (kg/m3) in LIM3 |
---|
32 | ICE_rho_sno = 330.0 #-- Snow volumic mass (kg/m3) in LIM3 |
---|
33 | OCE_rho_liq = 1026.0 #-- Ocean water volumic mass (kg/m3) in NEMO |
---|
34 | ATM_rho = 1000.0 #-- Water volumic mass in atmosphere (kg/m^3) |
---|
35 | SRF_rho = 1000.0 #-- Water volumic mass in surface reservoir (kg/m^3) |
---|
36 | RUN_rho = 1000.0 #-- Water volumic mass of rivers (kg/m^3) |
---|
37 | ICE_rho_pnd = 1000. #-- Water volumic mass in ice ponds (kg/m^3) |
---|
38 | YearLength = 365.25 * 24. * 60. * 60. #-- Year length (s) |
---|
39 | |
---|
40 | def setBool (chars) : |
---|
41 | '''Convert specific char string in boolean if possible''' |
---|
42 | setBool = chars |
---|
43 | if type (chars) == str : |
---|
44 | for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : |
---|
45 | if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] |
---|
46 | return setBool |
---|
47 | |
---|
48 | def setNum (chars) : |
---|
49 | '''Convert specific char string in integer or real if possible''' |
---|
50 | if type (chars) == str : |
---|
51 | realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") |
---|
52 | isNumber = realnum.match(chars.strip()) != None |
---|
53 | isInt = chars.strip().isdigit() |
---|
54 | #print (chars , isNumber , isInt) |
---|
55 | if isNumber : |
---|
56 | if isInt : setNum = int (chars) |
---|
57 | else : setNum = float (chars) |
---|
58 | else : setNum = chars |
---|
59 | else : setNum = chars |
---|
60 | return setNum |
---|
61 | |
---|
62 | def setNone (chars) : |
---|
63 | '''Convert specific char string to None if possible''' |
---|
64 | if type (chars) == str : |
---|
65 | if chars in ['None', 'NONE', 'none'] : setNone = None |
---|
66 | else : setNone = chars |
---|
67 | else : setNone = chars |
---|
68 | return setNone |
---|
69 | |
---|
70 | def unDefined (char) : |
---|
71 | '''True if a variable is not defined, or set to None''' |
---|
72 | if char in globals () : |
---|
73 | if globals()[char] == None : |
---|
74 | unDefined = True |
---|
75 | else : unDefined = False |
---|
76 | else : unDefined = True |
---|
77 | return unDefined |
---|
78 | |
---|
79 | def Defined (char) : |
---|
80 | '''True if a variable is defined and not equal to None''' |
---|
81 | if char in globals () : |
---|
82 | if globals()[char] == None : |
---|
83 | Defined = False |
---|
84 | else : Defined = True |
---|
85 | else : Defined = False |
---|
86 | return Defined |
---|
87 | |
---|
88 | def Ksum (tab) : |
---|
89 | ''' |
---|
90 | Kahan summation algorithm, also known as compensated summation |
---|
91 | https://en.wikipedia.org/wiki/Kahan_summation_algorithm |
---|
92 | ''' |
---|
93 | Ksum = 0.0 # Prepare the accumulator. |
---|
94 | comp = 0.0 # A running compensation for lost low-order bits. |
---|
95 | |
---|
96 | for xx in np.where ( np.isnan(tab), 0., tab ) : |
---|
97 | yy = xx - comp # comp is zero the first time around. |
---|
98 | tt = Ksum + yy # Alas, sum is big, y small, so low-order digits of y are lost. |
---|
99 | comp = (tt - Ksum) - yy # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy) |
---|
100 | Ksum = tt # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers ! |
---|
101 | # Next time around, the lost low part will be added to y in a fresh attempt. |
---|
102 | return Ksum |
---|
103 | |
---|
104 | def Ssum (tab) : |
---|
105 | ''' |
---|
106 | Precision summation by sorting by absolute values |
---|
107 | ''' |
---|
108 | Ssum = np.sum ( tab[np.argsort(np.abs(tab))] ) |
---|
109 | return Ssum |
---|
110 | |
---|
111 | def KSsum (tab) : |
---|
112 | ''' |
---|
113 | Precision summation by sorting by absolute value, and applying Kahan summation |
---|
114 | ''' |
---|
115 | KSsum = Ksum ( tab[np.argsort(np.abs(tab))] ) |
---|
116 | return KSsum |
---|
117 | |
---|
118 | # Choosing algorithm for precise summation |
---|
119 | Psum = KSsum |
---|
120 | |
---|
121 | def IsLeapYear ( Year, CalendarType="Gregorian" ) : |
---|
122 | ''' True is Year is a leap year ''' |
---|
123 | |
---|
124 | # What is the calendar : |
---|
125 | if CalendarType in [ '360d', '360_day', 'noleap', '365_day'] : |
---|
126 | IsLeapYear = False |
---|
127 | return IsLeapYear |
---|
128 | |
---|
129 | if CalendarType in [ 'all_leap', '366_day' ] : |
---|
130 | IsLeapYear = True |
---|
131 | return IsLeapYear |
---|
132 | |
---|
133 | Year = int ( Year ) |
---|
134 | |
---|
135 | # a year is a leap year if it is even divisible by 4 |
---|
136 | # but not evenly divisible by 100 |
---|
137 | # unless it is evenly divisible by 400 |
---|
138 | |
---|
139 | # if it is evenly divisible by 400 it must be a leap year |
---|
140 | if np.mod ( Year, 400 ) == 0 : |
---|
141 | IsLeapYear = True |
---|
142 | return IsLeapYear |
---|
143 | |
---|
144 | # if it is evenly divisible by 100 it must not be a leap year |
---|
145 | if np.mod ( Year, 100 ) == 0 : |
---|
146 | IsLeapYear = False |
---|
147 | return IsLeapYear |
---|
148 | |
---|
149 | # if it is evenly divisible by 4 it must be a leap year |
---|
150 | if np.mod ( Year, 4 ) == 0 : |
---|
151 | IsLeapYear = True |
---|
152 | return IsLeapYear |
---|
153 | |
---|
154 | IsLeapYear = False |
---|
155 | return IsLeapYear |
---|
156 | |
---|
157 | def DateFormat ( Date ) : |
---|
158 | ''' |
---|
159 | Get date format : |
---|
160 | [yy]yymmdd is Gregorian |
---|
161 | [yy]yy-mm-dd is Human |
---|
162 | ''' |
---|
163 | if type (Date) == str : |
---|
164 | if '-' in Date : |
---|
165 | DateFormat = 'Human' |
---|
166 | else : |
---|
167 | DateFormat = 'Gregorian' |
---|
168 | if type (Date) == int : DateFormat = 'Gregorian' |
---|
169 | return DateFormat |
---|
170 | |
---|
171 | def PrintDate ( YE, MO, DA, Format ) : |
---|
172 | '''Return a date in the requested format''' |
---|
173 | if Format == 'Human' : PrintDate = f'{YE:04d}-{MO:02d}-{DA:02d}' |
---|
174 | if Format == 'Gregorian' : PrintDate = f'{YE:04d}{MO:02d}{DA:02d}' |
---|
175 | return PrintDate |
---|
176 | |
---|
177 | def FormatToGregorian ( Date ) : |
---|
178 | ''' from a yyyy-mm-dd or yyymmdd date format returns |
---|
179 | a yyymmdd date format |
---|
180 | ''' |
---|
181 | YE, MO, DA = SplitDate ( Date ) |
---|
182 | FormatToGregorian = f'{YE:04d}{MO:02d}{DA:02d}' |
---|
183 | return FormatToGregorian |
---|
184 | |
---|
185 | def FormatToHuman ( Date ) : |
---|
186 | ''' From a yyyymmdd or yyymmdd date format returns |
---|
187 | a yyy-mm-dd date format |
---|
188 | ''' |
---|
189 | YE, MO, DA = SplitDate ( Date ) |
---|
190 | FormatToHuman = f'{YE_new:04d}-{MO:02d}-{DA:02d}' |
---|
191 | return FormatToHuman |
---|
192 | |
---|
193 | def SplitDate ( Date, Out='int' ) : |
---|
194 | ''' Split Date in format [yy]yymmdd or [yy]yy-mm-dd |
---|
195 | to yy, mm, dd ''' |
---|
196 | if type (Date) == str : |
---|
197 | if '-' in Date : |
---|
198 | Dash = True |
---|
199 | YE, MO, DA = Date.split ('-') |
---|
200 | else : |
---|
201 | Dash = False |
---|
202 | YE = Date[:-4] ; MO = Date[-4:-2] ; DA = Date[-2:] |
---|
203 | YearDigits = len (YE) |
---|
204 | if type (Date) == int : |
---|
205 | DA = np.mod ( Date, 100) |
---|
206 | MO = np.mod ( Date//100, 100) |
---|
207 | YE = Date // 10000 |
---|
208 | |
---|
209 | YE = int(YE) ; MO = int(MO) ; DA=int(DA) |
---|
210 | return YE, MO, DA |
---|
211 | |
---|
212 | def DateAddYear ( Date, YearInc=1 ) : |
---|
213 | ''' Add on year to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
---|
214 | Format = DateFormat ( Date ) |
---|
215 | YE, MO, DA = SplitDate ( Date ) |
---|
216 | YE_new = YE + YearInc |
---|
217 | DateAddYear = PrintDate ( YE_new, MO, DA, Format) |
---|
218 | return DateAddYear |
---|
219 | |
---|
220 | def DateMinusOneDay ( Date ) : |
---|
221 | ''' Substracts one day to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
---|
222 | Format = DateFormat ( Date ) |
---|
223 | mth_length = np.array ( [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] ) |
---|
224 | YE, MO, DA = SplitDate ( Date ) |
---|
225 | if IsLeapYear ( YE) : mth_length[1] += 1 |
---|
226 | |
---|
227 | YE = int(YE) ; MO = int(MO) ; DA=int(DA) |
---|
228 | if DA == 1 : |
---|
229 | if MO == 1 : DA_new = mth_length[-1] ; MO_new = 12 ; YE_new = YE - 1 |
---|
230 | else : DA_new = mth_length[MO-2] ; MO_new = MO - 1 ; YE_new = YE |
---|
231 | else : DA_new = DA - 1 ; MO_new = MO ; YE_new = YE |
---|
232 | |
---|
233 | DateMinusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format) |
---|
234 | return DateMinusOneDay |
---|
235 | |
---|
236 | def DatePlusOneDay ( Date ) : |
---|
237 | ''' Add one day to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
---|
238 | Format = DateFormat ( Date ) |
---|
239 | mth_length = np.array ( [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] ) |
---|
240 | YE, MO, DA = SplitDate ( Date ) |
---|
241 | if IsLeapYear ( YE ) : mth_length[1] += 1 |
---|
242 | |
---|
243 | YE_new = YE ; MO_new=MO ; DA_new=DA+1 |
---|
244 | if DA_new > mth_length [MO_new-1] : |
---|
245 | DA_new = 1 ; MO_new = MO_new + 1 |
---|
246 | if MO_new == 13 : |
---|
247 | MO_new = 1 ; YE_new = YE_new+1 |
---|
248 | |
---|
249 | DatePlusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format ) |
---|
250 | return DatePlusOneDay |
---|