Changeset 6277 for TOOLS/WATER_BUDGET
- Timestamp:
- 12/08/22 10:24:05 (19 months ago)
- Location:
- TOOLS/WATER_BUDGET
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/WATER_BUDGET/ATM_waterbudget.py
r6271 r6277 21 21 ### 22 22 ## Import system modules 23 import sys, os, shutil, subprocess 23 import sys, os, shutil, subprocess, platform 24 24 import numpy as np 25 25 import configparser, re … … 29 29 config.optionxform = str # To keep capitals 30 30 31 config['Files'] = {}32 config['System'] = {}33 34 ##-- Some physical constants35 #-- Earth Radius36 Ra = 6366197.723675813537 #-- Gravity38 Grav = 9.8139 #-- Ice volumic mass (kg/m3) in LIM340 ICE_rho_ice = 917.041 #-- Snow volumic mass (kg/m3) in LIM342 ICE_rho_sno = 330.043 #-- Ocean water volumic mass (kg/m3) in NEMO44 OCE_rho_liq = 1026.45 #-- Water volumic mass in atmosphere46 ATM_rho = 1.0e347 #-- Water volumic mass in surface reservoirs48 SRF_rho = 1.0e349 #-- Water volumic mass of rivers50 RUN_rho = 1.0e351 52 31 ## Read experiment parameters 53 ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 32 ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ; TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 33 34 ## 35 ARCHIVE =None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 36 TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 37 file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ; file_ICE_his=None ; file_OCE_sca=None 38 file_restart_beg=None ; file_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None 39 file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_beg=None ; file_OCE_end=None 40 41 d_ATM_his=None ; d_SRF_his=None ; d_RUN_his=None ; d_OCE_his=None ; d_ICE_his=None ; d_OCE_sca=None 42 d_restart_beg=None ; d_restart_end=None ; d_ATM_beg=None ; d_ATM_end=None ; d_DYN_beg=None ; d_DYN_end=None ; d_SRF_beg=None ; d_SRF_end=None 43 d_RUN_beg=None ; d_RUN_end=None ; d_RUN_end=None ; d_OCE_beg=None ; d_ICE_beg=None ; d_OCE_beg=None ; d_OCE_end=None 54 44 55 45 # Arguments passed 56 46 print ( "Name of Python script:", sys.argv[0] ) 57 IniFile = 47 IniFile = sys.argv[1] 58 48 print ("Input file : ", IniFile ) 59 49 config.read (IniFile) 50 FullIniFile = 'full_' + IniFile 60 51 61 52 def setBool (chars) : … … 79 70 return setNum 80 71 81 print ('[Experiment]') 82 for VarName in config['Experiment'].keys() : 83 locals()[VarName] = config['Experiment'][VarName] 84 locals()[VarName] = setBool (locals()[VarName]) 85 locals()[VarName] = setNum (locals()[VarName]) 86 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 87 88 # ### 72 def setNone (chars) : 73 '''Convert specific char string to None if possible''' 74 if type (chars) == str : 75 if chars in ['None', 'NONE', 'none'] : setNone = None 76 else : setNone = chars 77 else : setNone = chars 78 return setNone 79 80 ## Reading config 81 for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 82 if Section in config.keys() : 83 print ( f'[{Section}]' ) 84 for VarName in config[Section].keys() : 85 locals()[VarName] = config[Section][VarName] 86 locals()[VarName] = setBool (locals()[VarName]) 87 locals()[VarName] = setNum (locals()[VarName]) 88 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 89 90 if not 'Files' in config.keys() : config['Files'] = {} 91 92 def unDefined (char) : 93 if char in globals () : 94 if char == None : return True 95 else : return False 96 else : return True 97 98 ##-- Some physical constants 99 #-- Earth Radius 100 if not 'Ra' in locals () : Ra = 6366197.7236758135 101 #-- Gravity 102 if not 'Grav' in locals () : Grav = 9.81 103 #-- Ice volumic mass (kg/m3) in LIM3 104 if not 'ICE_rho_ice' in locals () : ICE_rho_ice = 917.0 105 #-- Snow volumic mass (kg/m3) in LIM3 106 if not 'ICE_rho_sno' in locals () : ICE_rho_sno = 330.0 107 #-- Ocean water volumic mass (kg/m3) in NEMO 108 if not 'OCE_rho_liq' in locals () : OCE_rho_liq = 1026. 109 #-- Water volumic mass in atmosphere 110 if not 'ATM_rho' in locals () : ATM_rho = 1000. 111 #-- Water volumic mass in surface reservoirs 112 if not 'SRF_rho' in locals () : SRF_rho = 1000. 113 #-- Water volumic mass of rivers 114 if not 'RUN_rho' in locals () : RUN_rho = 1000. 115 #-- Year length 116 if not 'YearLength' in locals () : YearLength = 365.25 * 24. * 60. * 60. 117 118 config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 119 120 ## -- 89 121 ICO = ( 'ICO' in ATM ) 90 122 LMDZ = ( 'LMD' in ATM ) 91 123 92 ###93 ## Import system modules94 import sys, os, shutil, subprocess95 import configparser, re96 97 config = configparser.ConfigParser()98 config['Files'] = {}99 124 100 125 # Where do we run ? 101 SysName, NodeName, Release, Version, Machine = os.uname ()126 SysName, NodeName, Release, Version, Machine = os.uname () 102 127 TGCC = ( 'irene' in NodeName ) 103 128 IDRIS = ( 'jeanzay' in NodeName ) … … 108 133 if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene' 109 134 if "AMD" in CPU : Machine = 'irene-amd' 110 111 ARCHIVE = subprocess.getoutput ( f'ccc_home --cccstore -d {Project} -u {User}' ) 112 STORAGE = subprocess.getoutput ( f'ccc_home --cccwork -d {Project} -u {User}' ) 113 SCRATCHDIR = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 114 R_IN = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 115 rebuild = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 135 136 if libIGCM : 137 if ARCHIVE == None : ARCHIVE = subprocess.getoutput ( f'ccc_home --cccstore -d {Project} -u {User}' ) 138 if STORAGE == None : STORAGE = subprocess.getoutput ( f'ccc_home --cccwork -d {Project} -u {User}' ) 139 if SCRATCHDIR == None : SCRATCHDIR = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 140 if R_IN == None : R_IN = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 141 if rebuild == None : rebuild = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 116 142 117 143 ## Specific to run at TGCC. … … 121 147 122 148 ## Creates output directory 123 #TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 124 TmpDir = os.path.join ( '/ccc/scratch/cont003/drf/p86mart', f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 125 149 if TmpDir == None : 150 TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 151 config['Files']['TmpDir'] = TmpDir 152 126 153 if IDRIS : 127 154 raise Exception ("Pour IDRIS : repertoires et chemins a definir") 155 156 config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU, 157 'Program' : "Generated by : " + str(sys.argv), 158 'HOSTNAME' : platform.node (), 'LOGNAME' : os.getlogin (), 159 'Python' : f'{platform.python_version ()}', 160 'OS' : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}', 161 'SVN_Author' : "$Author$", 162 'SVN_Date' : "$Date$", 163 'SVN_Revision' : "$Revision$", 164 'SVN_Id' : "$Id$", 165 'SVN_HeadURL' : "$HeadURL$"} 166 167 if libIGCM : 168 config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 128 169 129 170 ## Import specific module … … 132 173 import xarray as xr 133 174 134 config['Files'][TmpDir] = TmpDir135 136 175 # Output file 137 FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 176 if FileOut == None : 177 FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 178 config['Files']['FileOut'] = FileOut 179 138 180 f_out = open ( FileOut, mode = 'w' ) 139 181 140 182 # Function to print to stdout *and* output file 141 183 def echo (string, end='\n') : 142 print ( str ing, end=end )184 print ( str(string), end=end ) 143 185 sys.stdout.flush () 144 f_out.write ( str ing+ end )186 f_out.write ( str(string) + end ) 145 187 f_out.flush () 146 188 return None … … 165 207 if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 166 208 209 echo (' ') 210 echo ( f'JobName : {JobName}' ) 211 echo (Comment) 167 212 echo ( f'Working in TMPDIR : {TmpDir}' ) 168 213 … … 172 217 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 173 218 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 174 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 175 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 176 219 if dir_ATM_his == None : 220 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 221 config['Files']['dir_ATM_his'] = dir_ATM_his 222 if dir_SRF_his == None : 223 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 224 config['Files']['dir_SRF_his'] = dir_SRF_his 225 177 226 echo ( f'The analysis relies on files from the following model output directories : ' ) 178 227 echo ( f'{dir_ATM_his}' ) … … 180 229 181 230 #-- Files Names 182 if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 183 if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 184 231 if Period == None : 232 if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 233 if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 234 config['Files']['Period'] = Period 235 if FileCommon == None : 236 FileCommon = f'{JobName}_{Period}' 237 config['Files']['FileCommon'] = FileCommon 238 239 if Title == None : 240 Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 241 config['Files']['Title'] = Title 242 185 243 echo ('\nOpen history files' ) 186 file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 187 file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 188 if Routing == 'ORCHIDEE' : 189 file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 244 if file_ATM_his == None : 245 file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 246 config['Files']['file_ATM_his'] = file_ATM_his 247 if file_SRF_his == None : 248 file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 249 config['Files']['file_SRF_his'] = file_SRF_his 250 #if Routing == 'SECHIBA' : 251 # file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 190 252 if Routing == 'SIMPLE' : 191 file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 253 if file_RUN_his == None : 254 file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 255 config['Files']['file_RUN_his'] = file_RUN_his 192 256 193 257 d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 194 258 d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 195 if Routing == ' ORCHIDEE' :259 if Routing == 'SECHIBA' : 196 260 d_RUN_his = d_SRF_his 197 261 if Routing == 'SIMPLE' : … … 200 264 echo ( file_ATM_his ) 201 265 echo ( file_SRF_his ) 202 echo ( file_RUN_his ) 203 204 205 config['Files']['file_ATM_his'] = file_ATM_his 206 config['Files']['file_SRF_his'] = file_SRF_his 207 config['Files']['file_RUN_his'] = file_SRF_his 266 if Routing == 'SIMPLE' : 267 echo ( file_RUN_his ) 208 268 209 269 ## Compute run length … … 218 278 dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 219 279 280 ## Number of years 281 NbYear = dtime_sec / YearLength 282 220 283 #-- Open restart files 284 221 285 YearRes = YearBegin - 1 # Year of the restart of beginning of simulation 222 286 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 223 287 224 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 225 226 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' ) 227 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' ) 228 288 if TarRestartPeriod_beg == None : 289 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 290 TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 291 config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 292 293 if TarRestartPeriod_end == None : 294 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 295 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 296 TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231' 297 config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 298 299 if file_restart_beg == None : 300 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 301 config['Files']['file_restart_beg'] = file_restart_beg 302 if file_restart_end == None : 303 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 304 config['Files']['file_restart_end'] = file_restart_end 305 229 306 echo ( f'{file_restart_beg}' ) 230 307 echo ( f'{file_restart_end}' ) 231 308 232 file_ATM_beg = f'ATM_{JobName}_{YearRes}1231_restartphy.nc' 233 file_ATM_end = f'ATM_{JobName}_{YearEnd}1231_restartphy.nc' 234 if LMDZ : 235 file_DYN_beg = f'ATM_{JobName}_{YearRes}1231_restart.nc' 236 file_DYN_end = f'ATM_{JobName}_{YearEnd}1231_restart.nc' 237 if ICO : 238 file_DYN_beg = f'ICO_{JobName}_{YearRes}1231_restart.nc' 239 file_DYN_end = f'ICO_{JobName}_{YearEnd}1231_restart.nc' 240 file_SRF_beg = f'SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 241 file_SRF_end = f'SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 242 liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg, ] 243 liste_end = [file_ATM_end, file_DYN_end, file_SRF_end, ] 309 if file_ATM_beg == None : 310 file_ATM_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc' 311 config['Files']['file_ATM_beg'] = file_ATM_beg 312 if file_ATM_end == None : 313 file_ATM_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc' 314 config['Files']['file_ATM_end'] = file_ATM_end 315 316 liste_beg = [file_ATM_beg, ] 317 liste_end = [file_ATM_end, ] 318 319 320 if file_DYN_beg == None : 321 if LMDZ : 322 file_DYN_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restart.nc' 323 if ICO : 324 file_DYN_beg = f'{TmpDir}/ICO_{JobName}_{YearRes}1231_restart.nc' 325 liste_beg.append (file_DYN_beg) 326 config['Files']['file_DYN_beg'] = file_DYN_beg 327 328 if file_DYN_end == None : 329 if LMDZ : 330 file_DYN_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restart.nc' 331 if ICO : 332 file_DYN_end = f'{TmpDir}/ICO_{JobName}_{YearEnd}1231_restart.nc' 333 liste_end.append ( file_DYN_end ) 334 config['Files']['file_DYN_end'] = file_DYN_end 335 336 if file_SRF_beg == None : 337 file_SRF_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 338 config['Files']['file_SRF_beg'] = file_SRF_beg 339 if file_SRF_end == None : 340 file_SRF_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 341 config['Files']['file_SRF_end'] = file_SRF_end 342 343 liste_beg.append ( file_SRF_beg ) 344 liste_end.append ( file_SRF_end ) 345 244 346 echo ( f'{file_ATM_beg}') 245 347 echo ( f'{file_ATM_end}') 348 echo ( f'{file_DYN_beg}') 349 echo ( f'{file_DYN_end}') 246 350 echo ( f'{file_SRF_beg}') 247 351 echo ( f'{file_SRF_end}') 248 352 249 if Routing == 'SIMPLE' : 250 file_RUN_beg = f'SRF_{JobName}_{YearRes}1231_routing_restart.nc' 251 file_RUN_end = f'SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 353 if Routing == 'SIMPLE' : 354 if file_RUN_beg == None : 355 file_RUN_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc' 356 config['Files']['file_RUN_beg'] = file_RUN_beg 357 if file_RUN_end == None : 358 file_RUN_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 359 config['Files']['file_RUN_end'] = file_RUN_end 360 252 361 liste_beg.append ( file_RUN_beg ) 253 362 liste_end.append ( file_RUN_end ) … … 257 366 echo ('\nExtract restart files from tar : ATM, ICO and SRF') 258 367 for resFile in liste_beg : 259 if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 260 command = f'cd {TmpDir} ; tar xf {file_restart_beg} {resFile}' 368 if os.path.exists ( os.path.join (TmpDir, resFile) ) : 369 echo ( f'file found : {resFile}' ) 370 else : 371 base_resFile = os.path.basename (resFile) 372 command = f'cd {TmpDir} ; tar xf {file_restart_beg} {base_resFile}' 261 373 echo ( command ) 262 374 os.system ( command ) 263 375 264 376 for resFile in liste_end : 265 if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 266 command = f'cd {TmpDir} ; tar xf {file_restart_end} {resFile}' 377 if os.path.exists ( os.path.join (TmpDir, resFile) ) : 378 echo ( f'file found : {resFile}' ) 379 else : 380 base_resFile = os.path.basename (resFile) 381 command = f'cd {TmpDir} ; tar xf {file_restart_end} {base_resFile}' 267 382 echo ( command ) 268 383 os.system ( command ) 269 270 config['Files']['file_ATM_beg'] = file_ATM_beg 271 config['Files']['file_ATM_end'] = file_ATM_end 272 config['Files']['file_SRF_beg'] = file_SRF_beg 273 config['Files']['file_SRF_end'] = file_SRF_end 274 if Routing == 'SIMPLE' : 275 config['Files']['file_RUN_beg'] = file_RUN_beg 276 config['Files']['file_RUN_end'] = file_RUN_end 277 config['Files']['file_DYN_beg'] = file_DYN_beg 278 config['Files']['file_DYN_end'] = file_DYN_end 279 384 280 385 echo ('\nOpening ATM SRF and ICO restart files') 281 386 d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze() … … 287 392 288 393 for var in d_SRF_beg.variables : 289 d_SRF_beg[var] = d_SRF_beg[var].where ( 290 d_SRF_end[var] = d_SRF_end[var].where ( 291 292 if ICO:394 d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 395 d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 396 397 if Routing == 'SIMPLE' : 293 398 d_RUN_beg = xr.open_dataset ( os.path.join (TmpDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze() 294 399 d_RUN_end = xr.open_dataset ( os.path.join (TmpDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze() … … 304 409 echo ( file_RUN_end ) 305 410 411 ## 412 config_out = open (FullIniFile, 'w') 413 config.write (config_out ) 414 config_out.close () 415 416 def Ksum (tab) : 417 ''' 418 Kahan summation algorithm, also known as compensated summation 419 https://en.wikipedia.org/wiki/Kahan_summation_algorithm 420 ''' 421 Ksum = 0.0 # Prepare the accumulator. 422 comp = 0.0 # A running compensation for lost low-order bits. 423 424 for xx in np.where ( np.isnan(tab), 0., tab ) : 425 yy = xx - comp # comp is zero the first time around. 426 tt = Ksum + yy # Alas, sum is big, y small, so low-order digits of y are lost. 427 comp = (tt - Ksum) - yy # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy) 428 Ksum = tt # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers! 429 # Next time around, the lost low part will be added to y in a fresh attempt. 430 return Ksum 431 432 def Ssum (tab) : 433 ''' 434 Precision summation by sorting by absolute value 435 ''' 436 Ssum = np.sum ( tab[np.argsort(np.abs(tab))] ) 437 return Ssum 438 439 def KSsum (tab) : 440 ''' 441 Precision summation by sorting by absolute value, and applying Kahan summation 442 ''' 443 KSsum = Ksum ( tab[np.argsort(np.abs(tab))] ) 444 return KSsum 445 446 # Choosing algorithm 447 Psum = KSsum 448 306 449 def ATM_stock_int (stock) : 307 450 '''Integrate stock on atmosphere grid''' 308 ATM_stock_int = np.sum ( np.sort ( (stock * DYN_aire).to_masked_array().ravel()) )451 ATM_stock_int = Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 309 452 return ATM_stock_int 310 453 311 454 def ATM_flux_int (flux) : 312 455 '''Integrate flux on atmosphere grid''' 313 ATM_stock_int = np.sum ( np.sort ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel()) ) 456 ATM_flux_int = Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 457 return ATM_flux_int 458 459 def SRF_stock_int (stock) : 460 '''Integrate stock on land grid''' 461 ATM_stock_int = Ksum ( ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 314 462 return ATM_stock_int 463 464 def SRF_flux_int (flux) : 465 '''Integrate flux on land grid''' 466 SRF_flux_int = Psum ( (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 467 return SRF_flux_int 315 468 316 469 def ONE_stock_int (stock) : 317 470 '''Sum stock ''' 318 ONE_stock_int = np.sum ( np.sort ( (stock).to_masked_array().ravel()) )471 ONE_stock_int = Psum ( stock.to_masked_array().ravel() ) 319 472 return ONE_stock_int 320 473 321 474 def ONE_flux_int (flux) : 322 475 '''Sum flux ''' 323 ONE_flux_int = np.sum ( np.sort ( (flux * dtime_per_sec ).to_masked_array().ravel()) )476 ONE_flux_int = Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 324 477 return ONE_flux_int 478 479 def kg2Sv (val, rho=ATM_rho) : 480 '''From kg to Sverdrup''' 481 return val/dtime_sec*1.0e-6/rho 482 483 def kg2myear (val, rho=ATM_rho) : 484 '''From kg to m/year''' 485 return val/ATM_aire_sea_tot/rho/NbYear 325 486 326 487 # ATM grid with cell surfaces … … 329 490 file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 330 491 config['Files']['file_ATM_aire'] = file_ATM_aire 331 echo ( 'Aire sur grille reguliere :', file_ATM_aire)492 echo ( f'Aire sur grille reguliere : {file_ATM_aire}' ) 332 493 d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 333 494 ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True ) 334 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 335 ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_lic'][0] ) 495 ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 496 ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] ) 497 ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] ) 498 ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] ) 499 ATM_fsea = ATM_foce + ATM_fsic 500 ATM_flnd = ATM_fter + ATM_flic 336 501 SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 502 #SRF_aire = ATM_aire * lmdz.geo2point (d_SRF_his ['Contfrac'] ) 503 #SRF_aire = ATM_aire * ATM_fter 337 504 338 505 if LMDZ : … … 340 507 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 341 508 ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 342 SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 509 #SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 510 SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'] ) 343 511 344 512 SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.) … … 353 521 DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic'] 354 522 DYN_flnd = 1.0 - DYN_fsea 523 DYN_fter = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'}) 524 DYN_flic = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'}) 355 525 356 526 if LMDZ : … … 358 528 DYN_fsea = ATM_fsea 359 529 DYN_flnd = ATM_flnd 530 DYN_fter = d_ATM_beg['FTER'] 531 DYN_flic = d_ATM_beg['FTER'] 532 533 534 DYN_aire_fter = DYN_aire * DYN_fter 360 535 361 536 #if LMDZ : … … 374 549 raise Exception ('Erreur surface interpolee sur grille reguliere') 375 550 376 echo ( '\n ------------------------------------------------------------------------------------' )377 echo ( '-- ATM changes in stores' )551 echo ( '\n====================================================================================' ) 552 echo ( f'-- ATM changes in stores -- {Title} ' ) 378 553 379 554 #-- Change in precipitable water from the atmosphere daily and monthly files … … 386 561 # Surface pressure 387 562 if ICO : 388 DYN_ps_beg = d_DYN_beg['ps'] 389 DYN_ps_end = d_DYN_end['ps'] 390 563 DYN_psol_beg = d_DYN_beg['ps'] 564 DYN_psol_end = d_DYN_end['ps'] 391 565 if LMDZ : 392 DYN_ps_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 393 DYN_ps_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 394 395 # 3D Pressure 396 DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg 397 DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end 398 399 if ICO : klevp1, cell_mesh = DYN_p_beg.shape 400 if LMDZ : klevp1, points_physiques = DYN_p_beg.shape 566 DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 567 DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 568 569 # 3D Pressure at the interface layers (not scalar points) 570 DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg 571 DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end 572 573 klevp1 = ATM_Bhyb.shape[-1] 574 if ICO : cell_mesh = DYN_psol_beg.shape[-1] 575 if LMDZ : points_physiques = DYN_psol_beg.shape[-1] 401 576 klev = klevp1 - 1 402 577 403 # Layer thickness 578 # Layer thickness (pressure) 404 579 if ICO : 405 DYN_d sigma_beg = xr.DataArray ( np.empty( (klev, cell_mesh )), dims=('sigs', 'cell_mesh' ), coords=(np.arange(klev), np.arange(cell_mesh)) )406 DYN_d sigma_end = xr.DataArray ( np.empty( (klev, cell_mesh )), dims=('sigs', 'cell_mesh' ), coords=(np.arange(klev), np.arange(cell_mesh)) )580 DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell_mesh )), dims=('sigs', 'cell_mesh' ), coords=(np.arange(klev), np.arange(cell_mesh) ) ) 581 DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell_mesh )), dims=('sigs', 'cell_mesh' ), coords=(np.arange(klev), np.arange(cell_mesh) ) ) 407 582 if LMDZ : 408 DYN_d sigma_beg = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) )409 DYN_d sigma_end = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) )583 DYN_dpres_beg = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 584 DYN_dpres_end = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 410 585 411 586 for k in np.arange (klevp1-1) : 412 DYN_d sigma_beg[k,:] = DYN_p_beg[k,:] - DYN_p_beg[k+1,:]413 DYN_d sigma_end[k,:] = DYN_p_end[k,:] - DYN_p_end[k+1,:]587 DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] 588 DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:] 414 589 415 590 ##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases … … 429 604 DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) ).rename ( {'lev':'sigs'} ) 430 605 431 # Integral 432 DYN_mas_wat_beg = ATM_stock_int (DYN_dsigma_beg * DYN_wat_beg) / Grav 433 DYN_mas_wat_end = ATM_stock_int (DYN_dsigma_end * DYN_wat_end) / Grav 434 606 # Compute water content : vertical and horizontal integral 607 DYN_mas_wat_beg = ATM_stock_int (DYN_dpres_beg * DYN_wat_beg) / Grav 608 DYN_mas_wat_end = ATM_stock_int (DYN_dpres_end * DYN_wat_end) / Grav 609 610 # Variation of water content 435 611 dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 436 612 437 echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' ) 613 # \([a-z,A-Z,_]*\)/dtime_sec\*1e-9 kg2Sv(\1) 614 # \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear kg2myear(\1) 615 616 def var2prt (var, small=False) : 617 if small : return var , kg2Sv(var)*1000., kg2myear(var)*1000. 618 else : return var , kg2Sv(var) , kg2myear(var) 619 620 def prtFlux (Desc, var, Form='F', small=False) : 621 if small : 622 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 623 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 624 else : 625 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 626 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 627 echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) ) 628 return None 629 630 echo ( f'\nVariation du contenu en eau atmosphere (dynamique) -- {Title} ' ) 631 echo ( '------------------------------------------------------------------------------------' ) 438 632 echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 439 echo ( 'dMass (atm) = {:12.3e} kg '.format (dDYN_mas_wat) ) 440 echo ( 'dMass (atm) = {:12.3e} Sv '.format (dDYN_mas_wat/dtime_sec*1.e-6/ATM_rho) ) 441 echo ( 'dMass (atm) = {:12.3e} m '.format (dDYN_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 633 prtFlux ( 'dMass (atm) ', dDYN_mas_wat, 'e', True ) 442 634 443 635 ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER']+d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']+d_ATM_beg['SNOW03']*d_ATM_beg['FOCE']+d_ATM_beg['SNOW04']*d_ATM_beg['FSIC'] … … 499 691 dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg 500 692 501 echo ( '\nVariation du contenu en neige atmosphere (calottes)' ) 693 dATM_mas_sno_2 = ATM_stock_int ( ATM_sno_end - ATM_sno_beg ) 694 695 echo ( f'\nVariation du contenu en neige atmosphere (calottes) -- {Title} ' ) 696 echo ( '------------------------------------------------------------------------------------' ) 502 697 echo ( 'ATM_mas_sno_beg = {:12.6e} kg | ATM_mas_sno_end = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 503 echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_sno ))504 echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1e-6/ICE_rho_ice))505 echo ( 'dMass (neige atm) = {:12.3e} m '.format (dATM_mas_sno/ATM_aire_sea_tot/ATM_rho) ) 506 507 echo ( ' \nVariation du contenu humidite du sol' )698 prtFlux ( 'dMass (neige atm) ', dATM_mas_sno , 'e', True ) 699 prtFlux ( 'dMass (neige atm) ', dATM_mas_sno_2, 'e', True ) 700 701 echo ( f'\nVariation du contenu humidite du sol -- {Title} ' ) 702 echo ( '------------------------------------------------------------------------------------' ) 508 703 echo ( 'ATM_mas_qs_beg = {:12.6e} kg | ATM_mas_qs_end = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) ) 509 echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_qs ) ) 510 echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_qs/dtime_sec*1e-6/ATM_rho) ) 511 echo ( 'dMass (neige atm) = {:12.3e} m '.format (dATM_mas_qs/ATM_aire_sea_tot/ATM_rho) ) 512 513 echo ( '\nVariation du contenu en eau+neige atmosphere ' ) 514 echo ( 'dMass (eau + neige atm) = {:12.3e} kg '.format ( dDYN_mas_wat + dATM_mas_sno) ) 515 echo ( 'dMass (eau + neige atm) = {:12.3e} Sv '.format ( (dDYN_mas_wat + dATM_mas_sno)/dtime_sec*1E-6/ATM_rho) ) 516 echo ( 'dMass (eau + neige atm) = {:12.3e} m '.format ( (dDYN_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot/ATM_rho) ) 517 518 echo ( '\n------------------------------------------------------------------------------------' ) 519 echo ( '-- SRF changes ' ) 520 521 # dSoilHum_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=(maxvegetfrac*DelSoilMoist_daily)*Contfrac' ${file} -gridarea ${file}` 522 # dInterce_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelIntercept_daily*Contfrac' ${file} -gridarea ${file}` 523 # dSWE_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelSWE_daily*Contfrac' ${file} -gridarea ${file}` 524 # dStream_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delstreamr_daily*Contfrac' ${file} -gridarea ${file}` 525 # dFastR_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfastr_daily*Contfrac' ${file} -gridarea ${file}` 526 # dSlowR_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delslowr_daily*Contfrac' ${file} -gridarea ${file}` 527 # dFlood_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfloodr_daily*Contfrac' ${file} -gridarea ${file}` 528 # dPond_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delpondr_daily*Contfrac' ${file} -gridarea ${file}` 529 # dLake_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=dellakevol_daily*Contfrac' ${file} -gridarea ${file}` 530 #echo 'dSTOCK (Sv) Soil Intercept SWE Stream FastR SlowR Lake Pond Flood=' $dSoilHum_in_Sv, $dInterce_in_Sv, $dSWE_in_Sv, $dStream_in_Sv, $dFastR_in_Sv, $dSlowR_in_Sv, $dLake_in_Sv, $dPond_in_Sv, $dFlood_in_Sv >> ${fileout} 531 #dSRF_tot=`python -c "print $dSoilHum_in_Sv+$dInterce_in_Sv+$dSWE_in_Sv+$dStream_in_Sv+$dFastR_in_Sv+$dSlowR_in_Sv+$dLake_in_Sv+$dPond_in_Sv+$dFlood_in_Sv"` 532 #echo 'dSTOCK (Sv) total='${dSRF_tot} >> ${fileout} 533 704 prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True ) 705 706 echo ( f'\nVariation du contenu en eau+neige atmosphere -- {Title} ' ) 707 echo ( '------------------------------------------------------------------------------------' ) 708 prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True) 709 710 echo ( '\n====================================================================================' ) 711 echo ( f'-- SRF changes -- {Title} ' ) 534 712 535 713 if Routing == 'SIMPLE' : 536 RUN_mas_wat_beg = ONE_stock_int ( d_RUN_beg ['fast_reservoir'] + d_RUN_beg ['slow_reservoir'] + d_RUN_beg ['stream_reservoir'] ) 537 RUN_mas_wat_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] + d_RUN_end ['slow_reservoir'] + d_RUN_end ['stream_reservoir'] ) 538 539 if Routing == 'ORCHIDEE' : 540 RUN_mas_wat_beg = ONE_stock_int ( d_SRF_beg['fastres'] + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \ 541 + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] ) 542 RUN_mas_wat_end = ONE_stock_int ( d_SRF_end['fastres'] + d_SRF_end['slowres'] + d_SRF_end['streamres'] \ 543 + d_SRF_end['floodres'] + d_SRF_end['lakeres'] + d_SRF_end['pondres'] ) 544 545 dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg 546 547 echo ( '\nWater content in routing ' ) 714 RUN_mas_wat_fast_beg = ONE_stock_int ( d_RUN_beg ['fast_reservoir'] ) 715 RUN_mas_wat_slow_beg = ONE_stock_int ( d_RUN_beg ['slow_reservoir'] ) 716 RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 717 RUN_mas_wat_flood_beg = 0.0 718 RUN_mas_wat_lake_beg = 0.0 719 RUN_mas_wat_pond_beg = 0.0 720 721 RUN_mas_wat_fast_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] ) 722 RUN_mas_wat_slow_end = ONE_stock_int ( d_RUN_end ['slow_reservoir'] ) 723 RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 724 RUN_mas_wat_flood_end = 0.0 725 RUN_mas_wat_lake_end = 0.0 726 RUN_mas_wat_pond_end = 0.0 727 728 if Routing == 'SECHIBA' : 729 RUN_mas_wat_fast_beg = ONE_stock_int ( d_SRF_beg ['fastres'] ) 730 RUN_mas_wat_slow_beg = ONE_stock_int ( d_SRF_beg ['slowres'] ) 731 RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 732 RUN_mas_wat_flood_beg = ONE_stock_int ( d_SRF_beg ['floodres'] ) 733 RUN_mas_wat_lake_beg = ONE_stock_int ( d_SRF_beg ['lakeres'] ) 734 RUN_mas_wat_pond_beg = ONE_stock_int ( d_SRF_beg ['pondres'] ) 735 736 737 RUN_mas_wat_fast_end = ONE_stock_int ( d_SRF_end ['fastres'] ) 738 RUN_mas_wat_slow_end = ONE_stock_int ( d_SRF_end ['slowres'] ) 739 RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 740 RUN_mas_wat_flood_end = ONE_stock_int ( d_SRF_end ['floodres'] ) 741 RUN_mas_wat_lake_end = ONE_stock_int ( d_SRF_end ['lakeres'] ) 742 RUN_mas_wat_pond_end = ONE_stock_int ( d_SRF_end ['pondres'] ) 743 744 RUN_mas_wat_beg = RUN_mas_wat_fast_beg + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg 745 RUN_mas_wat_end = RUN_mas_wat_fast_end + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end 746 747 dRUN_mas_wat_fast = RUN_mas_wat_fast_end - RUN_mas_wat_fast_beg 748 dRUN_mas_wat_slow = RUN_mas_wat_slow_end - RUN_mas_wat_slow_beg 749 dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg 750 dRUN_mas_wat_flood = RUN_mas_wat_flood_end - RUN_mas_wat_flood_beg 751 dRUN_mas_wat_lake = RUN_mas_wat_lake_end - RUN_mas_wat_lake_beg 752 dRUN_mas_wat_pond = RUN_mas_wat_pond_end - RUN_mas_wat_pond_beg 753 754 dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg 755 756 echo ( f'\nLes differents reservoirs -- {Title} ') 757 echo ( '------------------------------------------------------------------------------------' ) 758 echo ( 'RUN_mas_wat_fast_beg = {:12.6e} kg | RUN_mas_wat_fast_end = {:12.6e} kg '.format (RUN_mas_wat_fast_beg , RUN_mas_wat_fast_end ) ) 759 echo ( 'RUN_mas_wat_slow_beg = {:12.6e} kg | RUN_mas_wat_slow_end = {:12.6e} kg '.format (RUN_mas_wat_slow_beg , RUN_mas_wat_slow_end ) ) 760 echo ( 'RUN_mas_wat_stream_beg = {:12.6e} kg | RUN_mas_wat_stream_end = {:12.6e} kg '.format (RUN_mas_wat_stream_beg, RUN_mas_wat_stream_end ) ) 761 echo ( 'RUN_mas_wat_flood_beg = {:12.6e} kg | RUN_mas_wat_flood_end = {:12.6e} kg '.format (RUN_mas_wat_flood_beg , RUN_mas_wat_flood_end ) ) 762 echo ( 'RUN_mas_wat_lake_beg = {:12.6e} kg | RUN_mas_wat_lake_end = {:12.6e} kg '.format (RUN_mas_wat_lake_beg , RUN_mas_wat_lake_end ) ) 763 echo ( 'RUN_mas_wat_pond_beg = {:12.6e} kg | RUN_mas_wat_pond_end = {:12.6e} kg '.format (RUN_mas_wat_pond_beg , RUN_mas_wat_pond_end ) ) 764 765 echo ( '------------------------------------------------------------------------------------' ) 766 767 prtFlux ( 'dMass (fast) ', dRUN_mas_wat_fast , 'e', True ) 768 prtFlux ( 'dMass (slow) ', dRUN_mas_wat_slow , 'e', True ) 769 prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True ) 770 prtFlux ( 'dMass (flood) ', dRUN_mas_wat_flood , 'e', True ) 771 prtFlux ( 'dMass (lake) ', dRUN_mas_wat_lake , 'e', True ) 772 prtFlux ( 'dMass (pond) ', dRUN_mas_wat_pond , 'e', True ) 773 prtFlux ( 'dMass (all) ', dRUN_mas_wat , 'e', True ) 774 775 echo ( f'\nWater content in routing -- {Title} ' ) 776 echo ( '------------------------------------------------------------------------------------' ) 548 777 echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) ) 549 echo ( 'dMass (routing) = {:12.3e} kg '.format(dRUN_mas_wat) ) 550 echo ( 'dMass (routing) = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) ) 551 echo ( 'dMass (routing) = {:12.3e} m '.format(dRUN_mas_wat/ATM_aire_sea_tot*1E-3) ) 552 553 print ('Reading SRF restart') 554 tot_watveg_beg = d_SRF_beg['tot_watveg_beg'] ; tot_watveg_beg = tot_watveg_beg .where (tot_watveg_beg < 1E10, 0.) 555 tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; tot_watsoil_beg = tot_watsoil_beg.where (tot_watsoil_beg< 1E10, 0.) 556 snow_beg = d_SRF_beg['snow_beg'] ; snow_beg = snow_beg .where (snow_beg < 1E10, 0.) 557 558 tot_watveg_end = d_SRF_end['tot_watveg_beg'] ; tot_watveg_end = tot_watveg_end .where (tot_watveg_end < 1E10, 0.) 559 tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; tot_watsoil_end = tot_watsoil_end.where (tot_watsoil_end< 1E10, 0.) 560 snow_end = d_SRF_end['snow_beg'] ; snow_end = snow_end .where (snow_end < 1E10, 0.) 778 prtFlux ( 'dMass (routing) ', dRUN_mas_wat ) 779 780 echo ( '\n====================================================================================' ) 781 print (f'Reading SRF restart') 782 SRF_tot_watveg_beg = d_SRF_beg['tot_watveg_beg'] ; SRF_tot_watveg_beg = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg < 1E15, 0.) 783 SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg< 1E15, 0.) 784 SRF_snow_beg = d_SRF_beg['snow_beg'] ; SRF_snow_beg = SRF_snow_beg .where (SRF_snow_beg < 1E15, 0.) 785 SRF_lakeres_beg = d_SRF_beg['lakeres'] ; SRF_lakeres_beg = SRF_lakeres_beg .where (SRF_lakeres_beg < 1E15, 0.) 786 787 SRF_tot_watveg_end = d_SRF_end['tot_watveg_beg'] ; SRF_tot_watveg_end = SRF_tot_watveg_end .where (SRF_tot_watveg_end < 1E15, 0.) 788 SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end< 1E15, 0.) 789 SRF_snow_end = d_SRF_end['snow_beg'] ; SRF_snow_end = SRF_snow_end .where (SRF_snow_end < 1E15, 0.) 790 SRF_lakeres_end = d_SRF_end['lakeres'] ; SRF_lakeres_end = SRF_lakeres_end .where (SRF_lakeres_end < 1E15, 0.) 561 791 562 792 if LMDZ : 563 tot_watveg_beg = lmdz.geo2point (tot_watveg_beg) 564 tot_watsoil_beg = lmdz.geo2point (tot_watsoil_beg) 565 snow_beg = lmdz.geo2point (snow_beg) 566 tot_watveg_end = lmdz.geo2point (tot_watveg_end) 567 tot_watsoil_end = lmdz.geo2point (tot_watsoil_end) 568 snow_end = lmdz.geo2point (snow_end) 793 SRF_tot_watveg_beg = lmdz.geo2point (SRF_tot_watveg_beg ) 794 SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg) 795 SRF_snow_beg = lmdz.geo2point (SRF_snow_beg ) 796 SRF_lakeres_beg = lmdz.geo2point (SRF_lakeres_beg ) 797 SRF_tot_watveg_end = lmdz.geo2point (SRF_tot_watveg_end ) 798 SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end) 799 SRF_snow_end = lmdz.geo2point (SRF_snow_end ) 800 SRF_lakeres_end = lmdz.geo2point (SRF_lakeres_end ) 801 802 if ICO : 803 SRF_tot_watveg_beg = SRF_tot_watveg_beg .rename ( {'y':'cell_mesh'} ) 804 SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.rename ( {'y':'cell_mesh'} ) 805 SRF_snow_beg = SRF_snow_beg .rename ( {'y':'cell_mesh'} ) 806 SRF_lakeres_beg = SRF_lakeres_beg .rename ( {'y':'cell_mesh'} ) 807 SRF_tot_watveg_end = SRF_tot_watveg_end .rename ( {'y':'cell_mesh'} ) 808 SRF_tot_watsoil_end = SRF_tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 809 SRF_snow_end = SRF_snow_end .rename ( {'y':'cell_mesh'} ) 810 SRF_lakeres_end = SRF_lakeres_end .rename ( {'y':'cell_mesh'} ) 569 811 570 812 # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 571 813 572 SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 573 SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 574 575 if ICO : 576 tot_watveg_beg = tot_watveg_beg .rename ( {'y':'cell_mesh'} ) 577 tot_watsoil_beg = tot_watsoil_beg.rename ( {'y':'cell_mesh'} ) 578 snow_beg = snow_beg .rename ( {'y':'cell_mesh'} ) 579 tot_watveg_end = tot_watveg_end .rename ( {'y':'cell_mesh'} ) 580 tot_watsoil_end = tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 581 snow_end = snow_end .rename ( {'y':'cell_mesh'} ) 582 SRF_wat_beg = SRF_wat_beg .rename ( {'y':'cell_mesh'} ) 583 SRF_wat_end = SRF_wat_end .rename ( {'y':'cell_mesh'} ) 584 814 SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg 815 SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 816 817 echo ( '\n====================================================================================' ) 585 818 print ('Computing integrals') 586 819 587 print ( ' 1/6', end='' ) ; sys.stdout.flush () 588 SRF_mas_watveg_beg = ATM_stock_int ( tot_watveg_beg ) 589 print ( ' 2/6', end='' ) ; sys.stdout.flush () 590 SRF_mas_watsoil_beg = ATM_stock_int ( tot_watsoil_beg ) 591 print ( ' 3/6', end='' ) ; sys.stdout.flush () 592 SRF_mas_snow_beg = ATM_stock_int ( snow_beg ) 593 print ( ' 4/6', end='' ) ; sys.stdout.flush () 594 SRF_mas_watveg_end = ATM_stock_int ( tot_watveg_end ) 595 print ( ' 5/6', end='' ) ; sys.stdout.flush () 596 SRF_mas_watsoil_end = ATM_stock_int ( tot_watsoil_end ) 597 print ( ' 6/6', end='' ) ; sys.stdout.flush () 598 SRF_mas_snow_end = ATM_stock_int ( snow_end ) 820 print ( ' 1/8', end='' ) ; sys.stdout.flush () 821 SRF_mas_watveg_beg = SRF_stock_int ( SRF_tot_watveg_beg ) 822 print ( ' 2/8', end='' ) ; sys.stdout.flush () 823 SRF_mas_watsoil_beg = SRF_stock_int ( SRF_tot_watsoil_beg ) 824 print ( ' 3/8', end='' ) ; sys.stdout.flush () 825 SRF_mas_snow_beg = SRF_stock_int ( SRF_snow_beg ) 826 print ( ' 4/8', end='' ) ; sys.stdout.flush () 827 SRF_mas_lake_beg = ONE_stock_int ( SRF_lakeres_beg ) 828 print ( ' 5/8', end='' ) ; sys.stdout.flush () 829 830 SRF_mas_watveg_end = SRF_stock_int ( SRF_tot_watveg_end ) 831 print ( ' 6/8', end='' ) ; sys.stdout.flush () 832 SRF_mas_watsoil_end = SRF_stock_int ( SRF_tot_watsoil_end ) 833 print ( ' 7/8', end='' ) ; sys.stdout.flush () 834 SRF_mas_snow_end = SRF_stock_int ( SRF_snow_end ) 835 print ( ' 8/8', end='' ) ; sys.stdout.flush () 836 SRF_mas_lake_end = ONE_stock_int ( SRF_lakeres_end ) 837 599 838 print (' -- ') ; sys.stdout.flush () 600 839 601 dSRF_mas_watveg = SRF_mas_watveg_end - SRF_mas_watveg_beg 602 dSRF_mas_watsoil = SRF_mas_watsoil_end - SRF_mas_watsoil_beg 603 dSRF_mas_snow = SRF_mas_snow_end - SRF_mas_snow_beg 604 605 echo ('\nLes differents reservoirs') 840 dSRF_mas_watveg = SRF_mas_watveg_end - SRF_mas_watveg_beg 841 dSRF_mas_watsoil = SRF_mas_watsoil_end - SRF_mas_watsoil_beg 842 dSRF_mas_snow = SRF_mas_snow_end - SRF_mas_snow_beg 843 dSRF_mas_lake = SRF_mas_lake_end - SRF_mas_lake_beg 844 845 846 847 echo ( '------------------------------------------------------------------------------------' ) 848 echo ( f'\nLes differents reservoirs -- {Title} ') 606 849 echo ( 'SRF_mas_watveg_beg = {:12.6e} kg | SRF_mas_watveg_end = {:12.6e} kg '.format (SRF_mas_watveg_beg , SRF_mas_watveg_end ) ) 607 850 echo ( 'SRF_mas_watsoil_beg = {:12.6e} kg | SRF_mas_watsoil_end = {:12.6e} kg '.format (SRF_mas_watsoil_beg , SRF_mas_watsoil_end ) ) 608 851 echo ( 'SRF_mas_snow_beg = {:12.6e} kg | SRF_mas_snow_end = {:12.6e} kg '.format (SRF_mas_snow_beg , SRF_mas_snow_end ) ) 609 610 echo ( 'dMass (watveg) = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watveg , dSRF_mas_watveg /dtime_sec*1E-9, dSRF_mas_watveg /ATM_aire_sea_tot*1E-3) ) 611 echo ( 'dMass (watsoil) = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watsoil , dSRF_mas_watsoil /dtime_sec*1E-9, dSRF_mas_watsoil /ATM_aire_sea_tot*1E-3) ) 612 echo ( 'dMass (sno) = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_snow , dSRF_mas_snow /dtime_sec*1E-9, dSRF_mas_snow /ATM_aire_sea_tot*1E-3) ) 852 echo ( 'SRF_mas_lake_beg = {:12.6e} kg | SRF_mas_lake_end = {:12.6e} kg '.format (SRF_mas_lake_beg , SRF_mas_lake_end ) ) 853 854 prtFlux ('dMass (watveg) ', dSRF_mas_watveg , 'e' , True) 855 prtFlux ('dMass (watsoil)', dSRF_mas_watsoil, 'e' , True) 856 prtFlux ('dMass (snow) ', dSRF_mas_snow , 'e' , True) 857 prtFlux ('dMass (lake) ', dSRF_mas_lake , 'e' , True) 613 858 614 859 SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg 615 860 SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end 616 dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg617 618 echo ( ' \nWater content in surface' )619 echo ( 'SRF_mas_wat_beg = {:12.6e} kg | SRF_mas_wat_end = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end))620 echo ( ' dMass (water srf) = {:12.3e} kg '.format (dSRF_mas_wat) )621 echo ( 'dMass (water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-6/ATM_rho))622 echo ( 'dMass (water srf) = {:12.3e} m '.format (dSRF_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 623 624 echo ( ' \nWater content in ATM + SRF + RUN' )861 dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg 862 863 echo ( '------------------------------------------------------------------------------------' ) 864 echo ( f'Water content in surface -- {Title} ' ) 865 echo ( 'SRF_mas_wat_beg = {:12.6e} kg | SRF_mas_wat_end = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 866 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 867 868 echo ( '------------------------------------------------------------------------------------' ) 869 echo ( 'Water content in ATM + SRF + RUN + LAKE' ) 625 870 echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '. 626 format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg, 627 DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) ) 628 echo ( 'dMass (water atm+srf+run) = {:12.6e} kg '.format ( dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat) ) 629 echo ( 'dMass (water atm+srf+run) = {:12.3e} Sv '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/dtime_sec*1E-6/ATM_rho) ) 630 echo ( 'dMass (water atm+srf+run) = {:12.3e} m '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/ATM_aire_sea_tot/ATM_rho) ) 631 632 echo ( '\n------------------------------------------------------------------------------------' ) 633 echo ( '-- ATM Fluxes ' ) 871 format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg , 872 DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) ) 873 prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True) 874 875 echo ( '\n====================================================================================' ) 876 echo ( f'-- ATM Fluxes -- {Title} ' ) 634 877 635 878 ATM_wbilo_oce = lmdz.geo2point ( d_ATM_his ['wbilo_oce'] ) … … 643 886 ATM_snowf = lmdz.geo2point ( d_ATM_his ['snow'] ) 644 887 ATM_evap = lmdz.geo2point ( d_ATM_his ['evap'] ) 645 ATM_bilo = ATM_wbilo_oce + ATM_wbilo_ter + ATM_wbilo_lic + ATM_wbilo_sic 888 ATM_wbilo = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic 889 ATM_emp = ATM_evap - ATM_precip 646 890 647 891 RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'] ) … … 651 895 RUN_riversret = lmdz.geo2point ( d_RUN_his ['riversret'] ) 652 896 653 SRF_evap = lmdz.geo2point ( d_SRF_his ['evap'] ) 897 RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'] ) 898 RUN_riverflow_cpl = lmdz.geo2point ( d_RUN_his ['riverflow_cpl'] ) 899 900 SRF_rain = lmdz.geo2point ( d_SRF_his ['rain'] ) 901 SRF_evap = lmdz.geo2point ( d_SRF_his ['evap'] ) 654 902 SRF_snowf = lmdz.geo2point ( d_SRF_his ['snowf'] ) 655 SRF_TWS = lmdz.geo2point ( d_SRF_his ['TWS'] )656 903 SRF_subli = lmdz.geo2point ( d_SRF_his ['subli'] ) 657 904 SRF_transpir = lmdz.geo2point ( np.sum(d_SRF_his ['transpir'], axis=1) ) ; SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 658 659 def mmd2SI ( Var) : 905 SRF_emp = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 906 907 ## Correcting units of SECHIBA variables 908 def mmd2SI ( Var ) : 660 909 '''Change unit from mm/d or m^3/s to kg/s if needed''' 661 910 if 'units' in VarT.attrs : … … 665 914 VarT.values = VarT.values * ATM_rho * (1e-3/86400.) ; VarT.attrs['units'] = 'kg/s' 666 915 667 for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow'] :916 for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 668 917 VarT = locals()['RUN_' + var] 669 918 mmd2SI (VarT) 670 919 671 for var in ['evap', 'snowf', ' TWS', 'subli', 'transpir'] :920 for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] : 672 921 VarT = locals()['SRF_' + var] 673 922 mmd2SI (VarT) 674 923 924 675 925 RUN_input = RUN_runoff + RUN_drainage 676 926 RUN_output = RUN_coastalflow + RUN_riverflow … … 678 928 ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic 679 929 680 ATM_flx_oce = ATM_flux_int ( ATM_wbilo_oce ) 681 ATM_flx_sic = ATM_flux_int ( ATM_wbilo_sic ) 682 ATM_flx_sea = ATM_flux_int ( ATM_wbilo_sea ) 683 ATM_flx_ter = ATM_flux_int ( ATM_wbilo_ter ) 684 ATM_flx_lic = ATM_flux_int ( ATM_wbilo_lic ) 685 ATM_flx_calving = ATM_flux_int ( ATM_fqcalving ) 686 ATM_flx_qfonte = ATM_flux_int ( ATM_fqfonte ) 687 ATM_flx_precip = ATM_flux_int ( ATM_precip ) 688 ATM_flx_snowf = ATM_flux_int ( ATM_snowf ) 689 ATM_flx_evap = ATM_flux_int ( ATM_evap ) 690 ATM_flx_runlic = ATM_flux_int ( ATM_runofflic ) 691 692 RUN_flx_coastal = ONE_flux_int ( RUN_coastalflow) 693 RUN_flx_river = ONE_flux_int ( RUN_riverflow ) 694 RUN_flx_drainage = ATM_flux_int ( RUN_drainage ) 695 RUN_flx_riversret = ATM_flux_int ( RUN_riversret ) 696 RUN_flx_runoff = ATM_flux_int ( RUN_runoff ) 697 RUN_flx_input = ATM_flux_int ( RUN_input ) 698 RUN_flx_output = ONE_flux_int ( RUN_output ) 699 700 ATM_flx_emp = ATM_flx_evap - ATM_flx_precip 701 ATM_flx_bilo = ATM_flx_oce + ATM_flx_sic + ATM_flx_ter + ATM_flx_lic 930 ATM_flx_oce = ATM_flux_int ( ATM_wbilo_oce ) 931 ATM_flx_sic = ATM_flux_int ( ATM_wbilo_sic ) 932 ATM_flx_sea = ATM_flux_int ( ATM_wbilo_sea ) 933 ATM_flx_ter = ATM_flux_int ( ATM_wbilo_ter ) 934 ATM_flx_lic = ATM_flux_int ( ATM_wbilo_lic ) 935 ATM_flx_calving = ATM_flux_int ( ATM_fqcalving ) 936 ATM_flx_qfonte = ATM_flux_int ( ATM_fqfonte ) 937 ATM_flx_precip = ATM_flux_int ( ATM_precip ) 938 ATM_flx_snowf = ATM_flux_int ( ATM_snowf ) 939 ATM_flx_evap = ATM_flux_int ( ATM_evap ) 940 ATM_flx_runlic = ATM_flux_int ( ATM_runofflic ) 941 942 RUN_flx_coastal = ONE_flux_int ( RUN_coastalflow) 943 RUN_flx_river = ONE_flux_int ( RUN_riverflow ) 944 RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 945 RUN_flx_river_cpl = ONE_flux_int ( RUN_riverflow_cpl ) 946 RUN_flx_drainage = SRF_flux_int ( RUN_drainage ) 947 RUN_flx_riversret = SRF_flux_int ( RUN_riversret ) 948 RUN_flx_runoff = SRF_flux_int ( RUN_runoff ) 949 RUN_flx_input = SRF_flux_int ( RUN_input ) 950 RUN_flx_output = ONE_flux_int ( RUN_output ) 951 952 ATM_flx_emp = ATM_flux_int (ATM_emp) 953 ATM_flx_wbilo = ATM_flux_int (ATM_wbilo) 954 702 955 RUN_flx_bil = RUN_flx_input - RUN_flx_output 703 704 956 RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 705 957 706 ATM_flx_bilo2 = ATM_flux_int (ATM_bilo) 707 708 709 echo (' wbilo_oce {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_oce , ATM_flx_oce / dtime_sec*1E-6/ATM_rho, ATM_flx_oce /ATM_aire_sea_tot/ATM_rho )) 710 echo (' wbilo_sic {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_sic , ATM_flx_sic / dtime_sec*1E-6/ATM_rho, ATM_flx_sic /ATM_aire_sea_tot/ATM_rho )) 711 echo (' wbilo_sic+oce {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_sea , ATM_flx_sea / dtime_sec*1E-6/ATM_rho, ATM_flx_sea /ATM_aire_sea_tot/ATM_rho )) 712 echo (' wbilo_ter {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_ter , ATM_flx_ter / dtime_sec*1E-6/ATM_rho, ATM_flx_ter /ATM_aire_sea_tot/ATM_rho )) 713 echo (' wbilo_lic {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_lic , ATM_flx_lic / dtime_sec*1E-6/ATM_rho, ATM_flx_lic /ATM_aire_sea_tot/ATM_rho )) 714 echo (' Sum wbilo_* {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_bilo , ATM_flx_bilo / dtime_sec*1E-6/ATM_rho, ATM_flx_bilo /ATM_aire_sea_tot/ATM_rho )) 715 echo (' E-P {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_emp , ATM_flx_emp / dtime_sec*1E-6/ATM_rho, ATM_flx_emp /ATM_aire_sea_tot/ATM_rho )) 716 echo (' calving {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_calving , ATM_flx_calving / dtime_sec*1E-6/ATM_rho, ATM_flx_calving /ATM_aire_sea_tot/ATM_rho )) 717 echo (' fqfonte {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_qfonte , ATM_flx_qfonte / dtime_sec*1E-6/ATM_rho, ATM_flx_qfonte /ATM_aire_sea_tot/ATM_rho )) 718 echo (' precip {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_precip , ATM_flx_precip / dtime_sec*1E-6/ATM_rho, ATM_flx_precip /ATM_aire_sea_tot/ATM_rho )) 719 echo (' snowf {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_snowf , ATM_flx_snowf / dtime_sec*1E-6/ATM_rho, ATM_flx_snowf /ATM_aire_sea_tot/ATM_rho )) 720 echo (' evap {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_evap , ATM_flx_evap / dtime_sec*1E-6/ATM_rho, ATM_flx_evap /ATM_aire_sea_tot/ATM_rho )) 721 echo (' coastalflow {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_coastal , RUN_flx_coastal / dtime_sec*1E-6/ATM_rho, RUN_flx_coastal /ATM_aire_sea_tot/ATM_rho )) 722 echo (' riverflow {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_river , RUN_flx_river / dtime_sec*1E-6/ATM_rho, RUN_flx_river /ATM_aire_sea_tot/ATM_rho )) 723 echo (' river+coastal {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_rivcoa , RUN_flx_rivcoa / dtime_sec*1E-6/ATM_rho, RUN_flx_rivcoa /ATM_aire_sea_tot/ATM_rho )) 724 echo (' drainage {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_drainage , RUN_flx_drainage / dtime_sec*1E-6/ATM_rho, RUN_flx_drainage /ATM_aire_sea_tot/ATM_rho )) 725 echo (' riversret {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_riversret, RUN_flx_riversret/ dtime_sec*1E-6/ATM_rho, RUN_flx_riversret/ATM_aire_sea_tot/ATM_rho )) 726 echo (' runoff {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_runoff , RUN_flx_runoff / dtime_sec*1E-6/ATM_rho, RUN_flx_runoff /ATM_aire_sea_tot/ATM_rho )) 727 echo (' river in {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_input , RUN_flx_input / dtime_sec*1E-6/ATM_rho, RUN_flx_input /ATM_aire_sea_tot/ATM_rho )) 728 echo (' river out {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_output , RUN_flx_output / dtime_sec*1E-6/ATM_rho, RUN_flx_output /ATM_aire_sea_tot/ATM_rho )) 729 echo (' river bil {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_bil , RUN_flx_bil / dtime_sec*1E-6/ATM_rho, RUN_flx_bil /ATM_aire_sea_tot/ATM_rho )) 730 echo (' runoff lic {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_runlic , ATM_flx_runlic / dtime_sec*1E-6/ATM_rho, ATM_flx_runlic /ATM_aire_sea_tot/ATM_rho )) 731 732 ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river 733 echo ('\n Global {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_budget , ATM_flx_budget / dtime_sec*1E-9, ATM_flx_budget /ATM_aire_sea_tot/ATM_rho )) 734 735 #echo (' E-P-R {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_emp , ATM_flx_emp / dtime_sec*1E-6/ATM_rho, ATM_flx_emp /ATM_aire_sea_tot/ATM_rho )) 736 737 738 echo ( '------------------------------------------------------------------------------------' ) 739 echo ( '-- SECHIBA fluxes' ) 740 741 742 SRF_flx_evap = ATM_flux_int ( SRF_evap ) 743 SRF_flx_snowf = ATM_flux_int ( SRF_snowf ) 744 SRF_flx_subli = ATM_flux_int ( SRF_subli ) 745 SRF_flx_transpir = ATM_flux_int ( SRF_transpir ) 746 747 echo (' evap {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_evap , SRF_flx_evap / dtime_sec*1E-6/ATM_rho, SRF_flx_evap /ATM_aire_sea_tot/ATM_rho )) 748 echo (' snowf {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_snowf , SRF_flx_snowf / dtime_sec*1E-6/ATM_rho, SRF_flx_snowf /ATM_aire_sea_tot/ATM_rho )) 749 echo (' subli {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_subli , SRF_flx_subli / dtime_sec*1E-6/ATM_rho, SRF_flx_subli /ATM_aire_sea_tot/ATM_rho )) 750 echo (' transpir {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_transpir , SRF_flx_transpir / dtime_sec*1E-6/ATM_rho, SRF_flx_transpir /ATM_aire_sea_tot/ATM_rho )) 958 prtFlux ('wbilo_oce ', ATM_flx_oce, 'f' ) 959 prtFlux ('wbilo_oce ', ATM_flx_oce, 'f' ) 960 prtFlux ('wbilo_sic ', ATM_flx_sic, 'f' ) 961 prtFlux ('wbilo_sic+oce ', ATM_flx_sea, 'f' ) 962 prtFlux ('wbilo_ter ', ATM_flx_ter, 'f' ) 963 prtFlux ('wbilo_lic ', ATM_flx_lic, 'f' ) 964 prtFlux ('Sum wbilo_* ', ATM_flx_wbilo , 'f', True) 965 prtFlux ('E-P ', ATM_flx_emp , 'f', True) 966 prtFlux ('calving ', ATM_flx_calving , 'f' ) 967 prtFlux ('fqfonte ', ATM_flx_qfonte , 'f' ) 968 prtFlux ('precip ', ATM_flx_precip , 'f' ) 969 prtFlux ('snowf ', ATM_flx_snowf , 'f' ) 970 prtFlux ('evap ', ATM_flx_evap , 'f' ) 971 prtFlux ('coastalflow ', RUN_flx_coastal , 'f' ) 972 prtFlux ('riverflow ', RUN_flx_river , 'f' ) 973 prtFlux ('coastal_cpl ', RUN_flx_coastal_cpl, 'f' ) 974 prtFlux ('riverf_cpl ', RUN_flx_river_cpl , 'f' ) 975 prtFlux ('river+coastal ', RUN_flx_rivcoa , 'f' ) 976 prtFlux ('drainage ', RUN_flx_drainage , 'f' ) 977 prtFlux ('riversret ', RUN_flx_riversret , 'f' ) 978 prtFlux ('runoff ', RUN_flx_runoff , 'f' ) 979 prtFlux ('river in ', RUN_flx_input , 'f' ) 980 prtFlux ('river out ', RUN_flx_output , 'f' ) 981 prtFlux ('river bil ', RUN_flx_bil , 'f' ) 982 prtFlux ('runoff lic ', ATM_flx_runlic , 'f' ) 983 984 ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river 985 echo ('') 986 echo (' Global {:12.3e} kg | {:12.4f} Sv | {:12.4f} m '.format ( ATM_flx_budget , ATM_flx_budget / dtime_sec*1E-9, ATM_flx_budget /ATM_aire_sea_tot/ATM_rho )) 987 988 #echo (' E-P-R {:12.3e} kg | {:12.4e} Sv | {:12.4f} m '.format ( ATM_flx_emp , ATM_flx_emp / dtime_sec*1E-6/ATM_rho, ATM_flx_emp /ATM_aire_sea_tot/ATM_rho )) 989 990 echo (' ') 991 echo ( '\n====================================================================================' ) 992 echo ( f'-- Atmosphere -- {Title} ' ) 993 echo ( 'Mass begin = {:12.6e} kg | Mass end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 994 echo ( 'dmass (atm) = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( dDYN_mas_wat , kg2Sv(dDYN_mas_wat) , kg2myear(dDYN_mas_wat )*1000. )) 995 echo ( 'Sum wbilo_* = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( ATM_flx_wbilo, kg2Sv(ATM_flx_wbilo), kg2myear(ATM_flx_wbilo )*1000. )) 996 echo ( 'E-P = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( ATM_flx_emp , kg2Sv(ATM_flx_emp ) , kg2myear(ATM_flx_emp )*1000. )) 997 echo ( ' ' ) 998 prtFlux ( 'Perte eau atm ', ATM_flx_wbilo - dDYN_mas_wat, 'e', True ) 999 echo ( 'Perte eau atm = {:12.3e} (rel) '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat ) ) 1000 1001 echo (' ') 1002 prtFlux ( 'Perte eau atm', ATM_flx_emp - dDYN_mas_wat , 'e', True ) 1003 echo ( 'Perte eau atm = {:12.3e} (rel) '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat ) ) 1004 echo (' ') 1005 1006 echo ( '\n====================================================================================' ) 1007 echo ( f'-- SECHIBA fluxes -- {Title} ' ) 1008 1009 SRF_flx_rain = SRF_flux_int ( SRF_rain ) 1010 SRF_flx_evap = SRF_flux_int ( SRF_evap ) 1011 SRF_flx_snowf = SRF_flux_int ( SRF_snowf ) 1012 SRF_flx_subli = SRF_flux_int ( SRF_subli ) 1013 SRF_flx_transpir = SRF_flux_int ( SRF_transpir ) 1014 SRF_flx_emp = SRF_flux_int ( SRF_emp ) 1015 1016 RUN_flx_torouting = RUN_flx_runoff + RUN_flx_drainage 1017 RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river 1018 1019 SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage 1020 1021 prtFlux ('rain ', SRF_flx_rain , 'f' ) 1022 prtFlux ('evap ', SRF_flx_evap , 'f' ) 1023 prtFlux ('snowf ', SRF_flx_snowf , 'f' ) 1024 prtFlux ('E-P ', SRF_flx_emp , 'f' ) 1025 prtFlux ('subli ', SRF_flx_subli , 'f' ) 1026 prtFlux ('transpir ', SRF_flx_transpir , 'f' ) 1027 prtFlux ('to routing ', RUN_flx_torouting, 'f' ) 1028 prtFlux ('budget ', SRF_flx_all , 'f', True ) 1029 1030 echo ( '\n------------------------------------------------------------------------------------' ) 1031 echo ( 'Water content in surface ' ) 1032 echo ( 'SRF_mas_wat_beg = {:12.6e} kg | SRF_mas_wat_end = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 1033 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True) 1034 echo ( 'dMass (water srf) = {:12.4e} (rel) '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 1035 1036 echo ('') 1037 echo ( '\n====================================================================================' ) 1038 echo ( f'-- RUNOFF fluxes -- {Title} ' ) 1039 RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal 1040 prtFlux ('runoff ', RUN_flx_runoff , 'f' ) 1041 prtFlux ('drainage ', RUN_flx_drainage , 'f' ) 1042 prtFlux ('run+drain ', RUN_flx_torouting , 'f' ) 1043 prtFlux ('river ', RUN_flx_river , 'f' ) 1044 prtFlux ('coastal ', RUN_flx_coastal , 'f' ) 1045 prtFlux ('riv+coa ', RUN_flx_fromrouting , 'f' ) 1046 prtFlux ('budget ', RUN_flx_all , 'f' , True) 1047 1048 echo ( '\n------------------------------------------------------------------------------------' ) 1049 echo ( f'Water content in routing -- {Title} ' ) 1050 echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) ) 1051 prtFlux ( 'dMass (routing) ', dRUN_mas_wat+dSRF_mas_lake, 'f', True) 1052 1053 echo ( '\n------------------------------------------------------------------------------------' ) 1054 echo ( f'Water content in routing -- {Title} ' ) 1055 echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) ) 1056 prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', True ) 1057 1058 print ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 1059 1060 # PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &a 1061 # 1150 " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" -
TOOLS/WATER_BUDGET/CPL_waterbudget.py
r6271 r6277 692 692 RUN_mas_wat_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] + d_RUN_end ['slow_reservoir'] + d_RUN_end ['stream_reservoir']) 693 693 694 if Routing == ' ORCHIDEE' :694 if Routing == 'SECHIBA' : 695 695 RUN_mas_wat_beg = ONE_stock_int ( d_SRF_beg['fastres'] + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \ 696 696 + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] ) -
TOOLS/WATER_BUDGET/OCE_waterbudget.py
r6271 r6277 21 21 ### 22 22 ## Import system modules 23 import sys, os, shutil, subprocess 23 import sys, os, shutil, subprocess, platform 24 24 import numpy as np 25 25 import configparser, re … … 29 29 config.optionxform = str # To keep capitals 30 30 31 config['Files'] = {}32 config['System'] = {}33 34 ##-- Some physical constants35 #-- Earth Radius36 Ra = 40.e6 / (2. * np.pi)37 #-- Gravity38 g = 9.8139 #-- Ice density (kg/m3) in LIM340 ICE_rho_ice = 917.041 #-- Snow density (kg/m3) in LIM342 ICE_rho_sno = 330.043 #-- Ocean water density (kg/m3) in NEMO44 OCE_rho_liq = 1026.45 #-- Water density in ice pounds in SI346 ICE_rho_pnd = 1000.47 48 31 ## Read experiment parameters 49 ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 32 ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ; TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 33 34 ## 35 ARCHIVE=None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 36 TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 37 file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ; file_ICE_his=None ; file_OCE_sca=None 38 file_restart_beg=None ; file_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None 39 file_RUN_beg=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None 50 40 51 41 # Arguments passed 52 42 print ( "Name of Python script:", sys.argv[0] ) 53 IniFile = 43 IniFile = sys.argv[1] 54 44 print ("Input file : ", IniFile ) 55 45 config.read (IniFile) 46 if 'full' in IniFile : FullIniFile = IniFile 47 else : FullIniFile = 'full_' + IniFile 56 48 57 49 def setBool (chars) : … … 75 67 return setNum 76 68 77 print ('[Experiment]') 78 for VarName in config['Experiment'].keys() : 79 locals()[VarName] = config['Experiment'][VarName] 80 locals()[VarName] = setBool (locals()[VarName]) 81 locals()[VarName] = setNum (locals()[VarName]) 82 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 69 def setNone (chars) : 70 '''Convert specific char string to None if possible''' 71 if type (chars) == str : 72 if chars in ['None', 'NONE', 'none'] : setNone = None 73 else : setNone = chars 74 else : setNone = chars 75 return setNone 76 77 ## Reading config 78 for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 79 if Section in config.keys() : 80 print ( f'[{Section}]' ) 81 for VarName in config[Section].keys() : 82 locals()[VarName] = config[Section][VarName] 83 locals()[VarName] = setBool (locals()[VarName]) 84 locals()[VarName] = setNum (locals()[VarName]) 85 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 86 87 if not 'Files' in config.keys() : config['Files'] = {} 88 89 def unDefined (char) : 90 if char in globals () : 91 if char == None : return True 92 else : return False 93 else : return True 94 95 ##-- Some physical constants 96 #-- Earth Radius 97 if not 'Ra' in locals () : Ra = 6366197.7236758135 98 #-- Gravity 99 if not 'Grav' in locals () : Grav = 9.81 100 #-- Ice volumic mass (kg/m3) in LIM3 101 if not 'ICE_rho_ice' in locals () : ICE_rho_ice = 917.0 102 #-- Snow volumic mass (kg/m3) in LIM3 103 if not 'ICE_rho_sno' in locals () : ICE_rho_sno = 330.0 104 #-- Water density in ice pounds in SI3 105 if not 'ICE_rho_pnd' in locals () : ICE_rho_pnd = 1000. 106 #-- Ocean water volumic mass (kg/m3) in NEMO 107 if not 'OCE_rho_liq' in locals () : OCE_rho_liq = 1026. 108 #-- Water volumic mass in atmosphere 109 if not 'ATM_rho' in locals () : ATM_rho = 1000. 110 #-- Water volumic mass in surface reservoirs 111 if not 'SRF_rho' in locals () : SRF_rho = 1000. 112 #-- Water volumic mass of rivers 113 if not 'RUN_rho' in locals () : RUN_rho = 1000. 114 #-- Year length 115 if not 'YearLength' in locals () : YearLength = 365.25 * 24. * 60. * 60. 116 117 config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ICE_rho_pnd':ICE_rho_pnd, 118 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 119 120 83 121 84 122 # Where do we run ? … … 109 147 if IDRIS : 110 148 raise Exception ("Pour IDRIS : repertoires et chemins a definir") 111 149 150 config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU, 151 'Program' : "Generated by : " + str(sys.argv), 152 'HOSTNAME' : platform.node (), 'LOGNAME' : os.getlogin (), 153 'Python' : f'{platform.python_version ()}', 154 'OS' : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}', 155 'SVN_Author' : "$Author$", 156 'SVN_Date' : "$Date$", 157 'SVN_Revision' : "$Revision$", 158 'SVN_Id' : "$Id$", 159 'SVN_HeadURL' : "$HeadURL$"} 160 161 if libIGCM : 162 config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 112 163 113 164 ## Import specific module … … 117 168 118 169 # Output file 119 FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 170 if FileOut == None : 171 FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 172 config['Files']['FileOut'] = FileOut 173 120 174 f_out = open ( FileOut, mode = 'w' ) 121 175 122 176 # Function to print to stdout *and* output file 123 def echo (string) : 124 print ( string ) 125 f_out.write ( string + '\n' ) 177 def echo (string, end='\n') : 178 print ( str(string), end=end ) 179 sys.stdout.flush () 180 f_out.write ( str(string) + end ) 126 181 f_out.flush () 127 182 return None … … 146 201 if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 147 202 203 echo (' ') 204 echo ( f'JobName : {JobName}' ) 205 echo (Comment) 148 206 echo ( f'Working in TMPDIR : {TmpDir}' ) 149 207 … … 153 211 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 154 212 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 155 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 156 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 157 213 if dir_OCE_his == None : 214 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 215 config['Files']['dir_OCE_his'] = dir_OCE_his 216 if dir_ICE_his == None : 217 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 218 config['Files']['dir_OCE_his'] = dir_OCE_his 219 158 220 echo ( f'The analysis relies on files from the following model output directories : ' ) 159 221 echo ( f'{dir_OCE_his}' ) … … 161 223 162 224 #-- Files Names 163 if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 164 if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 165 225 if Period == None : 226 if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 227 if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 228 config['Files']['Period'] = Period 229 if FileCommon == None : 230 FileCommon = f'{JobName}_{Period}' 231 config['Files']['FileCommon'] = FileCommon 232 233 if Title == None : 234 Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 235 config['Files']['Title'] = Title 236 237 166 238 echo ('\nOpen history files' ) 167 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 168 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 169 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 239 if file_OCE_his == None : 240 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 241 file_OCE_his = file_OCE_his 242 if file_OCE_sca == None : 243 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 244 config['Files']['file_OCE_sca'] = file_OCE_sca 245 if file_ICE_his == None : 246 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 247 config['Files']['file_ICE_his'] = file_ICE_his 170 248 171 249 d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() … … 191 269 dtime_per_sec.attrs['unit'] = 's' 192 270 271 ## Number of years 272 NbYear = dtime_sec / YearLength 193 273 #-- Open restart files 194 YearRes = YearBegin - 1 # Year of the restart of beginning of simulation274 YearRes = YearBegin - 1 # Year of the restart of beginning of simulation 195 275 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 196 276 197 277 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 198 278 199 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' ) 200 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' ) 279 if TarRestartPeriod_beg == None : 280 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 281 TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 282 config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 283 284 if TarRestartPeriod_end == None : 285 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 286 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 287 TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231' 288 config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 289 290 if file_restart_beg == None : 291 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 292 config['Files']['file_restart_beg'] = file_restart_beg 293 if file_restart_end == None : 294 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 295 config['Files']['file_restart_end'] = file_restart_end 201 296 202 297 echo ( f'{file_restart_beg}' ) 203 298 echo ( f'{file_restart_end}' ) 204 299 205 file_OCE_beg = f'OCE_{JobName}_{YearRes}1231_restart.nc' 206 file_OCE_end = f'OCE_{JobName}_{YearEnd}1231_restart.nc' 207 file_ICE_beg = f'ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 208 file_ICE_end = f'ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 300 if file_OCE_beg == None : 301 file_OCE_beg = f'{TmpDir}/OCE_{JobName}_{YearRes}1231_restart.nc' 302 config['Files']['file_OCE_beg'] = file_OCE_beg 303 if file_OCE_end == None : 304 file_OCE_end = f'{TmpDir}/OCE_{JobName}_{YearEnd}1231_restart.nc' 305 config['Files']['file_OCE_end'] = file_OCE_end 306 if file_ICE_beg == None : 307 file_ICE_beg = f'{TmpDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 308 config['Files']['file_ICE_beg'] = file_ICE_beg 309 if file_ICE_end == None : 310 file_ICE_end = f'{TmpDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 311 config['Files']['file_ICE_end'] = file_ICE_end 209 312 210 313 echo ( f'{file_OCE_beg}' ) … … 221 324 return int (ndomain_opa) 222 325 223 if not os.path.exists ( os.path.join (TmpDir, file_OCE_beg) ) : 224 echo ( 'Extracting '+ os.path.join (TmpDir, file_OCE_beg) ) 326 if not os.path.exists ( file_OCE_beg) : 327 echo ( f'Extracting {file_OCE_beg}' ) 328 base_file_OCE_beg = os.path.basename (file_OCE_beg) 225 329 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart_0000.nc') ) : 226 330 command = f'cd {TmpDirOCE} ; tar xf {file_restart_beg} OCE_{JobName}_{YearRes}1231_restart_*.nc' … … 229 333 echo ('extract ndomain' ) 230 334 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') ) 231 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv { file_OCE_beg} {TmpDir}'335 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {base_file_OCE_beg} {TmpDir}' 232 336 echo ( command ) 233 337 os.system ( command ) 234 338 echo ( f'Rebuild done : {file_OCE_beg}') 235 339 236 if not os.path.exists ( os.path.join (TmpDir, file_OCE_end) ) : 340 if not os.path.exists ( file_OCE_end) : 341 echo ( f'Extracting {file_OCE_end}' ) 342 base_file_OCE_end = os.path.basename (file_OCE_end) 237 343 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart_0000.nc') ): 238 344 command = f'cd {TmpDirOCE} ; tar xf {file_restart_end} OCE_{JobName}_{YearEnd}1231_restart_*.nc' … … 241 347 echo ('extract ndomain' ) 242 348 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart') ) 243 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv { file_OCE_end} {TmpDir}'349 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {base_file_OCE_end} {TmpDir}' 244 350 echo ( command ) 245 351 os.system ( command ) 246 352 echo ( f'Rebuild done : {file_OCE_end}') 247 353 248 if not os.path.exists ( os.path.join (TmpDir, file_ICE_beg) ) : 354 if not os.path.exists ( file_ICE_beg) : 355 echo ( f'Extracting {file_ICE_beg}' ) 356 base_file_ICE_beg = os.path.basename (file_ICE_beg) 249 357 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod_0000.nc') ): 250 358 command = f'cd {TmpDirICE} ; tar xf {file_restart_beg} ICE_{JobName}_{YearRes}1231_restart_icemod_*.nc' … … 253 361 echo ('extract ndomain' ) 254 362 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 255 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv { file_ICE_beg} {TmpDir} '363 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_beg} {TmpDir} ' 256 364 echo ( command ) 257 365 os.system ( command ) 258 366 echo ( f'Rebuild done : {file_OCE_beg}') 259 367 260 if not os.path.exists ( os.path.join (TmpDir, file_ICE_end) ) : 368 if not os.path.exists ( file_ICE_end ) : 369 echo ( f'Extracting {file_ICE_end}' ) 370 base_file_ICE_end = os.path.basename (file_ICE_end) 261 371 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearEnd}1231_restart_icemod_0000.nc') ): 262 372 command = f'cd {TmpDirICE} ; tar xf {file_restart_end} ICE_{JobName}_{YearEnd}1231_restart_icemod_*.nc' … … 265 375 echo ('extract ndomain' ) 266 376 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 267 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv { file_ICE_end} {TmpDir}'377 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_end} {TmpDir}' 268 378 echo ( command ) 269 379 os.system ( command ) … … 282 392 d_ICE_end = xr.open_dataset ( os.path.join (TmpDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 283 393 394 ## 395 config_out = open (FullIniFile, 'w') 396 config.write (config_out ) 397 config_out.close () 398 399 400 def kg2Sv (val, rho=OCE_rho_liq) : 401 '''From kg to Sverdrup''' 402 return val/dtime_sec*1.0e-6/rho 403 404 def kg2myear (val, rho=OCE_rho_liq) : 405 '''From kg to m/year''' 406 return val/OCE_aire_tot/rho/NbYear 407 408 def var2prt (var, small=False) : 409 if small : return var , kg2Sv(var)*1000., kg2myear(var)*1000. 410 else : return var , kg2Sv(var) , kg2myear(var) 411 412 def prtFlux (Desc, var, Form='F', small=False) : 413 if small : 414 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 415 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 416 else : 417 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 418 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 419 echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) ) 420 return None 421 284 422 # Get mask and surfaces 285 423 sos = d_OCE_his ['sos'][0].squeeze() … … 300 438 def ONE_stock_int (stock) : 301 439 '''Sum stock''' 302 ONE_stock_int = np.sum ( np.sort ( (stock * OCE_aire).to_masked_array().ravel()) )440 ONE_stock_int = np.sum ( np.sort ( (stock ).to_masked_array().ravel()) ) 303 441 return ONE_stock_int 304 442 … … 306 444 '''Integrate flux on oce grid''' 307 445 OCE_flux_int = np.sum ( np.sort ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel()) ) 446 return OCE_flux_int 447 448 def ONE_flux_int (flux) : 449 '''Integrate flux on oce grid''' 450 OCE_flux_int = np.sum ( np.sort ( (flux * dtime_per_sec).to_masked_array().ravel()) ) 308 451 return OCE_flux_int 309 452 … … 343 486 echo ( 'dOCE vol = {:12.3e} m^3'.format ( dOCE_vol_wat) ) 344 487 echo ( 'dOCE ssh = {:12.3e} m '.format ( dOCE_vol_wat/OCE_aire_tot) ) 345 echo ( 'dOCE mass = {:12.3e} kg '.format ( dOCE_mas_wat) ) 346 echo ( 'dOCE mass = {:12.3e} Sv '.format ( dOCE_mas_wat/dtime_sec*1E-9) ) 488 prtFlux ( 'dOCE mass ', dOCE_mas_wat, 'e' ) 347 489 348 490 if NEMO == 3.6 : 349 491 echo ( 'dOCE e3tn vol = {:12.3e} m^3'.format ( dOCE_e3tn_vol) ) 350 echo ( 'dOCE e3tn mass = {:12.3e} kg '.format ( dOCE_e3tn_mas))492 prtFlux ( 'dOCE e3tn mass', dOCE_e3tn_mas, 'e' ) 351 493 352 494 ## Glace et neige 353 if NEMO == 3.6 : 495 if NEMO == 3.6 : 354 496 ICE_ice_beg = d_ICE_beg['v_i_htc1']+d_ICE_beg['v_i_htc2']+d_ICE_beg['v_i_htc3']+d_ICE_beg['v_i_htc4']+d_ICE_beg['v_i_htc5'] 355 497 ICE_ice_end = d_ICE_end['v_i_htc1']+d_ICE_end['v_i_htc2']+d_ICE_end['v_i_htc3']+d_ICE_end['v_i_htc4']+d_ICE_end['v_i_htc5'] … … 358 500 ICE_sno_end = d_ICE_end['v_s_htc1']+d_ICE_end['v_s_htc2']+d_ICE_end['v_s_htc3']+d_ICE_end['v_s_htc4']+d_ICE_end['v_s_htc5'] 359 501 360 ICE_pnd_beg = 0.0 361 ICE_pnd_end = 0.0 362 363 ICE_fzl_beg = 0.0 364 ICE_fzl_end = 0.0 502 ICE_pnd_beg = 0.0 ; ICE_pnd_end = 0.0 ; ICE_fzl_beg = 0.0 ; ICE_fzl_end = 0.0 365 503 366 504 ICE_mas_wat_beg = ICE_ice_beg*ICE_rho_ice + ICE_sno_beg*ICE_rho_sno … … 424 562 echo ( '\n------------------------------------------------------------') 425 563 echo ( 'Variation du contenu en eau ocean + glace ' ) 426 echo ( 'dMass (ocean) = {:12.6e} kg '.format(dSEA_mas_wat) ) 427 echo ( 'dVol (ocean) = {:12.3e} Sv '.format(dSEA_mas_wat/dtime_sec*1E-9) ) 428 echo ( 'dVol (ocean) = {:12.3e} m '.format(dSEA_mas_wat*1E-3/OCE_aire_tot) ) 564 prtFlux ( 'dMass (ocean)', dSEA_mas_wat, 'e', True ) 429 565 430 566 … … 504 640 505 641 echo ('\n-- Fields:' ) 506 echo ('OCE+ICE budget {:12.3e} (kg) {:12.3e} (Sv) '.format ( OCE_mas_all , OCE_mas_all / dtime_sec*1E-9 ))507 echo (' EMPMR {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_empmr , OCE_mas_empmr / dtime_sec*1E-9 ))508 echo (' WFOB {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfob , OCE_mas_wfob / dtime_sec*1E-9 ))509 echo (' EMP_OCE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_oce , OCE_mas_emp_oce / dtime_sec*1E-9 ))510 echo (' EMP_ICE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_ice , OCE_mas_emp_ice / dtime_sec*1E-9 ))511 echo (' EMP {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp , OCE_mas_emp / dtime_sec*1E-9 ))512 echo (' ICEBERG {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceberg , OCE_mas_iceberg / dtime_sec*1E-9 ))513 echo (' ICESHELF {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceshelf , OCE_mas_iceshelf / dtime_sec*1E-9 ))514 echo (' CALVING {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calving , OCE_mas_calving / dtime_sec*1E-9 ))515 echo (' FRIVER {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_friver , OCE_mas_friver / dtime_sec*1E-9 ))516 echo (' RUNOFFS {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_runoffs , OCE_mas_runoffs / dtime_sec*1E-9 ))517 echo (' WFXICE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxice , OCE_mas_wfxice / dtime_sec*1E-9 ))518 echo (' WFXSNW {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsnw , OCE_mas_wfxsnw / dtime_sec*1E-9 ))519 echo (' WFXSUB {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsub , OCE_mas_wfxsub / dtime_sec*1E-9 ))520 echo (' WFXPND {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxpnd , ICE_mas_wfxpnd / dtime_sec*1E-9 ))521 echo (' WFXSNW_SUB {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_sub, ICE_mas_wfxsnw_sub / dtime_sec*1E-9 ))522 echo (' WFXSNW_PRE {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_pre, ICE_mas_wfxsnw_pre / dtime_sec*1E-9 ))523 echo (' WFXSUB_ERR {:12.3e} (kg) {:12.3e} (Sv) '.format ( ICE_mas_wfxsub_err, ICE_mas_wfxsub_err / dtime_sec*1E-9 ))524 echo (' EVAP_OCE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_evap_oce , OCE_mas_evap_oce / dtime_sec*1E-9 ))525 echo (' EVAP_ICE {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_evap_ice , ICE_mas_evap_ice / dtime_sec*1E-9 ))526 echo (' SNOW_OCE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_oce , OCE_mas_snow_oce / dtime_sec*1E-9 ))527 echo (' SNOW_ICE {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_ice , OCE_mas_snow_ice / dtime_sec*1E-9 ))528 echo (' RAIN {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_rain , OCE_mas_rain / dtime_sec*1E-9 ))529 echo (' FWB {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_fwb , OCE_mas_fwb / dtime_sec*1E-9 ))530 echo (' SSR {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_ssr , OCE_mas_ssr / dtime_sec*1E-9 ))531 echo (' WFCORR {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfcorr , OCE_mas_wfcorr / dtime_sec*1E-9 ))532 echo (' BERG_ICB {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_berg_icb , OCE_mas_berg_icb / dtime_sec*1E-9 ))533 echo (' CALV_ICB {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calv_icb , OCE_mas_calv_icb / dtime_sec*1E-9 ))642 prtFlux ('OCE+ICE budget ', OCE_mas_all , 'e', True) 643 prtFlux (' EMPMR ', OCE_mas_empmr , 'e', True) 644 prtFlux (' WFOB ', OCE_mas_wfob , 'e', True) 645 prtFlux (' EMP_OCE ', OCE_mas_emp_oce , 'e', True) 646 prtFlux (' EMP_ICE ', OCE_mas_emp_ice , 'e', True) 647 prtFlux (' EMP ', OCE_mas_emp , 'e', True) 648 prtFlux (' ICEBERG ', OCE_mas_iceberg , 'e', True) 649 prtFlux (' ICESHELF ', OCE_mas_iceshelf , 'e', True) 650 prtFlux (' CALVING ', OCE_mas_calving , 'e', True) 651 prtFlux (' FRIVER ', OCE_mas_friver , 'e', True) 652 prtFlux (' RUNOFFS ', OCE_mas_runoffs , 'e', True) 653 prtFlux (' WFXICE ', OCE_mas_wfxice , 'e', True) 654 prtFlux (' WFXSNW ', OCE_mas_wfxsnw , 'e', True) 655 prtFlux (' WFXSUB ', OCE_mas_wfxsub , 'e', True) 656 prtFlux (' WFXPND ', ICE_mas_wfxpnd , 'e', True) 657 prtFlux (' WFXSNW_SUB ', ICE_mas_wfxsnw_sub, 'e', True) 658 prtFlux (' WFXSNW_PRE ', ICE_mas_wfxsnw_pre, 'e', True) 659 prtFlux (' WFXSUB_ERR ', ICE_mas_wfxsub_err, 'e', True) 660 prtFlux (' EVAP_OCE ', OCE_mas_evap_oce , 'e' ) 661 prtFlux (' EVAP_ICE ', ICE_mas_evap_ice , 'e', True) 662 prtFlux (' SNOW_OCE ', OCE_mas_snow_oce , 'e', True) 663 prtFlux (' SNOW_ICE ', OCE_mas_snow_ice , 'e', True) 664 prtFlux (' RAIN ', OCE_mas_rain , 'e' ) 665 prtFlux (' FWB ', OCE_mas_fwb , 'e', True) 666 prtFlux (' SSR ', OCE_mas_ssr , 'e', True) 667 prtFlux (' WFCORR ', OCE_mas_wfcorr , 'e', True) 668 prtFlux (' BERG_ICB ', OCE_mas_berg_icb , 'e', True) 669 prtFlux (' CALV_ICB ', OCE_mas_calv_icb , 'e', True) 534 670 535 671 echo ('\n===> Comparing ===>' )
Note: See TracChangeset
for help on using the changeset viewer.