| 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 locals () : |
|---|
| 73 | if locals()[char] == None : unDefined = True |
|---|
| 74 | else : unDefined = False |
|---|
| 75 | else : unDefined = True |
|---|
| 76 | return unDefined |
|---|
| 77 | |
|---|
| 78 | def Defined (char) : |
|---|
| 79 | '''True if a variable is defined and not equal to None''' |
|---|
| 80 | if char in globals () : |
|---|
| 81 | if globals()[char] == None : Defined = False |
|---|
| 82 | else : Defined = True |
|---|
| 83 | else : Defined = False |
|---|
| 84 | return Defined |
|---|
| 85 | |
|---|
| 86 | def Ksum (tab) : |
|---|
| 87 | ''' |
|---|
| 88 | Kahan summation algorithm, also known as compensated summation |
|---|
| 89 | https://en.wikipedia.org/wiki/Kahan_summation_algorithm |
|---|
| 90 | ''' |
|---|
| 91 | Ksum = 0.0 # Prepare the accumulator. |
|---|
| 92 | comp = 0.0 # A running compensation for lost low-order bits. |
|---|
| 93 | |
|---|
| 94 | for xx in np.where ( np.isnan(tab), 0., tab ) : |
|---|
| 95 | yy = xx - comp # comp is zero the first time around. |
|---|
| 96 | tt = Ksum + yy # Alas, sum is big, y small, so low-order digits of y are lost. |
|---|
| 97 | comp = (tt - Ksum) - yy # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy) |
|---|
| 98 | Ksum = tt # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers ! |
|---|
| 99 | # Next time around, the lost low part will be added to y in a fresh attempt. |
|---|
| 100 | return Ksum |
|---|
| 101 | |
|---|
| 102 | def Ssum (tab) : |
|---|
| 103 | ''' |
|---|
| 104 | Precision summation by sorting by absolute values |
|---|
| 105 | ''' |
|---|
| 106 | Ssum = np.sum ( tab[np.argsort(np.abs(tab))] ) |
|---|
| 107 | return Ssum |
|---|
| 108 | |
|---|
| 109 | def KSsum (tab) : |
|---|
| 110 | ''' |
|---|
| 111 | Precision summation by sorting by absolute value, and applying Kahan summation |
|---|
| 112 | ''' |
|---|
| 113 | KSsum = Ksum ( tab[np.argsort(np.abs(tab))] ) |
|---|
| 114 | return KSsum |
|---|
| 115 | |
|---|
| 116 | # Choosing algorithm for precise summation |
|---|
| 117 | Psum = KSsum |
|---|
| 118 | |
|---|
| 119 | def IsLeapYear ( Year, CalendarType="Gregorian" ) : |
|---|
| 120 | ''' True is Year is a leap year ''' |
|---|
| 121 | |
|---|
| 122 | # What is the calendar : |
|---|
| 123 | if CalendarType in [ '360d', '360_day', 'noleap', '365_day'] : |
|---|
| 124 | IsLeapYear = False |
|---|
| 125 | return IsLeapYear |
|---|
| 126 | |
|---|
| 127 | if CalendarType in [ 'all_leap', '366_day' ] : |
|---|
| 128 | IsLeapYear = True |
|---|
| 129 | return IsLeapYear |
|---|
| 130 | |
|---|
| 131 | Year = int ( Year ) |
|---|
| 132 | |
|---|
| 133 | # a year is a leap year if it is even divisible by 4 |
|---|
| 134 | # but not evenly divisible by 100 |
|---|
| 135 | # unless it is evenly divisible by 400 |
|---|
| 136 | |
|---|
| 137 | # if it is evenly divisible by 400 it must be a leap year |
|---|
| 138 | if np.mod ( Year, 400 ) == 0 : |
|---|
| 139 | IsLeapYear = True |
|---|
| 140 | return IsLeapYear |
|---|
| 141 | |
|---|
| 142 | # if it is evenly divisible by 100 it must not be a leap year |
|---|
| 143 | if np.mod ( Year, 100 ) == 0 : |
|---|
| 144 | IsLeapYear = False |
|---|
| 145 | return IsLeapYear |
|---|
| 146 | |
|---|
| 147 | # if it is evenly divisible by 4 it must be a leap year |
|---|
| 148 | if np.mod ( Year, 4 ) == 0 : |
|---|
| 149 | IsLeapYear = True |
|---|
| 150 | return IsLeapYear |
|---|
| 151 | |
|---|
| 152 | IsLeapYear = False |
|---|
| 153 | return IsLeapYear |
|---|
| 154 | |
|---|
| 155 | def DateFormat ( Date ) : |
|---|
| 156 | ''' |
|---|
| 157 | Get date format : |
|---|
| 158 | [yy]yymmdd is Gregorian |
|---|
| 159 | [yy]yy-mm-dd is Human |
|---|
| 160 | ''' |
|---|
| 161 | if type (Date) == str : |
|---|
| 162 | if '-' in Date : |
|---|
| 163 | DateFormat = 'Human' |
|---|
| 164 | else : |
|---|
| 165 | DateFormat = 'Gregorian' |
|---|
| 166 | if type (Date) == int : DateFormat = 'Gregorian' |
|---|
| 167 | return DateFormat |
|---|
| 168 | |
|---|
| 169 | def PrintDate ( YE, MO, DA, Format ) : |
|---|
| 170 | '''Return a date in the requested format''' |
|---|
| 171 | if Format == 'Human' : PrintDate = f'{YE:04d}-{MO:02d}-{DA:02d}' |
|---|
| 172 | if Format == 'Gregorian' : PrintDate = f'{YE:04d}{MO:02d}{DA:02d}' |
|---|
| 173 | return PrintDate |
|---|
| 174 | |
|---|
| 175 | def FormatToGregorian ( Date ) : |
|---|
| 176 | ''' from a yyyy-mm-dd or yyymmdd date format returns |
|---|
| 177 | a yyymmdd date format |
|---|
| 178 | ''' |
|---|
| 179 | YE, MO, DA = SplitDate ( Date ) |
|---|
| 180 | FormatToGregorian = f'{YE:04d}{MO:02d}{DA:02d}' |
|---|
| 181 | return FormatToGregorian |
|---|
| 182 | |
|---|
| 183 | def FormatToHuman ( Date ) : |
|---|
| 184 | ''' From a yyyymmdd or yyymmdd date format returns |
|---|
| 185 | a yyy-mm-dd date format |
|---|
| 186 | ''' |
|---|
| 187 | YE, MO, DA = SplitDate ( Date ) |
|---|
| 188 | FormatToHuman = f'{YE_new:04d}-{MO:02d}-{DA:02d}' |
|---|
| 189 | return FormatToHuman |
|---|
| 190 | |
|---|
| 191 | def SplitDate ( Date, Out='int' ) : |
|---|
| 192 | ''' Split Date in format [yy]yymmdd or [yy]yy-mm-dd |
|---|
| 193 | to yy, mm, dd ''' |
|---|
| 194 | if type (Date) == str : |
|---|
| 195 | if '-' in Date : |
|---|
| 196 | Dash = True |
|---|
| 197 | YE, MO, DA = Date.split ('-') |
|---|
| 198 | else : |
|---|
| 199 | Dash = False |
|---|
| 200 | YE = Date[:-4] ; MO = Date[-4:-2] ; DA = Date[-2:] |
|---|
| 201 | YearDigits = len (YE) |
|---|
| 202 | if type (Date) == int : |
|---|
| 203 | DA = np.mod ( Date, 100) |
|---|
| 204 | MO = np.mod ( Date//100, 100) |
|---|
| 205 | YE = Date // 10000 |
|---|
| 206 | |
|---|
| 207 | YE = int(YE) ; MO = int(MO) ; DA=int(DA) |
|---|
| 208 | return YE, MO, DA |
|---|
| 209 | |
|---|
| 210 | def DateAddYear ( Date, YearInc=1 ) : |
|---|
| 211 | ''' Add on year to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
|---|
| 212 | Format = DateFormat ( Date ) |
|---|
| 213 | YE, MO, DA = SplitDate ( Date ) |
|---|
| 214 | YE_new = YE + YearInc |
|---|
| 215 | DateAddYear = PrintDate ( YE_new, MO, DA, Format) |
|---|
| 216 | return DateAddYear |
|---|
| 217 | |
|---|
| 218 | def DateMinusOneDay ( Date ) : |
|---|
| 219 | ''' Substracts one day to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
|---|
| 220 | Format = DateFormat ( Date ) |
|---|
| 221 | mth_length = np.array ( [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] ) |
|---|
| 222 | YE, MO, DA = SplitDate ( Date ) |
|---|
| 223 | if IsLeapYear ( YE) : mth_length[1] += 1 |
|---|
| 224 | |
|---|
| 225 | YE = int(YE) ; MO = int(MO) ; DA=int(DA) |
|---|
| 226 | if DA == 1 : |
|---|
| 227 | if MO == 1 : DA_new = mth_length[-1] ; MO_new = 12 ; YE_new = YE - 1 |
|---|
| 228 | else : DA_new = mth_length[MO-2] ; MO_new = MO - 1 ; YE_new = YE |
|---|
| 229 | else : DA_new = DA - 1 ; MO_new = MO ; YE_new = YE |
|---|
| 230 | |
|---|
| 231 | DateMinusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format) |
|---|
| 232 | return DateMinusOneDay |
|---|
| 233 | |
|---|
| 234 | def DatePlusOneDay ( Date ) : |
|---|
| 235 | ''' Add one day to date in format [yy]yymmdd or [yy]yy-mm-dd''' |
|---|
| 236 | Format = DateFormat ( Date ) |
|---|
| 237 | mth_length = np.array ( [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] ) |
|---|
| 238 | YE, MO, DA = SplitDate ( Date ) |
|---|
| 239 | if IsLeapYear ( YE ) : mth_length[1] += 1 |
|---|
| 240 | |
|---|
| 241 | YE_new = YE ; MO_new=MO ; DA_new=DA+1 |
|---|
| 242 | if DA_new > mth_length [MO_new-1] : |
|---|
| 243 | DA_new = 1 ; MO_new = MO_new + 1 |
|---|
| 244 | if MO_new == 13 : |
|---|
| 245 | MO_new = 1 ; YE_new = YE_new+1 |
|---|
| 246 | |
|---|
| 247 | DatePlusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format ) |
|---|
| 248 | return DatePlusOneDay |
|---|
