[6508] | 1 | #!/usr/bin/env python3 |
---|
[6620] | 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 |
---|
[6508] | 14 | |
---|
[6620] | 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 | |
---|
[6508] | 24 | import numpy as np |
---|
| 25 | import configparser, re |
---|
| 26 | |
---|
[6620] | 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) |
---|
[6508] | 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 : |
---|
[6620] | 51 | realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") |
---|
[6508] | 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) : |
---|
[6620] | 71 | '''True if a variable is not defined, or set to None''' |
---|
[6647] | 72 | if char in globals () : |
---|
| 73 | if globals()[char] == None : |
---|
| 74 | unDefined = True |
---|
[6508] | 75 | else : unDefined = False |
---|
| 76 | else : unDefined = True |
---|
| 77 | return unDefined |
---|
| 78 | |
---|
[6620] | 79 | def Defined (char) : |
---|
| 80 | '''True if a variable is defined and not equal to None''' |
---|
| 81 | if char in globals () : |
---|
[6647] | 82 | if globals()[char] == None : |
---|
| 83 | Defined = False |
---|
[6620] | 84 | else : Defined = True |
---|
| 85 | else : Defined = False |
---|
| 86 | return Defined |
---|
| 87 | |
---|
[6508] | 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 | ''' |
---|
[6620] | 106 | Precision summation by sorting by absolute values |
---|
[6508] | 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 | |
---|
[6620] | 121 | def IsLeapYear ( Year, CalendarType="Gregorian" ) : |
---|
| 122 | ''' True is Year is a leap year ''' |
---|
[6508] | 123 | |
---|
[6620] | 124 | # What is the calendar : |
---|
| 125 | if CalendarType in [ '360d', '360_day', 'noleap', '365_day'] : |
---|
| 126 | IsLeapYear = False |
---|
| 127 | return IsLeapYear |
---|
[6508] | 128 | |
---|
[6620] | 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 |
---|
[6508] | 138 | |
---|
[6620] | 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 |
---|