- Timestamp:
- 06/13/23 12:58:38 (11 months ago)
- Location:
- TOOLS/WATER_BUDGET
- Files:
-
- 2 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/WATER_BUDGET/ATM_waterbudget.py
r6277 r6508 25 25 import configparser, re 26 26 27 ## Creates parser 28 config = configparser.ConfigParser() 27 # Check python version 28 if sys.version_info < (3, 8, 0) : 29 print ( f'Python version : {platform.python_version()}' ) 30 raise Exception ( "Minimum Python version is 3.8" ) 31 32 ## Import local module 33 import WaterUtils as wu 34 35 ## Creates parser for reading .ini input file 36 config = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() ) 29 37 config.optionxform = str # To keep capitals 30 38 31 ## Read experiment parameters32 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=None39 ## Experiment parameters 40 ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; 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 41 34 42 ## 35 ARCHIVE =None ; STORAGE =None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None36 TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None43 ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 44 TmpDir=None ; RunDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 37 45 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 44 45 # Arguments passed 46 tar_restart_beg=None ; tar_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 47 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 ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None ; ContinueOnError=False ; ErrorCount=0 48 tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None ; tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None 49 tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None ; tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None 50 SortIco = False 51 52 # Read command line arguments 46 53 print ( "Name of Python script:", sys.argv[0] ) 47 54 IniFile = sys.argv[1] 55 # Text existence of IniFile 56 if not os.path.exists (IniFile ) : 57 raise FileExistsError ( f"File not found : {IniFile = }" ) 58 59 if 'full' in IniFile : FullIniFile = IniFile 60 else : FullIniFile = 'full_' + IniFile 61 48 62 print ("Input file : ", IniFile ) 49 63 config.read (IniFile) 50 64 FullIniFile = 'full_' + IniFile 51 65 52 def setBool (chars) : 53 '''Convert specific char string in boolean if possible''' 54 setBool = chars 55 for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : 56 if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] 57 return setBool 58 59 def setNum (chars) : 60 '''Convert specific char string in integer or real if possible''' 61 if type (chars) == str : 62 realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") 63 isReal = realnum.match(chars.strip()) != None 64 isInt = chars.strip().isdigit() 65 if isReal : 66 if isInt : setNum = int (chars) 67 else : setNum = float (chars) 68 else : setNum = chars 69 else : setNum = chars 70 return setNum 71 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}]' ) 66 ## Reading config file 67 for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] : 68 if Section in config.keys () : 69 print ( f'\nReading [{Section}]' ) 84 70 for VarName in config[Section].keys() : 85 71 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 72 exec ( f'{VarName} = wu.setBool ({VarName})' ) 73 exec ( f'{VarName} = wu.setNum ({VarName})' ) 74 exec ( f'{VarName} = wu.setNone ({VarName})' ) 75 exec ( f'wu.{VarName} = {VarName}' ) 76 print ( f'{VarName:25} set to : {locals()[VarName]:}' ) 77 #exec ( f'del {VarName}' ) 78 79 print ( f'\nConfig file readed : {IniFile} ' ) 80 98 81 ##-- 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. 82 if wu.unDefined ( 'Ra' ) : Ra = 6366197.7236758135 #-- Earth Radius (m) 83 if wu.unDefined ( 'Grav' ) : Grav = 9.81 #-- Gravity (m^2/s 84 if wu.unDefined ( 'ICE_rho_ice' ) : ICE_rho_ice = 917.0 #-- Ice volumic mass (kg/m3) in LIM3 85 if wu.unDefined ( 'ICE_rho_sno') : ICE_rho_sno = 330.0 #-- Snow volumic mass (kg/m3) in LIM3 86 if wu.unDefined ( 'OCE_rho_liq' ) : OCE_rho_liq = 1026.0 #-- Ocean water volumic mass (kg/m3) in NEMO 87 if wu.unDefined ( 'ATM_rho' ) : ATM_rho = 1000.0 #-- Water volumic mass in atmosphere (kg/m^3) 88 if wu.unDefined ( 'SRF_rho' ) : SRF_rho = 1000.0 #-- Water volumic mass in surface reservoir (kg/m^3) 89 if wu.unDefined ( 'RUN_rho' ) : RUN_rho = 1000.0 #-- Water volumic mass of rivers (kg/m^3) 90 if wu.unDefined ( 'ICE_rho_pnd' ) : ICE_rho_pnd = 1000. #-- Water volumic mass in ice ponds (kg/m^3) 91 if wu.unDefined ( 'YearLength' ) : YearLength = 365.25 * 24. * 60. * 60. #-- Year length (s) - Use only to compute drif in approximate unit 92 93 ##-- Set libIGCM and machine dependant values 94 if not 'Files' in config.keys () : config['Files'] = {} 117 95 118 96 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 97 120 ## -- 121 ICO = ( 'ICO' in ATM ) 122 LMDZ = ( 'LMD' in ATM ) 123 124 125 # Where do we run ? 126 SysName, NodeName, Release, Version, Machine = os.uname () 127 TGCC = ( 'irene' in NodeName ) 128 IDRIS = ( 'jeanzay' in NodeName ) 129 130 ## Set site specific libIGCM directories, and other specific stuff 131 if TGCC : 132 CPU = subprocess.getoutput ( 'lscpu | grep "Model name"' ) 133 if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene' 134 if "AMD" in CPU : Machine = 'irene-amd' 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' ) 142 143 ## Specific to run at TGCC. 144 # Needed before importing a NetCDF library (netCDF4, xarray, cmds, etc ...) 145 import mpi4py 146 mpi4py.rc.initialize = False 147 148 ## Creates output directory 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 153 if IDRIS : 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$"} 98 config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco} 99 100 ## -------------------------- 101 ICO = ( 'ICO' in wu.ATM ) 102 LMDZ = ( 'LMD' in wu.ATM ) 103 104 with open ('SetLibIGCM.py') as f: exec ( f.read() ) 105 config['Files']['TmpDir'] = TmpDir 166 106 167 107 if libIGCM : 168 108 config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 169 109 170 ## Import specific module110 ##-- Import specific module 171 111 import nemo, lmdz 172 ## Now import needed scientific modules112 ##-- Now import needed scientific modules 173 113 import xarray as xr 174 114 175 # Output file 176 if FileOut == None : 177 FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 178 config['Files']['FileOut'] = FileOut 115 ##- Output file with water budget diagnostics 116 if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 117 config['Files']['FileOut'] = FileOut 179 118 180 119 f_out = open ( FileOut, mode = 'w' ) 181 120 182 # Function to print to stdout *and* output file 121 ##- Useful functions 122 def kg2Sv (val, rho=ATM_rho) : 123 '''From kg to Sverdrup''' 124 return val/dtime_sec*1.0e-6/rho 125 126 def kg2myear (val, rho=ATM_rho) : 127 '''From kg to m/year''' 128 return val/ATM_aire_sea_tot/rho/NbYear 129 130 def var2prt (var, small=False, rho=ATM_rho) : 131 if small : return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000. 132 else : return var , kg2Sv(var, rho=rho) , kg2myear(var, rho=rho) 133 134 def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) : 135 if small : 136 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 137 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 138 else : 139 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 140 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 141 echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) ) 142 return None 143 183 144 def echo (string, end='\n') : 145 '''Function to print to stdout *and* output file''' 184 146 print ( str(string), end=end ) 185 147 sys.stdout.flush () … … 188 150 return None 189 151 190 ## Set libIGCM directories 191 R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT') 192 R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT') 193 194 L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 195 R_SAVE = os.path.join ( R_OUT, L_EXP ) 196 R_BUFR = os.path.join ( R_BUF, L_EXP ) 197 POST_DIR = os.path.join ( R_BUF, L_EXP, 'Out' ) 198 REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' ) 199 R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 200 R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 201 202 #if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir ) 203 if not os.path.isdir (TmpDir) : os.mkdir (TmpDir) 204 TmpDirOCE = os.path.join (TmpDir, 'OCE') 205 TmpDirICE = os.path.join (TmpDir, 'ICE') 206 if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE ) 207 if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 152 ##- Set libIGCM directories 153 if wu.unDefined ('R_OUT' ) : R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT' ) 154 if wu.unDefined ('R_BUF' ) : R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT' ) 155 if wu.unDefined ('L_EXP' ) : L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 156 if wu.unDefined ('R_SAVE' ) : R_SAVE = os.path.join ( R_OUT, L_EXP ) 157 if wu.unDefined ('R_BUFR' ) : R_BUFR = os.path.join ( R_BUF, L_EXP ) 158 if wu.unDefined ('POST_DIR' ) : POST_DIR = os.path.join ( R_BUFR, 'Out' ) 159 if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' ) 160 if wu.unDefined ('R_BUF_KSH' ) : R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 161 if wu.unDefined ('R_FIGR' ) : R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 162 163 config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR } 164 165 ##- Set directory to extract files 166 if wu.unDefined ( 'RunDir' ) : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 167 config['Files']['RunDir'] = RunDir 168 169 if not os.path.isdir ( RunDir ) : os.makedirs ( RunDir ) 170 171 ##- Set directories to rebuild ocean and ice restart files 172 if wu.unDefined ( 'RunDirOCE' ) : RunDirOCE = os.path.join ( RunDir, 'OCE' ) 173 if wu.unDefined ( 'RunDirICE' ) : RunDirICE = os.path.join ( RunDir, 'ICE' ) 174 if not os.path.exists ( RunDirOCE ) : os.mkdir ( RunDirOCE ) 175 if not os.path.exists ( RunDirICE ) : os.mkdir ( RunDirICE ) 208 176 209 177 echo (' ') 210 echo ( f'JobName : {JobName}' ) 211 echo (Comment) 212 echo ( f'Working in TMPDIR : {TmpDir}' ) 213 214 echo ( f'\nDealing with {L_EXP}' ) 215 216 #-- Model output directories 178 echo ( f'JobName : {JobName}' ) 179 echo ( f'Comment : {Comment}' ) 180 echo ( f'TmpDir : {TmpDir}' ) 181 echo ( f'RunDir : {RunDir}' ) 182 echo ( f'RunDirOCE : {RunDirOCE}' ) 183 echo ( f'RunDirICE : {RunDirICE}' ) 184 185 echo ( f'\nDealing with {L_EXP}' ) 186 187 ##-- Creates model output directory names 217 188 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 218 189 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 219 if dir_ATM_his == None:190 if wu.unDefined ('dir_ATM_his' ) : 220 191 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 221 192 config['Files']['dir_ATM_his'] = dir_ATM_his 222 if dir_SRF_his == None:193 if wu.unDefined ( 'dir_SRF_his' ) : 223 194 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 224 195 config['Files']['dir_SRF_his'] = dir_SRF_his 225 196 226 197 echo ( f'The analysis relies on files from the following model output directories : ' ) 227 echo ( f'{dir_ATM_his }' )228 echo ( f'{dir_SRF_his }' )229 230 # -- Files Names231 if Period == None:198 echo ( f'{dir_ATM_his = }' ) 199 echo ( f'{dir_SRF_his = }' ) 200 201 ##-- Creates files names 202 if wu.unDefined ( 'Period' ) : 232 203 if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 233 204 if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 234 205 config['Files']['Period'] = Period 235 if FileCommon == None:206 if wu.unDefined ( 'FileCommon' ) : 236 207 FileCommon = f'{JobName}_{Period}' 237 208 config['Files']['FileCommon'] = FileCommon 238 209 239 if Title == None:210 if wu.unDefined ( 'Title' ) : 240 211 Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 241 212 config['Files']['Title'] = Title 242 213 243 214 echo ('\nOpen history files' ) 244 if file_ATM_his == None:215 if wu.unDefined ( 'file_ATM_his' ) : 245 216 file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 246 217 config['Files']['file_ATM_his'] = file_ATM_his 247 if file_SRF_his == None:218 if wu.unDefined ( 'file_SRF_his' ) : 248 219 file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 249 220 config['Files']['file_SRF_his'] = file_SRF_his … … 257 228 d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 258 229 d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 259 if Routing == 'SECHIBA' : 260 d_RUN_his = d_SRF_his 261 if Routing == 'SIMPLE' : 262 d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 263 264 echo ( file_ATM_his ) 265 echo ( file_SRF_his ) 266 if Routing == 'SIMPLE' : 267 echo ( file_RUN_his ) 268 269 ## Compute run length 230 if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 231 if Routing == 'SIMPLE' : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 232 233 echo ( f'{file_ATM_his = }' ) 234 echo ( f'{file_SRF_his = }' ) 235 if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' ) 236 237 ##-- Compute run length 270 238 dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() ) 271 239 echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) ) 272 240 dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds 273 241 274 ## Compute length of each period242 ##-- Compute length of each period 275 243 dtime_per = (d_ATM_his.time_counter_bounds[:,-1] - d_ATM_his.time_counter_bounds[:,0] ) 276 244 echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) ) 277 245 dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds 278 246 dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 279 280 ## Number of years 247 dtime_per_sec.attrs['unit'] = 's' 248 249 ##-- Number of years (approximative) 281 250 NbYear = dtime_sec / YearLength 282 251 283 #-- Open restart files 284 252 ##-- Extract restart files from tar 285 253 YearRes = YearBegin - 1 # Year of the restart of beginning of simulation 286 254 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 287 255 288 if TarRestartPeriod_beg == None : 256 config['Files']['YearPre'] = f'{YearBegin}' 257 258 if wu.unDefined ( 'TarRestartPeriod_beg' ) : 289 259 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 290 260 TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 291 261 config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 292 262 293 if TarRestartPeriod_end == None:263 if wu.unDefined ( 'TarRestartPeriod_end ' ) : 294 264 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 295 265 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') … … 297 267 config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 298 268 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 306 echo ( f'{file_restart_beg}' ) 307 echo ( f'{file_restart_end}' ) 308 309 if file_ATM_beg == None : 310 file_ATM_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc' 269 if wu.unDefined ( 'tar_restart_beg' ) : 270 tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 271 config['Files']['tar_restart_beg'] = tar_restart_beg 272 if wu.unDefined ( 'tar_restart_end' ) : 273 tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 274 config['Files']['tar_restart_end'] = tar_restart_end 275 276 echo ( f'{tar_restart_beg = }' ) 277 echo ( f'{tar_restart_end = }' ) 278 279 ##-- Names of tar files with restarts 280 if wu.unDefined ( 'SRF_HIS' ) : SRF_HIS = ATM_HIS 281 282 if wu.unDefined ( 'tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg 283 if wu.unDefined ( 'tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg 284 if wu.unDefined ( 'tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg 285 if wu.unDefined ( 'tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg 286 if wu.unDefined ( 'tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg 287 if wu.unDefined ( 'tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg 288 289 if wu.unDefined ( 'tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end 290 if wu.unDefined ( 'tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end 291 if wu.unDefined ( 'tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end 292 if wu.unDefined ( 'tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end 293 if wu.unDefined ( 'tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end 294 if wu.unDefined ( 'tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end 295 296 if wu.unDefined ( 'file_ATM_beg' ) : 297 file_ATM_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc' 311 298 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'299 if wu.unDefined ( 'file_ATM_end' ) : 300 file_ATM_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc' 314 301 config['Files']['file_ATM_end'] = file_ATM_end 315 302 316 303 liste_beg = [file_ATM_beg, ] 317 304 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' 305 306 if wu.unDefined ( 'file_DYN_beg' ) : 307 if LMDZ : file_DYN_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restart.nc' 308 if ICO : file_DYN_beg = f'{RunDir}/ICO_{JobName}_{YearRes}1231_restart.nc' 325 309 liste_beg.append (file_DYN_beg) 326 310 config['Files']['file_DYN_beg'] = file_DYN_beg 327 311 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' 312 if wu.unDefined ( 'file_DYN_end' ) : 313 if LMDZ : file_DYN_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restart.nc' 314 if ICO : file_DYN_end = f'{RunDir}/ICO_{JobName}_{YearEnd}1231_restart.nc' 333 315 liste_end.append ( file_DYN_end ) 334 316 config['Files']['file_DYN_end'] = file_DYN_end 335 317 336 if file_SRF_beg == None :337 file_SRF_beg = f'{ TmpDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'318 if wu.unDefined ( 'file_SRF_beg' ) : 319 file_SRF_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 338 320 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'321 if wu.unDefined ( 'file_SRF_end' ) : 322 file_SRF_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 341 323 config['Files']['file_SRF_end'] = file_SRF_end 342 324 … … 344 326 liste_end.append ( file_SRF_end ) 345 327 346 echo ( f'{file_ATM_beg}') 347 echo ( f'{file_ATM_end}') 348 echo ( f'{file_DYN_beg}') 349 echo ( f'{file_DYN_end}') 350 echo ( f'{file_SRF_beg}') 351 echo ( f'{file_SRF_end}') 328 echo ( f'{file_ATM_beg = }') 329 echo ( f'{file_ATM_end = }') 330 echo ( f'{file_DYN_beg = }') 331 echo ( f'{file_DYN_end = }') 332 echo ( f'{file_SRF_beg = }') 333 echo ( f'{file_SRF_end = }') 334 335 if ICO : 336 if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 337 config['Files'] ['file_DYN_aire'] = file_DYN_aire 352 338 353 339 if Routing == 'SIMPLE' : 354 if file_RUN_beg == None:355 file_RUN_beg = f'{ TmpDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc'340 if wu.unDefined ( 'file_RUN_beg' ) : 341 file_RUN_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc' 356 342 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'343 if wu.unDefined ( 'file_RUN_end' ) : 344 file_RUN_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 359 345 config['Files']['file_RUN_end'] = file_RUN_end 360 346 361 347 liste_beg.append ( file_RUN_beg ) 362 348 liste_end.append ( file_RUN_end ) 363 echo ( f'{file_RUN_beg }')364 echo ( f'{file_RUN_end }')349 echo ( f'{file_RUN_beg = }' ) 350 echo ( f'{file_RUN_end = }' ) 365 351 366 352 echo ('\nExtract restart files from tar : ATM, ICO and SRF') 367 for resFile in liste_beg : 368 if os.path.exists ( os.path.join (TmpDir, resFile) ) : 369 echo ( f'file found : {resFile}' ) 353 354 def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, RunDirComp=RunDir, ErrorCount=ErrorCount ) : 355 echo ( f'----------') 356 echo ( f'file to extract : {file_name = }' ) 357 if os.path.exists ( os.path.join (RunDir, file_name) ) : 358 echo ( f'file found : {file_name = }' ) 370 359 else : 371 base_resFile = os.path.basename (resFile) 372 command = f'cd {TmpDir} ; tar xf {file_restart_beg} {base_resFile}' 373 echo ( command ) 374 os.system ( command ) 375 376 for resFile in liste_end : 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}' 382 echo ( command ) 383 os.system ( command ) 384 360 echo ( f'file not found : {file_name = }' ) 361 base_resFile = os.path.basename (file_name) 362 if os.path.exists ( tar_restart ) : 363 command = f'cd {RunDir} ; tar xf {tar_restart} {base_resFile}' 364 echo ( f'{command = }' ) 365 try : 366 os.system ( command ) 367 except : 368 if ContinueOnError : 369 ErrorCount += 1 370 echo ( f'****** Command failed : {command}' ) 371 echo ( f'****** Trying to continue' ) 372 echo ( ' ') 373 else : 374 raise Exception ( f'**** command failed : {command} - Stopping' ) 375 else : 376 echo ( f'tar done : {base_resFile}' ) 377 else : 378 echo ( f'****** Tar restart file {tar_restart = } not found ' ) 379 if ContinueOnError : 380 ErrorCount += 1 381 else : 382 raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' ) 383 return ErrorCount 384 385 ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, RunDirComp=RunDir ) 386 ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, RunDirComp=RunDir ) 387 ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, RunDirComp=RunDir ) 388 389 ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, RunDirComp=RunDir ) 390 ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, RunDirComp=RunDir ) 391 ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, RunDirComp=RunDir ) 392 393 if Routing == 'SIMPLE' : 394 ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, RunDirComp=RunDir ) 395 ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, RunDirComp=RunDir ) 396 397 ##-- Exit in case of error in the opening file phase 398 if ErrorCount > 0 : 399 echo ( ' ' ) 400 raise Exception ( f'**** Some files missing - Stopping - {ErrorCount = }' ) 401 402 ## 385 403 echo ('\nOpening ATM SRF and ICO restart files') 386 d_ATM_beg = xr.open_dataset ( os.path.join ( TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze()387 d_ATM_end = xr.open_dataset ( os.path.join ( TmpDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze()388 d_SRF_beg = xr.open_dataset ( os.path.join ( TmpDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze()389 d_SRF_end = xr.open_dataset ( os.path.join ( TmpDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze()390 d_DYN_beg = xr.open_dataset ( os.path.join ( TmpDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze()391 d_DYN_end = xr.open_dataset ( os.path.join ( TmpDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze()404 d_ATM_beg = xr.open_dataset ( os.path.join (RunDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze() 405 d_ATM_end = xr.open_dataset ( os.path.join (RunDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze() 406 d_SRF_beg = xr.open_dataset ( os.path.join (RunDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze() 407 d_SRF_end = xr.open_dataset ( os.path.join (RunDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze() 408 d_DYN_beg = xr.open_dataset ( os.path.join (RunDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze() 409 d_DYN_end = xr.open_dataset ( os.path.join (RunDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze() 392 410 393 411 for var in d_SRF_beg.variables : … … 396 414 397 415 if Routing == 'SIMPLE' : 398 d_RUN_beg = xr.open_dataset ( os.path.join (TmpDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze() 399 d_RUN_end = xr.open_dataset ( os.path.join (TmpDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze() 400 401 echo ( file_ATM_beg ) 402 echo ( file_ATM_end ) 403 echo ( file_DYN_beg ) 404 echo ( file_DYN_end ) 405 echo ( file_SRF_beg ) 406 echo ( file_SRF_end ) 416 d_RUN_beg = xr.open_dataset ( os.path.join (RunDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze() 417 d_RUN_end = xr.open_dataset ( os.path.join (RunDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze() 418 419 def to_cell ( dd, newname='cell' ) : 420 '''Set space dimension to 'cell' ''' 421 for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] : 422 try : dd = dd.rename ( {oldname : newname} ) 423 except : pass 424 return dd 425 426 d_ATM_beg = to_cell ( d_ATM_beg ) 427 d_ATM_end = to_cell ( d_ATM_end ) 428 d_SRF_beg = to_cell ( d_SRF_beg ) 429 d_SRF_end = to_cell ( d_SRF_end ) 430 d_DYN_beg = to_cell ( d_DYN_beg ) 431 d_DYN_end = to_cell ( d_DYN_end ) 432 407 433 if Routing == 'SIMPLE' : 408 echo ( file_RUN_beg ) 409 echo ( file_RUN_end ) 410 411 ## 434 d_RUN_beg = to_cell ( d_RUN_beg ) 435 d_RUN_end = to_cell ( d_RUN_end ) 436 437 d_ATM_his = to_cell ( d_ATM_his ) 438 d_SRF_his = to_cell ( d_SRF_his ) 439 440 echo ( f'{file_ATM_beg = }' ) 441 echo ( f'{file_ATM_end = }' ) 442 echo ( f'{file_DYN_beg = }' ) 443 echo ( f'{file_DYN_end = }' ) 444 echo ( f'{file_SRF_beg = }' ) 445 echo ( f'{file_SRF_end = }' ) 446 if Routing == 'SIMPLE' : 447 echo ( f'{file_RUN_beg = }' ) 448 echo ( f'{file_RUN_end = }' ) 449 450 ## Write the full configuration 412 451 config_out = open (FullIniFile, 'w') 413 config.write ( config_out )452 config.write ( config_out ) 414 453 config_out.close () 415 454 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 449 def ATM_stock_int (stock) : 450 '''Integrate stock on atmosphere grid''' 451 ATM_stock_int = Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 452 return ATM_stock_int 453 454 def ATM_flux_int (flux) : 455 '''Integrate flux on atmosphere grid''' 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()) ) 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 468 469 def ONE_stock_int (stock) : 470 '''Sum stock ''' 471 ONE_stock_int = Psum ( stock.to_masked_array().ravel() ) 472 return ONE_stock_int 473 474 def ONE_flux_int (flux) : 475 '''Sum flux ''' 476 ONE_flux_int = Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 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 486 455 if ICO : 456 if SortIco : 457 # Creation d'une clef de tri pour les restarts (pour chaque restart !) 458 DYN_beg_keysort = np.lexsort ( (d_DYN_beg['lat_mesh'], d_DYN_beg['lon_mesh'] ) ) 459 ATM_beg_keysort = np.lexsort ( (d_ATM_beg['latitude'], d_ATM_beg['longitude']) ) 460 SRF_beg_keysort = np.lexsort ( (d_SRF_beg['nav_lat' ], d_SRF_beg['nav_lon'] ) ) 461 DYN_end_keysort = np.lexsort ( (d_DYN_end['lat_mesh'], d_DYN_end['lon_mesh'] ) ) 462 ATM_end_keysort = np.lexsort ( (d_ATM_end['latitude'], d_ATM_end['longitude']) ) 463 SRF_end_keysort = np.lexsort ( (d_SRF_end['nav_lat' ], d_SRF_end['nav_lon'] ) ) 464 465 if ATM_HIS == 'ico' : 466 ATM_his_keysort = np.lexsort ( (d_ATM_his['lat'], d_ATM_his['lon'] ) ) 467 if SRF_HIS == 'ico' : 468 SRF_his_keysort = np.lexsort ( (d_SRF_his['lat'], d_SRF_his['lon'] ) ) 469 RUN_his_keysort = SRF_his_keysort 470 else : 471 DYN_beg_keysort = np.arange ( len ( d_DYN_beg['lat_mesh'] ) ) 472 ATM_beg_keysort = DYN_beg_keysort 473 SRF_beg_keysort = DYN_beg_keysort 474 475 DYN_end_keysort = DYN_beg_keysort 476 ATM_end_keysort = DYN_beg_keysort 477 SRF_end_keysort = DYN_beg_keysort 478 479 ATM_his_keysort = DYN_beg_keysort 480 SRF_his_keysort = DYN_beg_keysort 481 RUN_his_keysort = SRF_his_keysort 482 487 483 # ATM grid with cell surfaces 484 if LMDZ : 485 #ATM_lat = lmdz.geo2point ( d_ATM_his ['lat'] , dim1D='cell' ) 486 #ATM_lon = lmdz.geo2point ( d_ATM_his ['lon'] , dim1D='cell' ) 487 ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'] [0], cumulPoles=True, dim1D='cell' ) 488 ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell' ) 489 ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell' ) 490 ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell' ) 491 ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell' ) 492 #SRF_lat = lmdz.geo2point ( d_SRF_his ['lat'] [0], dim1D='cell' ) 493 #SRF_lon = lmdz.geo2point ( d_SRF_his ['lon'] [0], dim1D='cell' ) 494 495 #??SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell' ) 496 SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'], dim1D='cell' ) 497 488 498 if ICO : 489 jpja, jpia = d_ATM_his['aire'][0].shape 490 file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 491 config['Files']['file_ATM_aire'] = file_ATM_aire 492 echo ( f'Aire sur grille reguliere : {file_ATM_aire}' ) 493 d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 494 ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True ) 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 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 504 505 if LMDZ : 506 ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=True ) 507 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 508 ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 509 #SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 510 SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'] ) 511 512 SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.) 513 499 if ATM_HIS == 'latlon' : 500 jpja, jpia = d_ATM_his['aire'][0].shape 501 file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 502 config['Files']['file_ATM_aire'] = file_ATM_aire 503 echo ( f'Aire sur grille reguliere : {file_ATM_aire = }' ) 504 d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 505 #ATM_lat = d_ATM_his ['lat'] 506 #ATM_lon = d_ATM_his ['lon'] 507 ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True, dim1D='cell' ) 508 ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell' ) 509 ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell' ) 510 ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell' ) 511 ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell' ) 512 #SRF_lat = d_SRF_his ['lat'] 513 #SRF_lon = d_SRF_his ['lon'] 514 SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell' ) 515 516 if ATM_HIS == 'ico' : 517 echo ( f'Grille icosaedre' ) 518 ATM_aire = d_ATM_his ['aire'] [0][ATM_his_keysort].squeeze() 519 #ATM_lat = d_ATM_his ['lat'] [ATM_his_keysort] 520 #ATM_lon = d_ATM_his ['lon'] [ATM_his_keysort] 521 ATM_fter = d_ATM_his ['fract_ter'][0][ATM_his_keysort] 522 ATM_foce = d_ATM_his ['fract_oce'][0][ATM_his_keysort] 523 ATM_fsic = d_ATM_his ['fract_sic'][0][ATM_his_keysort] 524 ATM_flic = d_ATM_his ['fract_lic'][0][ATM_his_keysort] 525 #SRF_lat = d_SRF_his ['lat'] [SRF_his_keysort] 526 #SRF_lon = d_SRF_his ['lon'] [SRF_his_keysort] 527 SRF_aire = d_SRF_his ['Areas'][SRF_his_keysort] * d_SRF_his ['Contfrac'][SRF_his_keysort] 528 529 ATM_fsea = ATM_foce + ATM_fsic 530 ATM_flnd = ATM_fter + ATM_flic 531 ATM_aire_fter = ATM_aire * ATM_fter 532 ATM_aire_flic = ATM_aire * ATM_flic 533 ATM_aire_fsic = ATM_aire * ATM_fsic 534 ATM_aire_foce = ATM_aire * ATM_foce 535 ATM_aire_flnd = ATM_aire * ATM_flnd 536 ATM_aire_fsea = ATM_aire * ATM_fsea 537 538 SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0.) 539 514 540 if ICO : 515 541 # Area on icosahedron grid 516 file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )517 542 d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze() 518 d_DYN_aire = d_DYN_aire.rename ( {'cell':'cell_mesh'} ) 519 DYN_aire = d_DYN_aire['aire'] 520 521 DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic'] 543 544 if SortIco : 545 # Creation d'une clef de tri pour le fichier aire 546 DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) ) 547 else : 548 DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) ) 549 550 DYN_lat = d_DYN_aire['lat'][DYN_aire_keysort] 551 DYN_lon = d_DYN_aire['lon'][DYN_aire_keysort] 552 553 DYN_aire = d_DYN_aire['aire'][DYN_aire_keysort] 554 DYN_fsea = d_DYN_aire ['fract_oce'][DYN_aire_keysort] + d_DYN_aire['fract_sic'][DYN_aire_keysort] 522 555 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'}) 556 DYN_fter = d_ATM_beg['FTER'][ATM_beg_keysort] 557 DYN_flic = d_ATM_beg['FLIC'][ATM_beg_keysort] 558 DYN_aire_fter = DYN_aire * DYN_fter 525 559 526 560 if LMDZ : 561 # Area on lon/lat grid 527 562 DYN_aire = ATM_aire 528 563 DYN_fsea = ATM_fsea 529 564 DYN_flnd = ATM_flnd 530 565 DYN_fter = d_ATM_beg['FTER'] 531 DYN_flic = d_ATM_beg['FTER'] 532 533 534 DYN_aire_fter = DYN_aire * DYN_fter 566 DYN_flic = d_ATM_beg['FLIC'] 567 DYN_aire_fter = DYN_aire * DYN_fter 568 569 # Functions computing integrals and sum 570 def ATM_stock_int (stock) : 571 '''Integrate (* surface) stock on atmosphere grid''' 572 ATM_stock_int = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 573 return ATM_stock_int 574 575 def ATM_flux_int (flux) : 576 '''Integrate (* time * surface) flux on atmosphere grid''' 577 ATM_flux_int = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 578 return ATM_flux_int 579 580 def SRF_stock_int (stock) : 581 '''Integrate (* surface) stock on land grid''' 582 SRF_stock_int = wu.Ksum ( ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 583 return SRF_stock_int 584 585 def SRF_flux_int (flux) : 586 '''Integrate (* time * surface) flux on land grid''' 587 SRF_flux_int = wu.Psum ( (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 588 return SRF_flux_int 589 590 def ONE_stock_int (stock) : 591 '''Sum stock ''' 592 ONE_stock_int = wu.Psum ( stock.to_masked_array().ravel() ) 593 return ONE_stock_int 594 595 def ONE_flux_int (flux) : 596 '''Integrate (* time) flux on area=1 grid ''' 597 ONE_flux_int = wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 598 return ONE_flux_int 599 535 600 536 601 #if LMDZ : … … 538 603 539 604 ATM_aire_sea = ATM_aire * ATM_fsea 540 ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea )541 542 605 ATM_aire_tot = ONE_stock_int (ATM_aire) 543 606 ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea) 544 607 545 608 546 echo ( 'A ire atmosphere/4pi R^2: {:12.5f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )609 echo ( 'Area of atmosphere/(4pi R^2) : {:12.5f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) ) 547 610 548 611 if ( np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) : 549 raise Exception ('Err eur surface interpolee sur grille reguliere')612 raise Exception ('Error of atmosphere surface interpolated on lon/lat grid') 550 613 551 614 echo ( '\n====================================================================================' ) … … 555 618 #-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv 556 619 557 # ATM vertical grid 620 echo ( 'ATM vertical grid' ) 558 621 ATM_Ahyb = d_ATM_his['Ahyb'].squeeze() 559 622 ATM_Bhyb = d_ATM_his['Bhyb'].squeeze() 560 623 561 # Surface pressure 562 if ICO : 563 DYN_psol_beg = d_DYN_beg['ps'] 564 DYN_psol_end = d_DYN_end['ps'] 624 echo ( 'Surface pressure' ) 625 if ICO : 626 DYN_psol_beg = d_DYN_beg['ps'][DYN_beg_keysort] 627 DYN_psol_end = d_DYN_end['ps'][DYN_end_keysort] 565 628 if LMDZ : 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)629 DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 630 DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 631 632 echo ( '3D Pressure at the interface layers (not scalar points)' ) 570 633 DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg 571 634 DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end 572 635 573 636 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] 637 cell = DYN_psol_beg.shape[-1] 576 638 klev = klevp1 - 1 577 639 578 # Layer thickness (pressure) 579 if ICO : 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) ) ) 582 if LMDZ : 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) ) ) 585 640 echo ( 'Layer thickness (pressure)' ) 641 DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 642 DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 643 #DYN_dpres_beg = DYN_dpres_beg 644 #DYN_dpres_end = DYN_dpres_end 645 586 646 for k in np.arange (klevp1-1) : 587 647 DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] 588 648 DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:] 589 649 590 ##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases 650 echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' ) 591 651 if LMDZ : 592 652 try : 593 DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'] ).isel(rlonv=slice(0,-1) ) )594 DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'] ).isel(rlonv=slice(0,-1) ) )653 DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 654 DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 595 655 except : 596 DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) )597 DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) )656 DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 657 DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 598 658 if ICO : 599 659 try : 600 DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).rename ( {'lev':'sigs'} ) 601 DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).rename ( {'lev':'sigs'} ) 660 DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'])[..., DYN_beg_keysort] 661 DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'])[..., DYN_beg_keysort] 662 602 663 except : 603 DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) ).rename ( {'lev':'sigs'} ) 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'} ) 605 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 664 DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) )[..., DYN_beg_keysort] 665 DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) )[..., DYN_beg_keysort] 666 667 DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} ) 668 DYN_wat_beg = DYN_wat_end.rename ( {'lev':'sigs'} ) 669 670 echo ( 'Compute water content : vertical and horizontal integral' ) 671 DYN_mas_wat_beg = ATM_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav 672 DYN_mas_wat_end = ATM_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav 673 674 echo ( 'Variation of water content' ) 611 675 dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 612 676 … … 614 678 # \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear kg2myear(\1) 615 679 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} ' ) 680 echo ( f'\nChange of atmosphere water content (dynamics) -- {Title} ' ) 631 681 echo ( '------------------------------------------------------------------------------------' ) 632 682 echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 633 683 prtFlux ( 'dMass (atm) ', dDYN_mas_wat, 'e', True ) 634 684 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']636 ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] +d_ATM_end['SNOW02']*d_ATM_end['FLIC']+d_ATM_end['SNOW03']*d_ATM_end['FOCE']+d_ATM_end['SNOW04']*d_ATM_end['FSIC']637 638 ATM_qs_beg = d_ATM_beg['QS01'] *d_ATM_beg['FTER']+d_ATM_beg['QS02']*d_ATM_beg['FLIC']+d_ATM_beg['QS03']*d_ATM_beg['FOCE']+d_ATM_beg['QS04']*d_ATM_beg['FSIC']639 ATM_qs_end = d_ATM_end['QS01'] *d_ATM_end['FTER']+d_ATM_end['QS02']*d_ATM_end['FLIC']+d_ATM_end['QS03']*d_ATM_end['FOCE']+d_ATM_end['QS04']*d_ATM_end['FSIC']685 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'] 686 ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC'] 687 688 ATM_qs_beg = d_ATM_beg['QS01'] *d_ATM_beg['FTER'] + d_ATM_beg['QS02'] *d_ATM_beg['FLIC'] + d_ATM_beg['QS03'] *d_ATM_beg['FOCE'] + d_ATM_beg['QS04'] *d_ATM_beg['FSIC'] 689 ATM_qs_end = d_ATM_end['QS01'] *d_ATM_end['FTER'] + d_ATM_end['QS02'] *d_ATM_end['FLIC'] + d_ATM_end['QS03'] *d_ATM_end['FOCE'] + d_ATM_end['QS04'] *d_ATM_end['FSIC'] 640 690 641 691 ATM_qsol_beg = d_ATM_beg['QSOL'] 642 692 ATM_qsol_end = d_ATM_end['QSOL'] 693 694 LIC_sno_beg = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] 695 LIC_sno_end = d_ATM_end['SNOW02']*d_ATM_end['FLIC'] 696 697 LIC_runlic0_beg = d_ATM_beg['RUNOFFLIC0'] 698 LIC_runlic0_end = d_ATM_end['RUNOFFLIC0'] 643 699 644 700 ATM_qs01_beg = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] … … 651 707 ATM_qs04_end = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 652 708 709 # if LMDZ : 710 # ATM_sno_beg = lmdz.geo2point ( ATM_sno_beg , dim1D='cell' ) 711 # ATM_sno_end = lmdz.geo2point ( ATM_sno_end , dim1D='cell' ) 712 # ATM_qs_beg = lmdz.geo2point ( ATM_qs_beg , dim1D='cell' ) 713 # ATM_qs_end = lmdz.geo2point ( ATM_qs_end , dim1D='cell' ) 714 # ATM_qsol_beg = lmdz.geo2point ( ATM_qsol_beg , dim1D='cell' ) 715 # ATM_qsol_end = lmdz.geo2point ( ATM_qsol_end , dim1D='cell' ) 716 # ATM_qs01_beg = lmdz.geo2point ( ATM_qs01_beg , dim1D='cell' ) 717 # ATM_qs01_end = lmdz.geo2point ( ATM_qs01_end , dim1D='cell' ) 718 # ATM_qs02_beg = lmdz.geo_point ( ATM_qs02_beg , dim1D='cell' ) 719 # ATM_qs02_end = lmdz.geo2point ( ATM_qs02_end , dim1D='cell' ) 720 # ATM_qs03_beg = lmdz.geo2point ( ATM_qs03_beg , dim1D='cell' ) 721 # ATM_qs03_end = lmdz.geo2point ( ATM_qs03_end , dim1D='cell' ) 722 # ATM_qs04_beg = lmdz.geo2point ( ATM_qs04_beg , dim1D='cell' ) 723 # ATM_qs04_end = lmdz.geo2point ( ATM_qs04_end , dim1D='cell' ) 724 # LIC_sno_beg = lmdz.geo2point ( LIC_sno_beg , dim1D='cell' ) 725 # LIC_sno_end = lmdz.geo2point ( LIC_sno_end , dim1D='cell' ) 726 # LIC_runlic0_beg = lmdz.geo2point ( LIC_runlic0_beg, dim1D='cell' ) 727 # LIC_runlic0_end = lmdz.geo2point ( LIC_runlic0_end, dim1D='cell' ) 653 728 if ICO : 654 ATM_sno_beg = ATM_sno_beg .rename ( {'points_physiques':'cell_mesh'} ) 655 ATM_sno_end = ATM_sno_end .rename ( {'points_physiques':'cell_mesh'} ) 656 ATM_qs_beg = ATM_qs_beg .rename ( {'points_physiques':'cell_mesh'} ) 657 ATM_qs_end = ATM_qs_end .rename ( {'points_physiques':'cell_mesh'} ) 658 ATM_qsol_beg = ATM_qsol_beg.rename ( {'points_physiques':'cell_mesh'} ) 659 ATM_qsol_end = ATM_qsol_end.rename ( {'points_physiques':'cell_mesh'} ) 660 ATM_qs01_beg = ATM_qs01_beg.rename ( {'points_physiques':'cell_mesh'} ) 661 ATM_qs01_end = ATM_qs01_end.rename ( {'points_physiques':'cell_mesh'} ) 662 ATM_qs02_beg = ATM_qs02_beg.rename ( {'points_physiques':'cell_mesh'} ) 663 ATM_qs02_end = ATM_qs02_end.rename ( {'points_physiques':'cell_mesh'} ) 664 ATM_qs03_beg = ATM_qs03_beg.rename ( {'points_physiques':'cell_mesh'} ) 665 ATM_qs03_end = ATM_qs03_end.rename ( {'points_physiques':'cell_mesh'} ) 666 ATM_qs04_beg = ATM_qs04_beg.rename ( {'points_physiques':'cell_mesh'} ) 667 ATM_qs04_end = ATM_qs04_end.rename ( {'points_physiques':'cell_mesh'} ) 668 669 ATM_mas_sno_beg = ATM_stock_int ( ATM_sno_beg ) 670 ATM_mas_sno_end = ATM_stock_int ( ATM_sno_end ) 671 ATM_mas_qs_beg = ATM_stock_int ( ATM_qs_beg ) 672 ATM_mas_qs_end = ATM_stock_int ( ATM_qs_end ) 673 ATM_mas_qsol_beg = ATM_stock_int ( ATM_qsol_beg ) 674 ATM_mas_qsol_end = ATM_stock_int ( ATM_qsol_end ) 675 ATM_mas_qs01_beg = ATM_stock_int ( ATM_qs01_beg ) 676 ATM_mas_qs01_end = ATM_stock_int ( ATM_qs01_end ) 677 ATM_mas_qs02_beg = ATM_stock_int ( ATM_qs02_beg ) 678 ATM_mas_qs02_end = ATM_stock_int ( ATM_qs02_end ) 679 ATM_mas_qs03_beg = ATM_stock_int ( ATM_qs03_beg ) 680 ATM_mas_qs03_end = ATM_stock_int ( ATM_qs03_end ) 681 ATM_mas_qs04_beg = ATM_stock_int ( ATM_qs04_beg ) 682 ATM_mas_qs04_end = ATM_stock_int ( ATM_qs04_end ) 729 ATM_sno_beg = ATM_sno_beg [ATM_beg_keysort] 730 ATM_sno_end = ATM_sno_end [ATM_end_keysort] 731 ATM_qs_beg = ATM_qs_beg [ATM_beg_keysort] 732 ATM_qs_end = ATM_qs_end [ATM_end_keysort] 733 ATM_qsol_beg = ATM_qsol_beg [ATM_beg_keysort] 734 ATM_qsol_end = ATM_qsol_end [ATM_end_keysort] 735 ATM_qs01_beg = ATM_qs01_beg [ATM_beg_keysort] 736 ATM_qs01_end = ATM_qs01_end [ATM_end_keysort] 737 ATM_qs02_beg = ATM_qs02_beg [ATM_beg_keysort] 738 ATM_qs02_end = ATM_qs02_end [ATM_end_keysort] 739 ATM_qs03_beg = ATM_qs03_beg [ATM_beg_keysort] 740 ATM_qs03_end = ATM_qs03_end [ATM_end_keysort] 741 ATM_qs04_beg = ATM_qs04_beg [ATM_beg_keysort] 742 ATM_qs04_end = ATM_qs04_end [ATM_end_keysort] 743 LIC_sno_beg = LIC_sno_beg [ATM_beg_keysort] 744 LIC_sno_end = LIC_sno_end [ATM_end_keysort] 745 LIC_runlic0_beg = LIC_runlic0_beg[ATM_beg_keysort] 746 LIC_runlic0_end = LIC_runlic0_end[ATM_end_keysort] 747 748 749 750 LIC_qs_beg = ATM_qs02_beg 751 LIC_qs_end = ATM_qs02_end 752 753 ATM_mas_sno_beg = ATM_stock_int ( ATM_sno_beg ) 754 ATM_mas_sno_end = ATM_stock_int ( ATM_sno_end ) 755 756 ATM_mas_qs_beg = ATM_stock_int ( ATM_qs_beg ) 757 ATM_mas_qs_end = ATM_stock_int ( ATM_qs_end ) 758 ATM_mas_qsol_beg = ATM_stock_int ( ATM_qsol_beg ) 759 ATM_mas_qsol_end = ATM_stock_int ( ATM_qsol_end ) 760 ATM_mas_qs01_beg = ATM_stock_int ( ATM_qs01_beg ) 761 ATM_mas_qs01_end = ATM_stock_int ( ATM_qs01_end ) 762 ATM_mas_qs02_beg = ATM_stock_int ( ATM_qs02_beg ) 763 ATM_mas_qs02_end = ATM_stock_int ( ATM_qs02_end ) 764 ATM_mas_qs03_beg = ATM_stock_int ( ATM_qs03_beg ) 765 ATM_mas_qs03_end = ATM_stock_int ( ATM_qs03_end ) 766 ATM_mas_qs04_beg = ATM_stock_int ( ATM_qs04_beg ) 767 ATM_mas_qs04_end = ATM_stock_int ( ATM_qs04_end ) 768 769 LIC_mas_sno_beg = ATM_stock_int ( LIC_sno_beg ) 770 LIC_mas_sno_end = ATM_stock_int ( LIC_sno_end ) 771 LIC_mas_runlic0_beg = ATM_stock_int ( LIC_runlic0_beg ) 772 LIC_mas_runlic0_end = ATM_stock_int ( LIC_runlic0_end ) 773 774 LIC_mas_qs_beg = ATM_mas_qs02_beg 775 LIC_mas_qs_end = ATM_mas_qs02_end 776 LIC_mas_wat_beg = LIC_mas_qs_beg + LIC_mas_sno_beg 777 LIC_mas_wat_end = LIC_mas_qs_end + LIC_mas_sno_end 683 778 684 779 dATM_mas_sno = ATM_mas_sno_end - ATM_mas_sno_beg … … 691 786 dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg 692 787 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} ' ) 788 dLIC_mas_qs = LIC_mas_qs_end - LIC_mas_qs_beg 789 dLIC_mas_sno = LIC_mas_sno_end - LIC_mas_sno_beg 790 dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg 791 dLIC_mas_wat = dLIC_mas_qs + dLIC_mas_sno 792 793 echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' ) 696 794 echo ( '------------------------------------------------------------------------------------' ) 697 795 echo ( 'ATM_mas_sno_beg = {:12.6e} kg | ATM_mas_sno_end = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 698 796 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} ' ) 797 798 echo ( f'\nChange of soil humidity content -- {Title} ' ) 702 799 echo ( '------------------------------------------------------------------------------------' ) 703 800 echo ( 'ATM_mas_qs_beg = {:12.6e} kg | ATM_mas_qs_end = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) ) 704 801 prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True ) 705 802 706 echo ( f'\n Variation du contenu en eau+neige atmosphere-- {Title} ' )803 echo ( f'\nChange of atmosphere water+snow content -- {Title} ' ) 707 804 echo ( '------------------------------------------------------------------------------------' ) 708 805 prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True) … … 715 812 RUN_mas_wat_slow_beg = ONE_stock_int ( d_RUN_beg ['slow_reservoir'] ) 716 813 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 814 815 RUN_mas_wat_flood_beg = ONE_stock_int ( d_SRF_beg ['floodres'] ) 816 RUN_mas_wat_lake_beg = ONE_stock_int ( d_SRF_beg ['lakeres'] ) 817 RUN_mas_wat_pond_beg = ONE_stock_int ( d_SRF_beg ['pondres'] ) 720 818 721 819 RUN_mas_wat_fast_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] ) 722 820 RUN_mas_wat_slow_end = ONE_stock_int ( d_RUN_end ['slow_reservoir'] ) 723 821 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'] ) 822 740 823 RUN_mas_wat_flood_end = ONE_stock_int ( d_SRF_end ['floodres'] ) 741 824 RUN_mas_wat_lake_end = ONE_stock_int ( d_SRF_end ['lakeres'] ) 742 825 RUN_mas_wat_pond_end = ONE_stock_int ( d_SRF_end ['pondres'] ) 826 827 if Routing == 'SECHIBA' : 828 RUN_mas_wat_fast_beg = ONE_stock_int ( d_SRF_beg ['fastres'] ) 829 RUN_mas_wat_slow_beg = ONE_stock_int ( d_SRF_beg ['slowres'] ) 830 RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 831 RUN_mas_wat_flood_beg = ONE_stock_int ( d_SRF_beg ['floodres'] ) 832 RUN_mas_wat_lake_beg = ONE_stock_int ( d_SRF_beg ['lakeres'] ) 833 RUN_mas_wat_pond_beg = ONE_stock_int ( d_SRF_beg ['pondres'] ) 834 835 RUN_mas_wat_fast_end = ONE_stock_int ( d_SRF_end ['fastres'] ) 836 RUN_mas_wat_slow_end = ONE_stock_int ( d_SRF_end ['slowres'] ) 837 RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 838 RUN_mas_wat_flood_end = ONE_stock_int ( d_SRF_end ['floodres'] ) 839 RUN_mas_wat_lake_end = ONE_stock_int ( d_SRF_end ['lakeres'] ) 840 RUN_mas_wat_pond_end = ONE_stock_int ( d_SRF_end ['pondres'] ) 743 841 744 842 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 … … 754 852 dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg 755 853 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 ) ) 854 echo ( f'\nRunoff reservoirs -- {Title} ') 855 echo ( f'------------------------------------------------------------------------------------' ) 856 echo ( f'RUN_mas_wat_fast_beg = {RUN_mas_wat_fast_beg :12.6e} kg | RUN_mas_wat_fast_end = {RUN_mas_wat_fast_end :12.6e} kg ' ) 857 echo ( f'RUN_mas_wat_slow_beg = {RUN_mas_wat_slow_beg :12.6e} kg | RUN_mas_wat_slow_end = {RUN_mas_wat_slow_end :12.6e} kg ' ) 858 echo ( f'RUN_mas_wat_stream_beg = {RUN_mas_wat_stream_beg:12.6e} kg | RUN_mas_wat_stream_end = {RUN_mas_wat_stream_end:12.6e} kg ' ) 859 echo ( f'RUN_mas_wat_flood_beg = {RUN_mas_wat_flood_beg :12.6e} kg | RUN_mas_wat_flood_end = {RUN_mas_wat_flood_end :12.6e} kg ' ) 860 echo ( f'RUN_mas_wat_lake_beg = {RUN_mas_wat_lake_beg :12.6e} kg | RUN_mas_wat_lake_end = {RUN_mas_wat_lake_end :12.6e} kg ' ) 861 echo ( f'RUN_mas_wat_pond_beg = {RUN_mas_wat_pond_beg :12.6e} kg | RUN_mas_wat_pond_end = {RUN_mas_wat_pond_end :12.6e} kg ' ) 862 echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_beg :12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end :12.6e} kg ' ) 764 863 765 864 echo ( '------------------------------------------------------------------------------------' ) … … 776 875 echo ( '------------------------------------------------------------------------------------' ) 777 876 echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) ) 778 prtFlux ( 'dMass (routing) ', dRUN_mas_wat 877 prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True ) 779 878 780 879 echo ( '\n====================================================================================' ) 781 880 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.)881 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.) 882 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.) 883 SRF_snow_beg = d_SRF_beg['snow_beg'] ; SRF_snow_beg = SRF_snow_beg .where (SRF_snow_beg < 1E15, 0.) 884 SRF_lakeres_beg = d_SRF_beg['lakeres'] ; SRF_lakeres_beg = SRF_lakeres_beg .where (SRF_lakeres_beg < 1E15, 0.) 885 886 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.) 887 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.) 888 SRF_snow_end = d_SRF_end['snow_beg'] ; SRF_snow_end = SRF_snow_end .where (SRF_snow_end < 1E15, 0.) 889 SRF_lakeres_end = d_SRF_end['lakeres'] ; SRF_lakeres_end = SRF_lakeres_end .where (SRF_lakeres_end < 1E15, 0.) 791 890 792 891 if LMDZ : 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 )892 SRF_tot_watveg_beg = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell') 893 SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell') 894 SRF_snow_beg = lmdz.geo2point (SRF_snow_beg , dim1D='cell') 895 SRF_lakeres_beg = lmdz.geo2point (SRF_lakeres_beg , dim1D='cell') 896 SRF_tot_watveg_end = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell') 897 SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell') 898 SRF_snow_end = lmdz.geo2point (SRF_snow_end , dim1D='cell') 899 SRF_lakeres_end = lmdz.geo2point (SRF_lakeres_end , dim1D='cell') 801 900 802 901 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'} )902 SRF_tot_watveg_beg = SRF_tot_watveg_beg [SRF_beg_keysort] 903 SRF_tot_watsoil_beg = SRF_tot_watsoil_beg [SRF_beg_keysort] 904 SRF_snow_beg = SRF_snow_beg [SRF_beg_keysort] 905 SRF_lakeres_beg = SRF_lakeres_beg [SRF_beg_keysort] 906 SRF_tot_watveg_end = SRF_tot_watveg_end [SRF_end_keysort] 907 SRF_tot_watsoil_end = SRF_tot_watsoil_end [SRF_end_keysort] 908 SRF_snow_end = SRF_snow_end [SRF_end_keysort] 909 SRF_lakeres_end = SRF_lakeres_end [SRF_end_keysort] 811 910 812 911 # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood … … 844 943 845 944 846 847 945 echo ( '------------------------------------------------------------------------------------' ) 848 echo ( f'\n Les differentsreservoirs -- {Title} ')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 ))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 ))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 ))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 ))946 echo ( f'\nSurface reservoirs -- {Title} ') 947 echo ( f'SRF_mas_watveg_beg = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end = {SRF_mas_watveg_end :12.6e} kg ' ) 948 echo ( f'SRF_mas_watsoil_beg = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end = {SRF_mas_watsoil_end:12.6e} kg ' ) 949 echo ( f'SRF_mas_snow_beg = {SRF_mas_snow_beg :12.6e} kg | SRF_mas_snow_end = {SRF_mas_snow_end :12.6e} kg ' ) 950 echo ( f'SRF_mas_lake_beg = {SRF_mas_lake_beg :12.6e} kg | SRF_mas_lake_end = {SRF_mas_lake_end :12.6e} kg ' ) 853 951 854 952 prtFlux ('dMass (watveg) ', dSRF_mas_watveg , 'e' , True) … … 863 961 echo ( '------------------------------------------------------------------------------------' ) 864 962 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))963 echo ( f'SRF_mas_wat_beg = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end = {SRF_mas_wat_end:12.6e} kg ') 866 964 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 867 965 … … 876 974 echo ( f'-- ATM Fluxes -- {Title} ' ) 877 975 878 ATM_wbilo_oce = lmdz.geo2point ( d_ATM_his ['wbilo_oce'] ) 879 ATM_wbilo_sic = lmdz.geo2point ( d_ATM_his ['wbilo_sic'] ) 880 ATM_wbilo_ter = lmdz.geo2point ( d_ATM_his ['wbilo_ter'] ) 881 ATM_wbilo_lic = lmdz.geo2point ( d_ATM_his ['wbilo_lic'] ) 882 ATM_runofflic = lmdz.geo2point ( d_ATM_his ['runofflic'] ) 883 ATM_fqcalving = lmdz.geo2point ( d_ATM_his ['fqcalving'] ) 884 ATM_fqfonte = lmdz.geo2point ( d_ATM_his ['fqfonte'] ) 885 ATM_precip = lmdz.geo2point ( d_ATM_his ['precip'] ) 886 ATM_snowf = lmdz.geo2point ( d_ATM_his ['snow'] ) 887 ATM_evap = lmdz.geo2point ( d_ATM_his ['evap'] ) 976 print ( 'Step 1', end='' ) 977 if ATM_HIS == 'latlon' : 978 ATM_wbilo_oce = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell' ) 979 ATM_wbilo_sic = lmdz.geo2point ( d_ATM_his ['wbilo_sic'], dim1D='cell' ) 980 ATM_wbilo_ter = lmdz.geo2point ( d_ATM_his ['wbilo_ter'], dim1D='cell' ) 981 ATM_wbilo_lic = lmdz.geo2point ( d_ATM_his ['wbilo_lic'], dim1D='cell' ) 982 ATM_runofflic = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell' ) 983 ATM_fqcalving = lmdz.geo2point ( d_ATM_his ['fqcalving'], dim1D='cell' ) 984 ATM_fqfonte = lmdz.geo2point ( d_ATM_his ['fqfonte'] , dim1D='cell' ) 985 ATM_precip = lmdz.geo2point ( d_ATM_his ['precip'] , dim1D='cell' ) 986 ATM_snowf = lmdz.geo2point ( d_ATM_his ['snow'] , dim1D='cell' ) 987 ATM_evap = lmdz.geo2point ( d_ATM_his ['evap'] , dim1D='cell' ) 988 ATM_wevap_ter = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell' ) 989 ATM_wevap_oce = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell' ) 990 ATM_wevap_lic = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell' ) 991 ATM_wevap_sic = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell' ) 992 ATM_runofflic = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell' ) 993 994 print ( ' - Step 2', end='' ) 995 if ATM_HIS == 'ico' : 996 ATM_wbilo_oce = d_ATM_his ['wbilo_oce'][..., ATM_his_keysort] 997 ATM_wbilo_sic = d_ATM_his ['wbilo_sic'][..., ATM_his_keysort] 998 ATM_wbilo_ter = d_ATM_his ['wbilo_ter'][..., ATM_his_keysort] 999 ATM_wbilo_lic = d_ATM_his ['wbilo_lic'][..., ATM_his_keysort] 1000 ATM_runofflic = d_ATM_his ['runofflic'][..., ATM_his_keysort] 1001 ATM_fqcalving = d_ATM_his ['fqcalving'][..., ATM_his_keysort] 1002 ATM_fqfonte = d_ATM_his ['fqfonte'] [..., ATM_his_keysort] 1003 ATM_precip = d_ATM_his ['precip'] [..., ATM_his_keysort] 1004 ATM_snowf = d_ATM_his ['snow'] [..., ATM_his_keysort] 1005 ATM_evap = d_ATM_his ['evap'] [..., ATM_his_keysort] 1006 ATM_wevap_ter = d_ATM_his ['wevap_ter'][..., ATM_his_keysort] 1007 ATM_wevap_oce = d_ATM_his ['wevap_oce'][..., ATM_his_keysort] 1008 ATM_wevap_lic = d_ATM_his ['wevap_lic'][..., ATM_his_keysort] 1009 ATM_wevap_sic = d_ATM_his ['wevap_sic'][..., ATM_his_keysort] 1010 ATM_runofflic = d_ATM_his ['runofflic'][..., ATM_his_keysort] 1011 1012 print ( ' - Step 3', end='' ) 888 1013 ATM_wbilo = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic 889 1014 ATM_emp = ATM_evap - ATM_precip 890 891 RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'] ) 892 RUN_riverflow = lmdz.geo2point ( d_RUN_his ['riverflow'] ) 893 RUN_runoff = lmdz.geo2point ( d_RUN_his ['runoff'] ) 894 RUN_drainage = lmdz.geo2point ( d_RUN_his ['drainage'] ) 895 RUN_riversret = lmdz.geo2point ( d_RUN_his ['riversret'] ) 896 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'] ) 902 SRF_snowf = lmdz.geo2point ( d_SRF_his ['snowf'] ) 903 SRF_subli = lmdz.geo2point ( d_SRF_his ['subli'] ) 904 SRF_transpir = lmdz.geo2point ( np.sum(d_SRF_his ['transpir'], axis=1) ) ; SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 1015 ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic 1016 ATM_wevap_sea = ATM_wevap_sic + ATM_wevap_oce 1017 1018 ATM_wemp_ter = ATM_wevap_ter - ATM_precip * ATM_fter 1019 ATM_wemp_oce = ATM_wevap_oce - ATM_precip * ATM_foce 1020 ATM_wemp_sic = ATM_wevap_sic - ATM_precip * ATM_fsic 1021 ATM_wemp_lic = ATM_wevap_lic - ATM_precip * ATM_flic 1022 ATM_wemp_sea = ATM_wevap_sic + ATM_wevap_oce - ATM_precip * ATM_fsea 1023 1024 print ( ' - Step 4', end='' ) 1025 if RUN_HIS == 'latlon' : 1026 RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell' ) 1027 RUN_riverflow = lmdz.geo2point ( d_RUN_his ['riverflow'] , dim1D='cell' ) 1028 RUN_runoff = lmdz.geo2point ( d_RUN_his ['runoff'] , dim1D='cell' ) 1029 RUN_drainage = lmdz.geo2point ( d_RUN_his ['drainage'] , dim1D='cell' ) 1030 RUN_riversret = lmdz.geo2point ( d_RUN_his ['riversret'] , dim1D='cell' ) 1031 1032 RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'], dim1D='cell' ) 1033 RUN_riverflow_cpl = lmdz.geo2point ( d_RUN_his ['riverflow_cpl'] , dim1D='cell' ) 1034 1035 print ( ' - Step 5', end='' ) 1036 if RUN_HIS == 'ico' : 1037 RUN_coastalflow = d_RUN_his ['coastalflow'][..., RUN_his_keysort] 1038 RUN_riverflow = d_RUN_his ['riverflow'] [..., RUN_his_keysort] 1039 RUN_runoff = d_RUN_his ['runoff'] [..., RUN_his_keysort] 1040 RUN_drainage = d_RUN_his ['drainage'] [..., RUN_his_keysort] 1041 RUN_riversret = d_RUN_his ['riversret'] [..., RUN_his_keysort] 1042 1043 RUN_coastalflow_cpl = d_RUN_his ['coastalflow_cpl'][..., RUN_his_keysort] 1044 RUN_riverflow_cpl = d_RUN_his ['riverflow_cpl'] [..., RUN_his_keysort] 1045 1046 print ( ' - Step 6', end='' ) 1047 if SRF_HIS == 'latlon' : 1048 SRF_rain = lmdz.geo2point ( d_SRF_his ['rain'] , dim1D='cell') 1049 SRF_evap = lmdz.geo2point ( d_SRF_his ['evap'] , dim1D='cell') 1050 SRF_snowf = lmdz.geo2point ( d_SRF_his ['snowf'] , dim1D='cell') 1051 SRF_subli = lmdz.geo2point ( d_SRF_his ['subli'] , dim1D='cell') 1052 SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell' ) 1053 1054 print ( '- Step 7', end='' ) 1055 if SRF_HIS == 'ico' : 1056 SRF_rain = d_SRF_his ['rain'] [..., SRF_his_keysort] 1057 SRF_evap = d_SRF_his ['evap'] [..., SRF_his_keysort] 1058 SRF_snowf = d_SRF_his ['snowf'][..., SRF_his_keysort] 1059 SRF_subli = d_SRF_his ['subli'][..., SRF_his_keysort] 1060 SRF_transpir = np.sum (d_SRF_his ['transpir'], axis=1)[..., SRF_his_keysort] 1061 1062 print ( '- Step 8', end='' ) 1063 SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 905 1064 SRF_emp = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 906 1065 907 1066 ## Correcting units of SECHIBA variables 908 1067 def mmd2SI ( Var ) : 909 1068 '''Change unit from mm/d or m^3/s to kg/s if needed''' 910 1069 if 'units' in VarT.attrs : 911 if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s- 3'] :1070 if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] : 912 1071 VarT.values = VarT.values * ATM_rho ; VarT.attrs['units'] = 'kg/s' 913 1072 if VarT.attrs['units'] == 'mm/d' : 914 1073 VarT.values = VarT.values * ATM_rho * (1e-3/86400.) ; VarT.attrs['units'] = 'kg/s' 915 1074 if VarT.attrs['units'] in ['m^3', 'm3', ] : 1075 VarT.values = VarT.values * ATM_rho ; VarT.attrs['units'] = 'kg' 1076 916 1077 for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 917 1078 VarT = locals()['RUN_' + var] … … 922 1083 mmd2SI (VarT) 923 1084 924 1085 print ( '- Step 9', end='' ) 925 1086 RUN_input = RUN_runoff + RUN_drainage 926 1087 RUN_output = RUN_coastalflow + RUN_riverflow 927 1088 928 ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic 929 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 )1089 print ( '- Step 10', end='' ) 1090 ATM_flx_wbilo = ATM_flux_int ( ATM_wbilo ) 1091 ATM_flx_wbilo_oce = ATM_flux_int ( ATM_wbilo_oce ) 1092 ATM_flx_wbilo_sic = ATM_flux_int ( ATM_wbilo_sic ) 1093 ATM_flx_wbilo_sea = ATM_flux_int ( ATM_wbilo_sea ) 1094 ATM_flx_wbilo_ter = ATM_flux_int ( ATM_wbilo_ter ) 1095 ATM_flx_wbilo_lic = ATM_flux_int ( ATM_wbilo_lic ) 935 1096 ATM_flx_calving = ATM_flux_int ( ATM_fqcalving ) 936 ATM_flx_qfonte = ATM_flux_int ( ATM_fqfonte ) 1097 ATM_flx_fqfonte = ATM_flux_int ( ATM_fqfonte ) 1098 1099 print ( '- Step 11', end='' ) 937 1100 ATM_flx_precip = ATM_flux_int ( ATM_precip ) 938 1101 ATM_flx_snowf = ATM_flux_int ( ATM_snowf ) … … 940 1103 ATM_flx_runlic = ATM_flux_int ( ATM_runofflic ) 941 1104 1105 print ( '- Step 12', end='' ) 1106 ATM_flx_wevap_ter = ATM_flux_int ( ATM_wevap_ter ) 1107 ATM_flx_wevap_sea = ATM_flux_int ( ATM_wevap_sea ) 1108 ATM_flx_precip_ter = ATM_flux_int ( ATM_precip * ATM_fter ) 1109 ATM_flx_precip_sea = ATM_flux_int ( ATM_precip * ATM_fsea ) 1110 ATM_flx_emp = ATM_flux_int ( ATM_emp) 1111 ATM_flx_wemp_ter = ATM_flux_int ( ATM_wemp_ter ) 1112 ATM_flx_wemp_oce = ATM_flux_int ( ATM_wemp_oce ) 1113 ATM_flx_wemp_lic = ATM_flux_int ( ATM_wemp_lic ) 1114 ATM_flx_wemp_sea = ATM_flux_int ( ATM_wemp_sea ) 1115 1116 print ( '- Step 13', end='' ) 942 1117 RUN_flx_coastal = ONE_flux_int ( RUN_coastalflow) 943 1118 RUN_flx_river = ONE_flux_int ( RUN_riverflow ) … … 950 1125 RUN_flx_output = ONE_flux_int ( RUN_output ) 951 1126 952 ATM_flx_emp = ATM_flux_int (ATM_emp) 953 ATM_flx_wbilo = ATM_flux_int (ATM_wbilo) 954 955 RUN_flx_bil = RUN_flx_input - RUN_flx_output 1127 print ( '- Step 14' ) 1128 RUN_flx_bil = RUN_flx_input - RUN_flx_output 956 1129 RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 957 1130 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' ) 1131 prtFlux ('wbilo_oce ', ATM_flx_wbilo_oce , 'f' ) 1132 prtFlux ('wbilo_sic ', ATM_flx_wbilo_sic , 'f' ) 1133 prtFlux ('wbilo_sic+oce ', ATM_flx_wbilo_sea , 'f' ) 1134 prtFlux ('wbilo_ter ', ATM_flx_wbilo_ter , 'f' ) 1135 prtFlux ('wbilo_lic ', ATM_flx_wbilo_lic , 'f' ) 964 1136 prtFlux ('Sum wbilo_* ', ATM_flx_wbilo , 'f', True) 965 1137 prtFlux ('E-P ', ATM_flx_emp , 'f', True) 966 1138 prtFlux ('calving ', ATM_flx_calving , 'f' ) 967 prtFlux ('fqfonte ', ATM_flx_ qfonte, 'f' )1139 prtFlux ('fqfonte ', ATM_flx_fqfonte , 'f' ) 968 1140 prtFlux ('precip ', ATM_flx_precip , 'f' ) 969 1141 prtFlux ('snowf ', ATM_flx_snowf , 'f' ) … … 982 1154 prtFlux ('runoff lic ', ATM_flx_runlic , 'f' ) 983 1155 984 ATM_flx_budget = -ATM_flx_ sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river1156 ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river 985 1157 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 ))1158 #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 1159 988 1160 #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 )) 1161 1162 ATM_flx_toSRF = -ATM_flx_wbilo_ter 989 1163 990 1164 echo (' ') 991 1165 echo ( '\n====================================================================================' ) 992 1166 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. ))1167 echo ( f'Mass begin = {DYN_mas_wat_beg:12.6e} kg | Mass end = {DYN_mas_wat_end:12.6e} kg' ) 1168 prtFlux ( 'dmass (atm) = ', dDYN_mas_wat , 'e', True ) 1169 prtFlux ( 'Sum wbilo_* = ', ATM_flx_wbilo, 'e', True ) 1170 prtFlux ( 'E-P = ', ATM_flx_emp , 'e', True ) 997 1171 echo ( ' ' ) 998 prtFlux ( ' Perte eau atm ', ATM_flx_wbilo - dDYN_mas_wat, 'e', True )999 echo ( ' Perte eauatm = {:12.3e} (rel) '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat ) )1172 prtFlux ( 'Water loss atm', ATM_flx_wbilo - dDYN_mas_wat, 'f', True ) 1173 echo ( 'Water loss atm = {:12.3e} (rel) '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat ) ) 1000 1174 1001 1175 echo (' ') 1002 prtFlux ( ' Perte eau atm', ATM_flx_emp - dDYN_mas_wat , 'e', True )1003 echo ( ' Perte eauatm = {:12.3e} (rel) '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat ) )1176 prtFlux ( 'Water loss atm', ATM_flx_emp - dDYN_mas_wat , 'f', True ) 1177 echo ( 'Water loss atm = {:12.3e} (rel) '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat ) ) 1004 1178 echo (' ') 1179 1180 echo (' ') 1181 echo ( '\n====================================================================================' ) 1182 1183 ATM_flx_runofflic_lic = ATM_flux_int ( ATM_runofflic * ATM_flic ) 1184 1185 LIC_flx_budget1 = -ATM_flx_wemp_lic - ATM_flx_calving - ATM_flx_fqfonte 1186 LIC_flx_budget2 = -ATM_flx_wbilo_lic - ATM_flx_calving - ATM_flx_fqfonte 1187 LIC_flx_budget3 = -ATM_flx_wbilo_lic - ATM_flx_runofflic_lic 1188 1189 echo ( f'-- LIC -- {Title} ' ) 1190 echo ( f'Mass total begin = {LIC_mas_wat_beg :12.6e} kg | Mass end = {LIC_mas_wat_end :12.6e} kg' ) 1191 echo ( f'Mass snow begin = {LIC_mas_sno_beg :12.6e} kg | Mass end = {LIC_mas_sno_end :12.6e} kg' ) 1192 echo ( f'Mass qs begin = {LIC_mas_qs_beg :12.6e} kg | Mass end = {LIC_mas_qs_end :12.6e} kg' ) 1193 echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' ) 1194 prtFlux ( 'dmass (LIC sno) ', dLIC_mas_sno , 'f', True, width=65 ) 1195 prtFlux ( 'dmass (LIC qs) ', dLIC_mas_qs , 'e', True, width=65 ) 1196 prtFlux ( 'dmass (LIC wat) ', dLIC_mas_wat , 'f', True, width=65 ) 1197 prtFlux ( 'dmass (LIC runlic0) ', dLIC_mas_runlic0 , 'e', True, width=65 ) 1198 prtFlux ( 'dmass (LIC total) ', dLIC_mas_wat , 'e', True, width=65 ) 1199 prtFlux ( 'LIC ATM_flx_wemp_lic ', ATM_flx_wemp_lic , 'f', True, width=65 ) 1200 prtFlux ( 'LIC ATM_flx_fqfonte ', ATM_flx_fqfonte , 'f', True, width=65 ) 1201 prtFlux ( 'LIC ATM_flx_calving ', ATM_flx_calving , 'f', True, width=65 ) 1202 prtFlux ( 'LIC ATM_flx_runofflic ', ATM_flx_runofflic_lic , 'f', True, width=65 ) 1203 prtFlux ( 'LIC fqfonte + calving ', ATM_flx_calving+ATM_flx_fqfonte , 'f', True, width=65 ) 1204 prtFlux ( 'LIC fluxes 1 (wevap - precip*frac_lic - fqcalving - fqfonte) ', LIC_flx_budget1 , 'f', True, width=65 ) 1205 prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte) ', LIC_flx_budget2 , 'f', True, width=65 ) 1206 prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic) ', LIC_flx_budget3 , 'f', True, width=65 ) 1207 prtFlux ( 'LIC error 1 ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=65 ) 1208 prtFlux ( 'LIC error 2 ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=65 ) 1209 prtFlux ( 'LIC error 3 ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=65 ) 1210 echo ( 'LIC error (wevap - precip*frac_lic - fqcalving - fqfonte) = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) ) 1211 echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte) = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) ) 1212 echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic) = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) ) 1005 1213 1006 1214 echo ( '\n====================================================================================' ) … … 1019 1227 SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage 1020 1228 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 )1229 prtFlux ('rain ', SRF_flx_rain , 'f' ) 1230 prtFlux ('evap ', SRF_flx_evap , 'f' ) 1231 prtFlux ('snowf ', SRF_flx_snowf , 'f' ) 1232 prtFlux ('E-P ', SRF_flx_emp , 'f' ) 1233 prtFlux ('subli ', SRF_flx_subli , 'f' ) 1234 prtFlux ('transpir ', SRF_flx_transpir , 'f' ) 1235 prtFlux ('to routing ', RUN_flx_torouting , 'f' ) 1236 prtFlux ('budget ', SRF_flx_all , 'f', small=True ) 1029 1237 1030 1238 echo ( '\n------------------------------------------------------------------------------------' ) 1031 1239 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) 1240 echo ( f'SRF_mas_wat_beg = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end = {SRF_mas_wat_end:12.6e} kg ' ) 1241 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True) 1242 prtFlux ( 'Error ', SRF_flx_all-dSRF_mas_wat, 'e', small=True ) 1034 1243 echo ( 'dMass (water srf) = {:12.4e} (rel) '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 1244 1245 echo ( '\n====================================================================================' ) 1246 echo ( f'-- Check ATM vs. SRF -- {Title} ' ) 1247 prtFlux ('E-P ATM ', ATM_flx_wemp_ter , 'f' ) 1248 prtFlux ('wbilo ter ', ATM_flx_wbilo_ter , 'f' ) 1249 prtFlux ('E-P SRF ', SRF_flx_emp , 'f' ) 1250 prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True) 1251 echo ( 'SRF/ATM error {:12.3e} (rel) '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) ) 1035 1252 1036 1253 echo ('') … … 1044 1261 prtFlux ('coastal ', RUN_flx_coastal , 'f' ) 1045 1262 prtFlux ('riv+coa ', RUN_flx_fromrouting , 'f' ) 1046 prtFlux ('budget ', RUN_flx_all , 'f' , True) 1263 prtFlux ('budget ', RUN_flx_all , 'f' , small=True) 1264 1265 echo ( '\n------------------------------------------------------------------------------------' ) 1266 echo ( f'Water content in routing+lake -- {Title} ' ) 1267 echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 1268 prtFlux ( 'dMass (routing) ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True) 1269 prtFlux ( 'Routing error ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True ) 1270 echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) ) 1047 1271 1048 1272 echo ( '\n------------------------------------------------------------------------------------' ) 1049 1273 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 ) ) 1274 echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 1275 prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True ) 1276 prtFlux ( 'Routing error ', RUN_flx_all-dRUN_mas_wat, 'e', small=True) 1277 echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 1059 1278 1060 1279 # PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &a 1061 1280 # 1150 " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" 1281 1282 echo ( ' ' ) 1283 echo ( f'{Title = }' ) -
TOOLS/WATER_BUDGET/CPL_waterbudget.py
r6277 r6508 22 22 ### 23 23 ## Import system modules 24 import sys, os, shutil, subprocess 24 import sys, os, shutil, subprocess, platform 25 25 import numpy as np 26 26 import configparser, re 27 28 ## Creates parser 29 config = configparser.ConfigParser() 27 from pathlib import Path 28 29 # Check python version 30 if sys.version_info < (3, 8, 0) : 31 print ( f'Python version : {platform.python_version()}' ) 32 raise Exception ( "Minimum Python version is 3.8" ) 33 34 ## Import local module 35 import WaterUtils as wu 36 37 ## Creates parser for reading .ini input file 38 config = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() ) 30 39 config.optionxform = str # To keep capitals 31 40 32 config['Files'] = {} 33 config['System'] = {} 34 35 ##-- Some physical constants 36 #-- Earth Radius 37 Ra = 6366197.7236758135 38 #-- Gravity 39 grav = 9.81 40 #-- Ice volumic mass (kg/m3) in LIM3/SI3 41 ICE_rho_ice = 917.0 42 #-- Snow volumic mass (kg/m3) in LIM3/SI3 43 ICE_rho_sno = 330.0 44 #-- Water in ponds coluclic mass (kg/m3) in LIM3/SI3 45 ICE_rho_pnd = 1000. 46 #-- Ocean water volumic mass (kg/m3) in NEMO 47 OCE_rho_liq = 1026. 48 #-- Water volumic mass in atmosphere 49 ATM_rho = 1.0e3 50 #-- Water volumic mass in surface reservoirs 51 SRF_rho = 1.0e3 52 #-- Water volumic mass of rivers 53 RUN_rho = 1.0e3 54 55 ## Read experiment parameters 56 ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 57 58 # Arguments passed 41 ## Experiment parameters 42 ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; 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 43 44 ## 45 ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 46 TmpDir=None ; RunDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 47 file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ; file_ICE_his=None ; file_OCE_sca=None 48 tar_restart_beg=None ; tar_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 49 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 ; file_OCE_srf=None ; file_OCE_sca=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None 50 51 # Read command line arguments 59 52 print ( "Name of Python script:", sys.argv[0] ) 60 IniFile = sys.argv[1] 53 IniFile = sys.argv[1] 54 # Text existence of IniFile 55 if not os.path.exists (IniFile ) : 56 raise FileExistsError ( f"File not found : {IniFile = }" ) 57 58 if 'full' in IniFile : FullIniFile = IniFile 59 else : FullIniFile = 'full_' + IniFile 60 61 61 print ("Input file : ", IniFile ) 62 62 config.read (IniFile) 63 64 def setBool (chars) : 65 '''Convert specific char string in boolean if possible''' 66 setBool = chars 67 for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : 68 if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] 69 return setBool 70 71 def setNum (chars) : 72 '''Convert specific char string in integer or real if possible''' 73 if type (chars) == str : 74 realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") 75 isReal = realnum.match(chars.strip()) != None 76 isInt = chars.strip().isdigit() 77 if isReal : 78 if isInt : setNum = int (chars) 79 else : setNum = float (chars) 80 else : setNum = chars 81 else : setNum = chars 82 return setNum 83 84 print ('[Experiment]') 85 for VarName in config['Experiment'].keys() : 86 locals()[VarName] = config['Experiment'][VarName] 87 locals()[VarName] = setBool (locals()[VarName]) 88 locals()[VarName] = setNum (locals()[VarName]) 89 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 90 91 ### 92 ICO = ( 'ICO' in ATM ) 93 LMDZ = ( 'LMD' in ATM ) 94 95 ### 96 ## Import system modules 97 import sys, os, shutil, subprocess 98 99 # Where do we run ? 100 SysName, NodeName, Release, Version, Machine = os.uname() 101 TGCC = ( 'irene' in NodeName ) 102 IDRIS = ( 'jeanzay' in NodeName ) 103 104 ## Set site specific libIGCM directories, and other specific stuff 105 if TGCC : 106 CPU = subprocess.getoutput ( 'lscpu | grep "Model name"' ) 107 if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene' 108 if "AMD" in CPU : Machine = 'irene-amd' 109 110 ARCHIVE = subprocess.getoutput ( f'ccc_home --cccstore -d {Project} -u {User}' ) 111 STORAGE = subprocess.getoutput ( f'ccc_home --cccwork -d {Project} -u {User}' ) 112 SCRATCHDIR = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 113 R_IN = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 114 rebuild = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 115 116 ## Specific to run at TGCC. 117 # Needed before importing a NetCDF library (netCDF4, xarray, cmds, etc ...) 118 import mpi4py 119 mpi4py.rc.initialize = False 120 121 ## Creates output directory 122 #TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 123 TmpDir = os.path.join ( '/ccc/scratch/cont003/drf/p86mart', f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 124 125 if IDRIS : 126 raise Exception("Pour IDRIS : repertoires et chemins a definir") 127 128 ## Import specific module 63 FullIniFile = 'full_' + IniFile 64 65 ## Reading config file 66 for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 67 if Section in config.keys () : 68 print ( f'Reading [{Section}]' ) 69 for VarName in config[Section].keys() : 70 locals()[VarName] = config[Section][VarName] 71 exec ( f'{VarName} = wu.setBool ({VarName})' ) 72 exec ( f'{VarName} = wu.setNum ({VarName})' ) 73 exec ( f'{VarName} = wu.setNone ({VarName})' ) 74 exec ( f'wu.{VarName} = {VarName}' ) 75 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 76 #exec ( f'del {VarName}' ) 77 78 #-- Some physical constants 79 if wu.unDefined ( 'Ra' ) : Ra = 6366197.7236758135 #-- Earth Radius (m) 80 if wu.unDefined ( 'Grav' ) : Grav = 9.81 #-- Gravity (m^2/s 81 if wu.unDefined ( 'ICE_rho_ice' ) : ICE_rho_ice = 917.0 #-- Ice volumic mass (kg/m3) in LIM3 82 if wu.unDefined ( 'ICE_rho_sno') : ICE_rho_sno = 330.0 #-- Snow volumic mass (kg/m3) in LIM3 83 if wu.unDefined ( 'OCE_rho_liq' ) : OCE_rho_liq = 1026.0 #-- Ocean water volumic mass (kg/m3) in NEMO 84 if wu.unDefined ( 'ATM_rho' ) : ATM_rho = 1000.0 #-- Water volumic mass in atmosphere (kg/m^3) 85 if wu.unDefined ( 'SRF_rho' ) : SRF_rho = 1000.0 #-- Water volumic mass in surface reservoir (kg/m^3) 86 if wu.unDefined ( 'RUN_rho' ) : RUN_rho = 1000.0 #-- Water volumic mass of rivers (kg/m^3) 87 if wu.unDefined ( 'ICE_rho_pnd' ) : ICE_rho_pnd = 1000. #-- Water volumic mass in ice ponds (kg/m^3) 88 if wu.unDefined ( 'YearLength' ) : YearLength = 365.25 * 24. * 60. * 60. #-- Year length (s) - Use only to compute drif in approximate unit 89 90 # Set libIGCM and machine dependant values 91 if not 'Files' in config.keys() : config['Files'] = {} 92 93 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} 94 95 ## -------------------------- 96 ICO = ( 'ICO' in wu.ATM ) 97 LMDZ = ( 'LMD' in wu.ATM ) 98 99 with open ('SetLibIGCM.py') as f: exec ( f.read() ) 100 config['Files']['TmpDir'] = TmpDir 101 102 if libIGCM : 103 config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 104 105 # Import specific module 129 106 import nemo, lmdz 130 107 ## Now import needed scientific modules 131 108 import xarray as xr 132 109 133 # Output file 134 FileOut = f'CPL_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 110 # Output file with water budget diagnostics 111 if wu.unDefined ( 'FileOut' ) : FileOut = f'CPL_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 112 config['Files']['FileOut'] = FileOut 113 135 114 f_out = open ( FileOut, mode = 'w' ) 136 137 # Function to print to stdout *and* output file 115 116 # Useful functions 117 def kg2Sv (val, rho=ATM_rho) : 118 '''From kg to Sverdrup''' 119 return val/dtime_sec*1.0e-6/rho 120 121 def kg2myear (val, rho=ATM_rho) : 122 '''From kg to m/year''' 123 return val/ATM_aire_sea_tot/rho/NbYear 124 125 def var2prt (var, small=False, rho=ATM_rho) : 126 if small : return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000. 127 else : return var , kg2Sv(var, rho=rho) , kg2myear(var, rho=rho) 128 129 def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) : 130 if small : 131 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 132 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 133 else : 134 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 135 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 136 echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) ) 137 return None 138 138 139 def echo (string, end='\n') : 139 print ( string, end=end ) 140 '''Function to print to stdout *and* output file''' 141 print ( str(string), end=end ) 140 142 sys.stdout.flush () 141 f_out.write ( str ing+ end )143 f_out.write ( str(string) + end ) 142 144 f_out.flush () 143 145 return None 144 146 145 147 ## Set libIGCM directories 146 R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT') 147 R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT') 148 149 L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 150 R_SAVE = os.path.join ( R_OUT, L_EXP ) 151 R_BUFR = os.path.join ( R_BUF, L_EXP ) 152 POST_DIR = os.path.join ( R_BUF, L_EXP, 'Out' ) 153 REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' ) 154 R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 155 R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 156 157 #if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir ) 158 if not os.path.isdir (TmpDir) : os.mkdir (TmpDir) 159 TmpDirOCE = os.path.join (TmpDir, 'OCE') 160 TmpDirICE = os.path.join (TmpDir, 'ICE') 161 if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE ) 162 if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 163 164 echo ( f'Working in TMPDIR : {TmpDir}' ) 165 166 echo ( f'\nDealing with {L_EXP}' ) 148 if wu.unDefined ('R_OUT' ) : R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT' ) 149 if wu.unDefined ('R_BUF' ) : R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT' ) 150 if wu.unDefined ('L_EXP' ) : L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 151 if wu.unDefined ('R_SAVE' ) : R_SAVE = os.path.join ( R_OUT, L_EXP ) 152 if wu.unDefined ('R_BUFR' ) : R_BUFR = os.path.join ( R_BUF, L_EXP ) 153 if wu.unDefined ('POST_DIR' ) : POST_DIR = os.path.join ( R_BUFR, 'Out' ) 154 if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' ) 155 if wu.unDefined ('R_BUF_KSH' ) : R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 156 if wu.unDefined ('R_FIGR' ) : R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 157 158 config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR } 159 160 # Set directory to extract filesa 161 if RunDir == None : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 162 config['Files']['RunDir'] = RunDir 163 164 if not os.path.isdir ( RunDir ) : os.makedirs ( RunDir ) 165 166 # Set directories to rebuild ocean and ice restart files 167 RunDirOCE = os.path.join ( RunDir, 'OCE' ) 168 RunDirICE = os.path.join ( RunDir, 'ICE' ) 169 #RunDirATM = os.path.join ( RunDir, 'ATM' ) 170 #RunDirSRF = os.path.join ( RunDir, 'SRF' ) 171 #RunDirRUN = os.path.join ( RunDir, 'RUN' ) 172 #RunDirDYN = os.path.join ( RunDir, 'DYN' ) 173 if not os.path.exists ( RunDirOCE ) : os.mkdir ( RunDirOCE ) 174 if not os.path.exists ( RunDirICE ) : os.mkdir ( RunDirICE ) 175 176 echo (' ') 177 echo ( f'JobName : {JobName}' ) 178 echo ( f'Comment : {Comment}' ) 179 echo ( f'TmpDir : {TmpDir}' ) 180 echo ( f'RunDirOCE : {RunDirOCE}' ) 181 echo ( f'RunDirICE : {RunDirICE}' ) 182 183 echo ( f'\nDealing with {L_EXP}' ) 167 184 168 185 #-- Model output directories 169 186 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 170 187 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 171 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 172 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 173 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 174 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 188 if dir_ATM_his == None : 189 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 190 config['Files']['dir_ATM_his'] = dir_ATM_his 191 if dir_SRF_his == None : 192 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 193 config['Files']['dir_SRF_his'] = dir_SRF_his 194 if dir_OCE_his == None : 195 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 196 config['Files']['dir_OCE_his'] = dir_OCE_his 197 if dir_ICE_his == None : 198 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 199 config['Files']['dir_ICE_his'] = dir_ICE_his 175 200 176 201 echo ( f'The analysis relies on files from the following model output directories : ' ) … … 180 205 echo ( f'{dir_SRF_his}' ) 181 206 182 #-- Files Names 183 if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 184 if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 207 #-- Creates files names 208 if Period == None : 209 if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 210 if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 211 config['Files']['Period'] = Period 212 if FileCommon == None : 213 FileCommon = f'{JobName}_{Period}' 214 config['Files']['FileCommon'] = FileCommon 215 216 if Title == None : 217 Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 218 config['Files']['Title'] = Title 185 219 186 220 echo ('\nOpen history files' ) 187 file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 188 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 189 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 190 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 191 file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 192 file_OCE_srf = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 221 if file_ATM_his == None : 222 file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 223 config['Files']['file_ATM_his'] = file_ATM_his 224 if file_SRF_his == None : 225 file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 226 config['Files']['file_SRF_his'] = file_SRF_his 227 #if Routing == 'SECHIBA' : 228 # file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 229 if Routing == 'SIMPLE' : 230 if file_RUN_his == None : 231 file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 232 config['Files']['file_RUN_his'] = file_RUN_his 233 if file_OCE_his == None : 234 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 235 config['Files']['file_OCE_his'] = file_OCE_his 236 if file_OCE_sca == None : 237 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 238 config['Files']['file_OCE_sca'] = file_OCE_sca 239 if file_ICE_his == None : 240 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 241 config['Files']['file_ICE_his'] = file_ICE_his 242 if file_OCE_srf == None : 243 file_OCE_srf = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 244 config['Files']['file_OCE_srf'] = file_OCE_srf 193 245 194 246 d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() … … 199 251 d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 200 252 d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 201 202 echo ( file_ATM_his ) 203 echo ( file_OCE_his ) 204 echo ( file_OCE_sca ) 205 echo ( file_ICE_his ) 206 echo ( file_SRF_his ) 253 if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 254 if Routing == 'SIMPLE' : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 255 256 257 echo ( f'{file_ATM_his = }' ) 258 echo ( f'{file_SRF_his = }' ) 259 if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' ) 260 echo ( f'{file_OCE_his = }' ) 261 echo ( f'{file_ICE_his = }' ) 262 echo ( f'{file_OCE_sca = }' ) 263 echo ( f'{file_OCE_srf = }' ) 207 264 208 265 ## Compute run length … … 216 273 dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds 217 274 dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 275 dtime_per_sec.attrs['unit'] = 's' 276 277 # Number of years 278 NbYear = dtime_sec / YearLength 218 279 219 280 #-- Open restart files … … 221 282 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 222 283 284 config['Files']['YearPre'] = f'{YearBegin}' 285 223 286 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 224 287 225 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' ) 226 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' ) 227 228 echo ( f'{file_restart_beg}' ) 229 echo ( f'{file_restart_end}' ) 230 231 file_OCE_beg = f'OCE_{JobName}_{YearRes}1231_restart.nc' 232 file_OCE_end = f'OCE_{JobName}_{YearEnd}1231_restart.nc' 233 file_ICE_beg = f'ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 234 file_ICE_end = f'ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 235 236 echo ( f'{file_OCE_beg}' ) 237 echo ( f'{file_OCE_end}' ) 238 echo ( f'{file_ICE_beg}' ) 239 echo ( f'{file_ICE_end}' ) 240 241 file_ATM_beg = f'ATM_{JobName}_{YearRes}1231_restartphy.nc' 242 file_ATM_end = f'ATM_{JobName}_{YearEnd}1231_restartphy.nc' 243 if LMDZ : 244 file_DYN_beg = f'ATM_{JobName}_{YearRes}1231_restart.nc' 245 file_DYN_end = f'ATM_{JobName}_{YearEnd}1231_restart.nc' 246 if ICO : 247 file_DYN_beg = f'ICO_{JobName}_{YearRes}1231_restart.nc' 248 file_DYN_end = f'ICO_{JobName}_{YearEnd}1231_restart.nc' 249 file_SRF_beg = f'SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 250 file_SRF_end = f'SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 251 liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg, ] 252 liste_end = [file_ATM_end, file_DYN_end, file_SRF_end, ] 253 echo ( f'{file_ATM_beg}') 254 echo ( f'{file_ATM_end}') 255 echo ( f'{file_SRF_beg}') 256 echo ( f'{file_SRF_end}') 257 258 if Routing == 'SIMPLE' : 259 file_RUN_beg = f'SRF_{JobName}_{YearRes}1231_routing_restart.nc' 260 file_RUN_end = f'SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 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 tar_restart_beg == None : 300 tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 301 config['Files']['tar_restart_beg'] = tar_restart_beg 302 if tar_restart_end == None : 303 tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 304 config['Files']['tar_restart_end'] = tar_restart_end 305 306 echo ( f'{tar_restart_beg}' ) 307 echo ( f'{tar_restart_end}' ) 308 309 if file_ATM_beg == None : 310 file_ATM_beg = f'{RunDir}/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'{RunDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc' 314 config['Files']['file_ATM_end'] = file_ATM_end 315 316 if file_DYN_beg == None : 317 if LMDZ : file_DYN_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restart.nc' 318 if ICO : file_DYN_beg = f'{RunDir}/ICO_{JobName}_{YearRes}1231_restart.nc' 319 config['Files']['file_DYN_beg'] = file_DYN_beg 320 321 if file_DYN_end == None : 322 if LMDZ : file_DYN_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restart.nc' 323 if ICO : file_DYN_end = f'{RunDir}/ICO_{JobName}_{YearEnd}1231_restart.nc' 324 config['Files']['file_DYN_end'] = file_DYN_end 325 326 if file_SRF_beg == None : 327 file_SRF_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 328 config['Files']['file_SRF_beg'] = file_SRF_beg 329 if file_SRF_end == None : 330 file_SRF_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 331 config['Files']['file_SRF_end'] = file_SRF_end 332 333 if file_OCE_beg == None : 334 file_OCE_beg = f'{RunDir}/OCE_{JobName}_{YearRes}1231_restart.nc' 335 config['Files']['file_OCE_beg'] = file_OCE_beg 336 if file_OCE_end == None : 337 file_OCE_end = f'{RunDir}/OCE_{JobName}_{YearEnd}1231_restart.nc' 338 config['Files']['file_OCE_end'] = file_OCE_end 339 if file_ICE_beg == None : 340 file_ICE_beg = f'{RunDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 341 config['Files']['file_ICE_beg'] = file_ICE_beg 342 if file_ICE_end == None : 343 file_ICE_end = f'{RunDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 344 config['Files']['file_ICE_end'] = file_ICE_end 345 346 liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg] 347 liste_end = [file_ATM_end, file_DYN_end, file_SRF_end] 348 349 if Routing == 'SIMPLE' : 350 if file_RUN_beg == None : 351 file_RUN_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc' 352 config['Files']['file_RUN_beg'] = file_RUN_beg 353 if file_RUN_end == None : 354 file_RUN_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 355 config['Files']['file_RUN_end'] = file_RUN_end 356 261 357 liste_beg.append ( file_RUN_beg ) 262 358 liste_end.append ( file_RUN_end ) 263 echo ( f'{file_RUN_beg}') 264 echo ( f'{file_RUN_end}') 359 echo ( f'{file_RUN_beg = }' ) 360 echo ( f'{file_RUN_end = }' ) 361 362 echo ( f'{file_ATM_beg = }' ) 363 echo ( f'{file_ATM_end = }' ) 364 echo ( f'{file_DYN_beg = }' ) 365 echo ( f'{file_DYN_end = }' ) 366 echo ( f'{file_SRF_beg = }' ) 367 echo ( f'{file_SRF_end = }' ) 368 echo ( f'{file_RUN_beg = }' ) 369 echo ( f'{file_RUN_end = }' ) 370 echo ( f'{file_OCE_beg = }' ) 371 echo ( f'{file_OCE_end = }' ) 372 echo ( f'{file_ICE_beg = }' ) 373 echo ( f'{file_ICE_end = }' ) 265 374 266 375 echo ('\nExtract restart files from tar : ATM, ICO and SRF') 267 for resFile in liste_beg : 268 if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 269 command = f'cd {TmpDir} ; tar xf {file_restart_beg} {resFile}' 376 for resFile in liste_beg + liste_end : 377 if os.path.exists ( os.path.join (RunDir, resFile) ) : 378 echo ( f'file found : {resFile = }' ) 379 else : 380 base_file = Path (file_name).stem # basename, and remove suffix 381 command = f'cd {RunDir} ; tar xf {tar_restart_beg} {base_resFile}.nc' 270 382 echo ( command ) 271 os.system ( command ) 272 273 for resFile in liste_end : 274 if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 275 command = f'cd {TmpDir} ; tar xf {file_restart_end} {resFile}' 276 echo ( command ) 277 os.system ( command ) 383 ierr = os.system ( command ) 384 if ierr == 0 : echo ( f'tar done : {base_resFile}') 385 else : raise Exception ( f'command failed : {command}' ) 278 386 279 387 echo ('\nOpening ATM SRF and ICO restart files') 280 d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze() 281 d_ATM_end = xr.open_dataset ( os.path.join (TmpDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze() 282 d_SRF_beg = xr.open_dataset ( os.path.join (TmpDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze() 283 d_SRF_end = xr.open_dataset ( os.path.join (TmpDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze() 284 d_DYN_beg = xr.open_dataset ( os.path.join (TmpDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze() 285 d_DYN_end = xr.open_dataset ( os.path.join (TmpDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze() 286 287 for var in d_SRF_beg.variables : 288 d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 289 d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 290 291 if ICO : 292 d_RUN_beg = xr.open_dataset ( os.path.join (TmpDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze() 293 d_RUN_end = xr.open_dataset ( os.path.join (TmpDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze() 388 d_ATM_beg = xr.open_dataset ( os.path.join (RunDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze() 389 d_ATM_end = xr.open_dataset ( os.path.join (RunDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze() 390 d_SRF_beg = xr.open_dataset ( os.path.join (RunDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze() 391 d_SRF_end = xr.open_dataset ( os.path.join (RunDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze() 392 d_DYN_beg = xr.open_dataset ( os.path.join (RunDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze() 393 d_DYN_end = xr.open_dataset ( os.path.join (RunDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze() 294 394 295 395 echo ('\nExtract and rebuild OCE and ICE restarts') … … 298 398 #ndomain_opa = d_zfile.attrs['DOMAIN_number_total'] 299 399 #d_zfile.close () 300 ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) .format()400 ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) #.format() 301 401 return int (ndomain_opa) 302 402 303 304 if not os.path.exists ( os.path.join (TmpDir, file_OCE_beg) ) : 305 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart_0000.nc') ): 306 command = f'cd {TmpDirOCE} ; tar xf {file_restart_beg} OCE_{JobName}_{YearRes}1231_restart_*.nc' 403 def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, RunDirComp=RunDirOCE ) : 404 '''Extract restart file from tar. Rebuild ocean files if needed''' 405 echo ( f'file to extract : {file_name} ' ) 406 if os.path.exists ( file_name ) : 407 echo ( f'-- File ready : {file_name}' ) 408 else : 409 echo ( f'-- Extracting {file_name}' ) 410 base_resFile = Path (file_name).stem # basename, and remove suffix 411 # Try to extract the rebuilded file 412 command = f'cd {RunDirComp} ; tar xf {tar_restart} {base_resFile}.nc' 307 413 echo ( command ) 308 os.system ( command ) 309 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') ) 310 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {file_OCE_beg} {TmpDir}' 311 echo ( command ) 312 os.system ( command ) 313 echo ( f'Rebuild done : {file_OCE_beg}') 314 315 if not os.path.exists ( os.path.join (TmpDir, file_OCE_end) ) : 316 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart_0000.nc') ): 317 command = f'cd {TmpDirOCE} ; tar xf {file_restart_end} OCE_{JobName}_{YearEnd}1231_restart_*.nc' 318 echo ( command ) 319 os.system ( command ) 320 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart') ) 321 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {file_OCE_end} {TmpDir}' 322 echo ( command ) 323 os.system ( command ) 324 echo ( f'Rebuild done : {file_OCE_end}') 325 326 if not os.path.exists ( os.path.join (TmpDir, file_ICE_beg) ) : 327 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ): 328 command = f'cd {TmpDirICE} ; tar xf {file_restart_beg} ICE_{JobName}_{YearRes}1231_restart_icemod_*.nc' 329 echo ( command ) 330 os.system ( command ) 331 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 332 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {file_ICE_beg} {TmpDir} ' 333 echo ( command ) 334 os.system ( command ) 335 echo ( f'Rebuild done : {file_OCE_beg}') 336 337 if not os.path.exists ( os.path.join (TmpDir, file_ICE_end) ) : 338 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearEnd}1231_restart_icemod') ): 339 command = f'cd {TmpDirICE} ; tar xf {file_restart_end} ICE_{JobName}_{YearEnd}1231_restart_icemod_*.nc' 340 echo ( command ) 341 os.system ( command ) 342 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 343 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {file_ICE_end} {TmpDir}' 344 echo ( command ) 345 os.system ( command ) 346 echo ( f'Rebuild done : {file_ICE_end}') 414 ierr = os.system ( command ) 415 if ierr == 0 : 416 echo ( f'tar done : {base_resFile}') 417 command = f'cd {RunDirComp} ; mv {base_resFile}.nc {RunDir}' 418 ierr = os.system ( command ) 419 if ierr == 0 : echo ( f'command done : {command}' ) 420 else : raise Exception ( f'command failed : {command}' ) 421 else : 422 if not os.path.exists ( os.path.join (RunDir, f'{base_resFile}_0000.nc') ): 423 command = f'cd {RunDirComp} ; tar xf {tar_restart_end} {base_resFile}_*.nc' 424 echo ( command ) 425 ierr = os.system ( command ) 426 if ierr == 0 : echo ( f'tar done : {file_OCE_beg}') 427 else : raise Exception ( f'command failed : {command}' ) 428 echo ('extract ndomain' ) 429 ndomain_opa = get_ndomain ( os.path.join (RunDir, f'{base_resFile}') ) 430 command = f'cd {RunDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {RunDir}' 431 echo ( command ) 432 ierr = os.system ( command ) 433 if ierr == 0 : echo ( f'Rebuild done : {base_resFile}.nc') 434 else : raise Exception ( f'command failed : {command}' ) 435 436 extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirOCE ) 437 extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, RunDirComp=RunDirOCE ) 438 extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirICE ) 439 extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, RunDirComp=RunDirICE ) 347 440 348 441 echo ('Opening OCE and ICE restart files') 349 442 if NEMO == 3.6 : 350 d_OCE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()351 d_OCE_end = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()352 d_ICE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()353 d_ICE_end = xr.open_dataset ( os.path.join ( TmpDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()443 d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 444 d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze() 445 d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze() 446 d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze() 354 447 if NEMO == 4.0 or NEMO == 4.2 : 355 d_OCE_beg = xr.open_dataset ( os.path.join (TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 356 d_OCE_end = xr.open_dataset ( os.path.join (TmpDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 357 d_ICE_beg = xr.open_dataset ( os.path.join (TmpDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 358 d_ICE_end = xr.open_dataset ( os.path.join (TmpDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 359 360 echo ( file_ATM_beg ) 361 echo ( file_ATM_end ) 362 echo ( file_DYN_beg ) 363 echo ( file_DYN_end ) 364 echo ( file_SRF_beg ) 365 echo ( file_SRF_end ) 366 if Routing == 'SIMPLE' : 367 echo ( file_RUN_beg ) 368 echo ( file_RUN_end ) 369 echo ( file_OCE_beg ) 370 echo ( file_OCE_end ) 371 echo ( file_ICE_beg ) 372 echo ( file_ICE_end ) 373 374 def ATM_stock_int (stock) : 375 '''Integrate stock on atmosphere grid''' 376 ATM_stock_int = np.sum ( np.sort ( (stock * DYN_aire).to_masked_array().ravel()) ) 377 return ATM_stock_int 378 379 def OCE_stock_int (stock) : 380 '''Integrate stock on ocean grid''' 381 OCE_stock_int = np.sum ( np.sort ( (stock * OCE_aire).to_masked_array().ravel()) ) 382 return OCE_stock_int 383 384 def ONE_stock_int (stock) : 385 '''Sum stock for field already integrated by cell''' 386 ONE_stock_int = np.sum ( np.sort ( (stock).to_masked_array().ravel()) ) 387 return ONE_stock_int 388 389 # ATM grid with cell surfaces 448 d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 449 d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 450 d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 451 d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 452 453 ## Write the full configuration 454 config_out = open (FullIniFile, 'w') 455 config.write (config_out ) 456 config_out.close () 457 458 for var in d_SRF_beg.variables : 459 d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 460 d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 461 462 if ICO : 463 d_RUN_beg = xr.open_dataset ( os.path.join (RunDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze() 464 d_RUN_end = xr.open_dataset ( os.path.join (RunDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze() 465 466 def kg2Sv (val, rho=OCE_rho_liq) : 467 '''From kg to Sverdrup''' 468 return val/dtime_sec*1.0e-6/rho 469 470 def kg2myear (val, rho=OCE_rho_liq) : 471 '''From kg to m/year''' 472 return val/OCE_aire_tot/rho/NbYear 473 474 def var2prt (var, small=False) : 475 if small : return var , kg2Sv(var)*1000., kg2myear(var)*1000. 476 else : return var , kg2Sv(var) , kg2myear(var) 477 478 def prtFlux (Desc, var, Form='F', small=False) : 479 if small : 480 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 481 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 482 else : 483 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 484 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 485 echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) ) 486 return None 487 488 489 390 490 # ATM grid with cell surfaces 391 491 if ICO : 392 jpja, jpia = d_ATM_his['aire'][0].shape 393 file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 394 d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 395 ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True ) 396 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 397 ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_lic'][0] ) 398 SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 399 492 if ATM_HIS == 'latlon' : 493 jpja, jpia = d_ATM_his['aire'][0].shape 494 file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 495 config['Files']['file_ATM_aire'] = file_ATM_aire 496 echo ( f'Aire sur grille reguliere : {file_ATM_aire = }' ) 497 d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 498 ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True ) 499 ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 500 ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] ) 501 ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] ) 502 ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] ) 503 SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 504 #SRF_aire = ATM_aire * lmdz.geo2point (d_SRF_his ['Contfrac'] ) 505 #SRF_aire = ATM_aire * ATM_fter 506 if ATM_HIS == 'ico' : 507 echo ( f'Aire sur grille icosaedre : {file_ATM_aire = }' ) 508 400 509 if LMDZ : 401 510 ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=True ) 402 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 403 ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 404 SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 511 ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 512 ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] ) 513 ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] ) 514 ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] ) 515 #SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 516 SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'] ) 517 518 ATM_fsea = ATM_foce + ATM_fsic 519 ATM_flnd = ATM_fter + ATM_flic 520 ATM_aire_fter = ATM_aire * ATM_fter 521 ATM_aire_flic = ATM_aire * ATM_flic 522 ATM_aire_fsic = ATM_aire * ATM_fsic 523 ATM_aire_foce = ATM_aire * ATM_foce 524 ATM_aire_flnd = ATM_aire * ATM_flnd 525 ATM_aire_fsea = ATM_aire * ATM_fsea 526 527 SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.) 405 528 406 529 if ICO : … … 415 538 416 539 if LMDZ : 540 # Area on lon/lat grid 417 541 DYN_aire = ATM_aire 418 542 DYN_fsea = ATM_fsea 419 543 DYN_flnd = ATM_flnd 420 544 DYN_fter = d_ATM_beg['FTER'] 545 DYN_flic = d_ATM_beg['FTER'] 546 547 def ATM_stock_int (stock) : 548 '''Integrate (* surface) stock on atmosphere grid''' 549 ATM_stock_int = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 550 return ATM_stock_int 551 552 def ATM_flux_int (flux) : 553 '''Integrate (* time * surface) flux on atmosphere grid''' 554 ATM_flux_int = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 555 return ATM_flux_int 556 557 def SRF_stock_int (stock) : 558 '''Integrate (* surface) stock on land grid''' 559 ATM_stock_int = wu.Ksum ( ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 560 return ATM_stock_int 561 562 def SRF_flux_int (flux) : 563 '''Integrate (* time * surface) flux on land grid''' 564 SRF_flux_int = wu.Psum ( (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 565 return SRF_flux_int 566 567 def OCE_stock_int (stock) : 568 '''Integrate stock on ocean grid''' 569 OCE_stock_int = np.sum ( np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) ) 570 return OCE_stock_int 571 572 def ONE_stock_int (stock) : 573 '''Sum stock''' 574 ONE_stock_int = np.sum ( np.sort ( (stock ).to_masked_array().ravel()) ) 575 return ONE_stock_int 576 577 def OCE_flux_int (flux) : 578 '''Integrate flux on oce grid''' 579 OCE_flux_int = np.sum ( np.sort ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel()) ) 580 return OCE_flux_int 581 582 def ONE_flux_int (flux) : 583 '''Integrate flux on oce grid''' 584 OCE_flux_int = np.sum ( np.sort ( (flux * dtime_per_sec).to_masked_array().ravel()) ) 585 return OCE_flux_int 586 421 587 #if LMDZ : 422 588 # d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} ) … … 440 606 ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea ) 441 607 442 echo ( '\n ------------------------------------------------------------------------------------' )608 echo ( '\n====================================================================================' ) 443 609 echo ( '-- NEMO change in stores (for the records)' ) 444 610 #-- Note that the total number of days of the run should be diagnosed rather than assumed … … 522 688 dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat 523 689 524 echo ( 'ICE_vol_ice_beg = {:12.6e} m^3 | ICE_vol_ice_end = {:12.6e} m^3'.format (ICE_vol_ice_beg, ICE_vol_ice_end))525 echo ( 'ICE_vol_sno_beg = {:12.6e} m^3 | ICE_vol_sno_end = {:12.6e} m^3'.format (ICE_vol_sno_beg, ICE_vol_sno_end))526 echo ( 'ICE_vol_pnd_beg = {:12.6e} m^3 | ICE_vol_pnd_end = {:12.6e} m^3'.format (ICE_vol_pnd_beg, ICE_vol_pnd_end))527 echo ( 'ICE_vol_fzl_beg = {:12.6e} m^3 | ICE_vol_fzl_end = {:12.6e} m^3'.format (ICE_vol_fzl_beg, ICE_vol_fzl_end))528 529 echo ( 'dICE_vol_ice = {:12.3e} m^3'.format (dICE_vol_ice))530 echo ( 'dICE_vol_sno = {:12.3e} m^3'.format (dICE_vol_sno))531 echo ( 'dICE_vol_pnd = {:12.3e} m^3'.format (dICE_vol_pnd))532 echo ( 'dICE_mas_ice = {:12.3e} m^3'.format (dICE_mas_ice))533 echo ( 'dICE_mas_sno = {:12.3e} m^3'.format (dICE_mas_sno))534 echo ( 'dICE_mas_pnd = {:12.3e} m^3'.format (dICE_mas_pnd))535 echo ( 'dICE_mas_fzl = {:12.3e} m^3'.format (dICE_mas_fzl))536 echo ( 'dICE_mas_wat = {:12.3e} m^3'.format (dICE_mas_wat))690 echo ( f'ICE_vol_ice_beg = {ICE_vol_ice_beg:12.6e} m^3 | ICE_vol_ice_end = {ICE_vol_ice_end:12.6e} m^3' ) 691 echo ( f'ICE_vol_sno_beg = {ICE_vol_sno_beg:12.6e} m^3 | ICE_vol_sno_end = {ICE_vol_sno_end:12.6e} m^3' ) 692 echo ( f'ICE_vol_pnd_beg = {ICE_vol_pnd_beg:12.6e} m^3 | ICE_vol_pnd_end = {ICE_vol_pnd_end:12.6e} m^3' ) 693 echo ( f'ICE_vol_fzl_beg = {ICE_vol_fzl_beg:12.6e} m^3 | ICE_vol_fzl_end = {ICE_vol_fzl_end:12.6e} m^3' ) 694 695 echo ( f'dICE_vol_ice = {dICE_vol_ice:12.3e} m^3' ) 696 echo ( f'dICE_vol_sno = {dICE_vol_sno:12.3e} m^3' ) 697 echo ( f'dICE_vol_pnd = {dICE_vol_pnd:12.3e} m^3' ) 698 echo ( f'dICE_mas_ice = {dICE_mas_ice:12.3e} m^3' ) 699 echo ( f'dICE_mas_sno = {dICE_mas_sno:12.3e} m^3' ) 700 echo ( f'dICE_mas_pnd = {dICE_mas_pnd:12.3e} m^3' ) 701 echo ( f'dICE_mas_fzl = {dICE_mas_fzl:12.3e} m^3' ) 702 echo ( f'dICE_mas_wat = {dICE_mas_wat:12.3e} m^3' ) 537 703 538 704 … … 576 742 577 743 for k in np.arange (klevp1-1) : 578 DYN_sigma_beg[k,:] = (DYN_p_beg[k,:] - DYN_p_beg[k+1,:]) / grav579 DYN_sigma_end[k,:] = (DYN_p_end[k,:] - DYN_p_end[k+1,:]) / grav744 DYN_sigma_beg[k,:] = (DYN_p_beg[k,:] - DYN_p_beg[k+1,:]) / Grav 745 DYN_sigma_end[k,:] = (DYN_p_end[k,:] - DYN_p_end[k+1,:]) / Grav 580 746 581 747 DYN_sigma_beg = DYN_sigma_beg.rename ( {'klevp1':'sigs'} ) … … 769 935 770 936 echo ('\nLes differents reservoirs') 771 echo ( 'SRF_mas_watveg_beg = {:12.6e} kg | SRF_mas_watveg_end = {:12.6e} kg '.format (SRF_mas_watveg_beg , SRF_mas_watveg_end))772 echo ( 'SRF_mas_watsoil_beg = {:12.6e} kg | SRF_mas_watsoil_end = {:12.6e} kg '.format (SRF_mas_watsoil_beg, SRF_mas_watsoil_end))773 echo ( 'SRF_mas_snow_beg = {:12.6e} kg | SRF_mas_snow_end = {:12.6e} kg '.format (SRF_mas_snow_beg , SRF_mas_snow_end))937 echo ( f'SRF_mas_watveg_beg = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end = {SRF_mas_watveg_end :12.6e} kg ' ) 938 echo ( f'SRF_mas_watsoil_beg = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end = {SRF_mas_watsoil_end:12.6e} kg ' ) 939 echo ( f'SRF_mas_snow_beg = {SRF_mas_snow_beg :12.6e} kg | SRF_mas_snow_end = {SRF_mas_snow_end :12.6e} kg ' ) 774 940 775 941 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 /OCE_aire_tot*1E-3) ) -
TOOLS/WATER_BUDGET/OCE_waterbudget.py
r6277 r6508 24 24 import numpy as np 25 25 import configparser, re 26 27 ## Creates parser 26 from pathlib import Path 27 28 ## Import local module 29 import WaterUtils as wu 30 31 # Check python version 32 if sys.version_info < (3, 8, 0) : 33 print ( f'Python version : {platform.python_version()}' ) 34 raise Exception ( "Minimum Python version is 3.8" ) 35 36 ## Creates parser for reading .ini input file 28 37 config = configparser.ConfigParser() 29 38 config.optionxform = str # To keep capitals … … 33 42 34 43 ## 35 ARCHIVE=None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild= None36 TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None44 ARCHIVE=None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild='rebuild_nemo' 45 TmpDir=None ; RunDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 37 46 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=None47 tar_restart_beg=None ; tar_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 48 file_RUN_beg=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None 40 41 # Arguments passed 49 ContinueOnError=False ; ErrorCount=0 50 tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None ; tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None 51 tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None ; tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None 52 53 # Read command line arguments 42 54 print ( "Name of Python script:", sys.argv[0] ) 43 55 IniFile = sys.argv[1] 56 # Text existence of IniFile 57 print ("Input file : ", IniFile ) 58 59 if 'full' in IniFile : FullIniFile = IniFile 60 else : FullIniFile = 'full_' + IniFile 61 44 62 print ("Input file : ", IniFile ) 45 63 config.read (IniFile) 46 if 'full' in IniFile : FullIniFile = IniFile 47 else : FullIniFile = 'full_' + IniFile 48 49 def setBool (chars) : 50 '''Convert specific char string in boolean if possible''' 51 setBool = chars 52 for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : 53 if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] 54 return setBool 55 56 def setNum (chars) : 57 '''Convert specific char string in integer or real if possible''' 58 if type (chars) == str : 59 realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") 60 isReal = realnum.match(chars.strip()) != None 61 isInt = chars.strip().isdigit() 62 if isReal : 63 if isInt : setNum = int (chars) 64 else : setNum = float (chars) 65 else : setNum = chars 66 else : setNum = chars 67 return setNum 68 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 64 FullIniFile = 'full_' + IniFile 65 66 ## Reading config file 78 67 for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 79 if Section in config.keys() :80 print ( f' [{Section}]' )68 if Section in config.keys () : 69 print ( f'Reading [{Section}]' ) 81 70 for VarName in config[Section].keys() : 82 71 locals()[VarName] = config[Section][VarName] 83 locals()[VarName] = setBool (locals()[VarName]) 84 locals()[VarName] = setNum (locals()[VarName]) 72 exec ( f'{VarName} = wu.setBool ({VarName})' ) 73 exec ( f'{VarName} = wu.setNum ({VarName})' ) 74 exec ( f'{VarName} = wu.setNone ({VarName})' ) 75 exec ( f'wu.{VarName} = {VarName}' ) 85 76 print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 86 77 #exec ( f'del {VarName}' ) 78 79 #-- Some physical constants 80 if wu.unDefined ( 'Ra' ) : Ra = 6366197.7236758135 #-- Earth Radius (m) 81 if wu.unDefined ( 'Grav' ) : Grav = 9.81 #-- Gravity (m^2/s 82 if wu.unDefined ( 'ICE_rho_ice' ) : ICE_rho_ice = 917.0 #-- Ice volumic mass (kg/m3) in LIM3 83 if wu.unDefined ( 'ICE_rho_sno') : ICE_rho_sno = 330.0 #-- Snow volumic mass (kg/m3) in LIM3 84 if wu.unDefined ( 'OCE_rho_liq' ) : OCE_rho_liq = 1026.0 #-- Ocean water volumic mass (kg/m3) in NEMO 85 if wu.unDefined ( 'ATM_rho' ) : ATM_rho = 1000.0 #-- Water volumic mass in atmosphere (kg/m^3) 86 if wu.unDefined ( 'SRF_rho' ) : SRF_rho = 1000.0 #-- Water volumic mass in surface reservoir (kg/m^3) 87 if wu.unDefined ( 'RUN_rho' ) : RUN_rho = 1000.0 #-- Water volumic mass of rivers (kg/m^3) 88 if wu.unDefined ( 'ICE_rho_pnd' ) : ICE_rho_pnd = 1000. #-- Water volumic mass in ice ponds (kg/m^3) 89 if wu.unDefined ( 'YearLength' ) : YearLength = 365.25 * 24. * 60. * 60. #-- Year length (s) 90 91 # Set libIGCM and machine dependant values 87 92 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} 93 94 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 95 120 121 122 # Where do we run ? 123 SysName, NodeName, Release, Version, Machine = os.uname() 124 TGCC = ( 'irene' in NodeName ) 125 IDRIS = ( 'jeanzay' in NodeName ) 126 127 ## Set site specific libIGCM directories, and other specific stuff 128 if TGCC : 129 CPU = subprocess.getoutput ( 'lscpu | grep "Model name"' ) 130 if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene' 131 if "AMD" in CPU : Machine = 'irene-amd' 132 133 ARCHIVE = subprocess.getoutput ( f'ccc_home --cccstore -d {Project} -u {User}' ) 134 STORAGE = subprocess.getoutput ( f'ccc_home --cccwork -d {Project} -u {User}' ) 135 SCRATCHDIR = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 136 R_IN = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg'), 'IGCM') 137 rebuild = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg'), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 138 139 ## Specific to run at TGCC. 140 # Needed before importing a NetCDF library (netCDF4, xarray, cmds, etc ...) 141 import mpi4py 142 mpi4py.rc.initialize = False 143 144 ## Creates output directory 145 TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 146 147 if IDRIS : 148 raise Exception ("Pour IDRIS : repertoires et chemins a definir") 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$"} 96 with open ('SetLibIGCM.py') as f: exec ( f.read() ) 97 config['Files']['TmpDir'] = TmpDir 160 98 161 99 if libIGCM : … … 167 105 import xarray as xr 168 106 169 # Output file 170 if FileOut == None : 171 FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 172 config['Files']['FileOut'] = FileOut 107 # Output file with water budget diagnostics 108 if wu.unDefined ( 'FileOut' ) : FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 109 config['Files']['FileOut'] = FileOut 173 110 174 111 f_out = open ( FileOut, mode = 'w' ) 175 176 # Function to print to stdout *and* output file 112 113 # Useful functions 114 def kg2Sv (val, rho=ATM_rho) : 115 '''From kg to Sverdrup''' 116 return val/dtime_sec*1.0e-6/rho 117 118 def kg2myear (val, rho=ATM_rho) : 119 '''From kg to m/year''' 120 return val/ATM_aire_sea_tot/rho/NbYear 121 122 def var2prt (var, small=False, rho=ATM_rho) : 123 if small : return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000. 124 else : return var , kg2Sv(var, rho=rho) , kg2myear(var, rho=rho) 125 126 def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) : 127 if small : 128 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 129 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 130 else : 131 if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv | {:12.4f} m/year " 132 if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv | {:12.4e} m/year " 133 echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) ) 134 return None 135 177 136 def echo (string, end='\n') : 137 '''Function to print to stdout *and* output file''' 178 138 print ( str(string), end=end ) 179 139 sys.stdout.flush () … … 182 142 return None 183 143 184 ## Set libIGCM directories 185 R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT') 186 R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT') 187 188 L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 189 R_SAVE = os.path.join ( R_OUT, L_EXP ) 190 R_BUFR = os.path.join ( R_BUF, L_EXP ) 191 POST_DIR = os.path.join ( R_BUF, L_EXP, 'Out' ) 192 REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' ) 193 R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 194 R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 195 196 #if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir ) 197 if not os.path.isdir (TmpDir) : os.mkdir (TmpDir) 198 TmpDirOCE = os.path.join (TmpDir, 'OCE') 199 TmpDirICE = os.path.join (TmpDir, 'ICE') 200 if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE ) 201 if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 144 ### Set libIGCM directories 145 if wu.unDefined ('R_OUT' ) : R_OUT = os.path.join ( ARCHIVE , 'IGCM_OUT' ) 146 if wu.unDefined ('R_BUF' ) : R_BUF = os.path.join ( SCRATCHDIR, 'IGCM_OUT' ) 147 if wu.unDefined ('L_EXP' ) : L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 148 if wu.unDefined ('R_SAVE' ) : R_SAVE = os.path.join ( R_OUT, L_EXP ) 149 if wu.unDefined ('R_BUFR' ) : R_BUFR = os.path.join ( R_BUF, L_EXP ) 150 if wu.unDefined ('POST_DIR' ) : POST_DIR = os.path.join ( R_BUFR, 'Out' ) 151 if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' ) 152 if wu.unDefined ('R_BUF_KSH' ) : R_BUF_KSH = os.path.join ( R_BUFR, 'Out' ) 153 if wu.unDefined ('R_FIGR' ) : R_FIGR = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 154 155 config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR } 156 157 # Set directory to extract files 158 if wu.unDefined ( 'RunDir' ) : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 159 config['Files']['RunDir'] = RunDir 160 161 if not os.path.isdir (RunDir) : os.makedirs (RunDir) 162 163 # Set directories to rebuild ocean and ice restart files 164 RunDirOCE = os.path.join (RunDir, 'OCE') 165 RunDirICE = os.path.join (RunDir, 'ICE') 166 if not os.path.exists (RunDirOCE) : os.mkdir (RunDirOCE ) 167 if not os.path.exists (RunDirICE) : os.mkdir (RunDirICE ) 202 168 203 169 echo (' ') 204 echo ( f'JobName : {JobName}' ) 205 echo (Comment) 206 echo ( f'Working in TMPDIR : {TmpDir}' ) 207 208 echo ( f'\nDealing with {L_EXP}' ) 209 210 #-- Model output directories 211 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 212 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 213 if dir_OCE_his == None : 170 echo ( f'JobName : {JobName}' ) 171 echo ( f'Comment : {Comment}' ) 172 echo ( f'TmpDir : {TmpDir}' ) 173 echo ( f'RunDirOCE : {RunDirOCE}' ) 174 echo ( f'RunDirICE : {RunDirICE}' ) 175 176 echo ( f'\nDealing with {L_EXP}' ) 177 178 #-- Creates model output directory names 179 if Freq == "MO" : FreqDir = os.path.join ( 'Output' , 'MO' ) 180 if Freq == "SE" : FreqDir = os.path.join ( 'Analyse', 'SE' ) 181 if wu.unDefined ( 'dir_OCE_his' ) : 214 182 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 215 183 config['Files']['dir_OCE_his'] = dir_OCE_his 216 if dir_ICE_his == None:184 if wu.unDefined ( 'dir_ICE_his' ) : 217 185 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 218 186 config['Files']['dir_OCE_his'] = dir_OCE_his … … 222 190 echo ( f'{dir_ICE_his}' ) 223 191 224 #-- Files Names225 if Period == None:192 #-- Creates files names 193 if wu.unDefined ( 'Period ' ) : 226 194 if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 227 195 if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 228 196 config['Files']['Period'] = Period 229 if FileCommon == None:197 if wu.unDefined ( 'FileCommon' ) : 230 198 FileCommon = f'{JobName}_{Period}' 231 199 config['Files']['FileCommon'] = FileCommon 232 200 233 if Title == None:201 if wu.unDefined ( 'Title' ) : 234 202 Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 235 203 config['Files']['Title'] = Title … … 237 205 238 206 echo ('\nOpen history files' ) 239 if file_OCE_his == None:207 if wu.unDefined ('file_OCE_his' ) : 240 208 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 241 209 file_OCE_his = file_OCE_his 242 if file_OCE_sca == None:210 if wu.unDefined ('file_OCE_sca' ) : 243 211 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 244 212 config['Files']['file_OCE_sca'] = file_OCE_sca 245 if file_ICE_his == None:213 if wu.unDefined ( 'file_ICE_hi' ) : 246 214 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 247 215 config['Files']['file_ICE_his'] = file_ICE_his … … 258 226 ## Compute run length 259 227 dtime = ( d_OCE_his.time_counter_bounds.max() - d_OCE_his.time_counter_bounds.min() ) 260 echo ( '\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ))228 echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' ) 261 229 dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds 262 230 263 231 ## Compute length of each period 264 232 dtime_per = ( d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] ) 265 echo ( '\nPeriods lengths (days) : ')266 echo ( ' {:}'.format ( (dtime_per/np.timedelta64(1, "D")).values ))233 echo ( f'\nPeriods lengths (days) : ') 234 echo ( f' {(dtime_per/np.timedelta64(1, "D")).values}' ) 267 235 dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds 268 236 dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_OCE_his.time_counter,] ) … … 271 239 ## Number of years 272 240 NbYear = dtime_sec / YearLength 241 273 242 #-- Open restart files 274 243 YearRes = YearBegin - 1 # Year of the restart of beginning of simulation … … 277 246 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 278 247 279 if TarRestartPeriod_beg == None:248 if wu.unDefined ( 'TarRestartPeriod_beg' ) : 280 249 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') 281 250 TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 282 251 config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 283 252 284 if TarRestartPeriod_end == None:253 if wu.unDefined ( 'TarRestartPeriod_end' ) : 285 254 YearPre = YearBegin - PackFrequency # Year to find the tarfile of the restart of beginning of simulation 286 255 echo (f'Restart dates - Start : {YearRes}-12-31 / End : {YearEnd}-12-31 ') … … 288 257 config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 289 258 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 296 297 echo ( f'{file_restart_beg}' ) 298 echo ( f'{file_restart_end}' ) 299 300 if file_OCE_beg == None : 301 file_OCE_beg = f'{TmpDir}/OCE_{JobName}_{YearRes}1231_restart.nc' 259 if wu.unDefined ( 'tar_restart_beg' ) : 260 tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 261 config['Files']['tar_restart_beg'] = tar_restart_beg 262 if wu.unDefined ( 'tar_restart_end' ) : 263 tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 264 config['Files']['tar_restart_end'] = tar_restart_end 265 266 echo ( f'{tar_restart_beg}' ) 267 echo ( f'{tar_restart_end}' ) 268 269 # 270 if wu.unDefined ('tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg 271 if wu.unDefined ('tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg 272 if wu.unDefined ('tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg 273 if wu.unDefined ('tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg 274 if wu.unDefined ('tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg 275 if wu.unDefined ('tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg 276 277 if wu.unDefined ('tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end 278 if wu.unDefined ('tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end 279 if wu.unDefined ('tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end 280 if wu.unDefined ('tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end 281 if wu.unDefined ('tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end 282 if wu.unDefined ('tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end 283 284 if wu.unDefined ( 'file_OCE_beg' ) : 285 file_OCE_beg = f'{RunDir}/OCE_{JobName}_{YearRes}1231_restart.nc' 302 286 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'287 if wu.unDefined ( 'file_OCE_end' ) : 288 file_OCE_end = f'{RunDir}/OCE_{JobName}_{YearEnd}1231_restart.nc' 305 289 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'290 if wu.unDefined ( 'file_ICE_beg' ) : 291 file_ICE_beg = f'{RunDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 308 292 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'293 if wu.unDefined ( 'file_ICE_end' ) : 294 file_ICE_end = f'{RunDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 311 295 config['Files']['file_ICE_end'] = file_ICE_end 312 296 … … 321 305 #ndomain_opa = d_zfile.attrs['DOMAIN_number_total'] 322 306 #d_zfile.close () 323 ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) .format()307 ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) #.format() 324 308 return int (ndomain_opa) 325 309 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) 329 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart_0000.nc') ) : 330 command = f'cd {TmpDirOCE} ; tar xf {file_restart_beg} OCE_{JobName}_{YearRes}1231_restart_*.nc' 331 echo ( command ) 332 os.system ( command ) 333 echo ('extract ndomain' ) 334 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') ) 335 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {base_file_OCE_beg} {TmpDir}' 336 echo ( command ) 337 os.system ( command ) 338 echo ( f'Rebuild done : {file_OCE_beg}') 339 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) 343 if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart_0000.nc') ): 344 command = f'cd {TmpDirOCE} ; tar xf {file_restart_end} OCE_{JobName}_{YearEnd}1231_restart_*.nc' 345 echo ( command ) 346 os.system ( command ) 347 echo ('extract ndomain' ) 348 ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart') ) 349 command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {base_file_OCE_end} {TmpDir}' 350 echo ( command ) 351 os.system ( command ) 352 echo ( f'Rebuild done : {file_OCE_end}') 353 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) 357 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod_0000.nc') ): 358 command = f'cd {TmpDirICE} ; tar xf {file_restart_beg} ICE_{JobName}_{YearRes}1231_restart_icemod_*.nc' 359 echo ( command ) 360 os.system ( command ) 361 echo ('extract ndomain' ) 362 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 363 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_beg} {TmpDir} ' 364 echo ( command ) 365 os.system ( command ) 366 echo ( f'Rebuild done : {file_OCE_beg}') 367 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) 371 if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearEnd}1231_restart_icemod_0000.nc') ): 372 command = f'cd {TmpDirICE} ; tar xf {file_restart_end} ICE_{JobName}_{YearEnd}1231_restart_icemod_*.nc' 373 echo ( command ) 374 os.system ( command ) 375 echo ('extract ndomain' ) 376 ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 377 command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_end} {TmpDir}' 378 echo ( command ) 379 os.system ( command ) 380 echo ( f'Rebuild done : {file_ICE_end}') 310 def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, RunDirComp=RunDirOCE, ErrorCount=ErrorCount ) : 311 '''Extract restart file from tar. Rebuild ocean files if needed''' 312 echo ( f'----------') 313 if os.path.exists ( file_name ) : 314 echo ( f'-- File ready : {file_name = }' ) 315 else : 316 echo ( f'-- Extracting {file_name = }' ) 317 base_resFile = Path (file_name).stem # basename, and remove suffix 318 # Try to extract the rebuilded file 319 if os.path.exists ( tar_restart ) : 320 command = f'cd {RunDirComp} ; tar xf {tar_restart} {base_resFile}.nc' 321 echo ( command ) 322 try : 323 os.system ( command ) 324 except : 325 if not os.path.exists ( os.path.join (RunDir, f'{base_resFile}_0000.nc') ): 326 command = f'cd {RunDirComp} ; tar xf {tar_restart_end} {base_file}_*.nc' 327 echo ( command ) 328 ierr = os.system ( command ) 329 if ierr == 0 : echo ( f'tar done : {base_resFile}.nc') 330 else : raise Exception ( f'command failed : {command}' ) 331 echo ( f'extract ndomain' ) 332 ndomain_opa = get_ndomain ( os.path.join (RunDir, f'{base_file}') ) 333 command = f'cd {RunDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {RunDir}' 334 echo ( command ) 335 ierr = os.system ( command ) 336 if ierr == 0 : echo ( f'Rebuild done : {base_resFile}.nc') 337 else : raise Exception ( f'command failed : {command}' ) 338 else : 339 echo ( f'tar done : {base_resFile}') 340 command = f'cd {RunDirComp} ; mv {base_resFile}.nc {RunDir}' 341 ierr = os.system ( command ) 342 if ierr == 0 : echo ( f'command done : {command}' ) 343 else : raise Exception ( f'command failed : {command = }' ) 344 else : 345 echo ( f'****** Tar restart file {tar_restart = } not found ' ) 346 if ContinueOnError : 347 ErrorCount += 1 348 else : 349 raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' ) 350 return ErrorCount 351 352 353 ErrorCount += extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirOCE ) 354 ErrorCount += extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, RunDirComp=RunDirOCE ) 355 ErrorCount += extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirICE ) 356 ErrorCount += extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, RunDirComp=RunDirICE ) 381 357 382 358 echo ('Opening OCE and ICE restart files') 383 359 if NEMO == 3.6 : 384 d_OCE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()385 d_OCE_end = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()386 d_ICE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()387 d_ICE_end = xr.open_dataset ( os.path.join ( TmpDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()360 d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 361 d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze() 362 d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze() 363 d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze() 388 364 if NEMO == 4.0 or NEMO == 4.2 : 389 d_OCE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()390 d_OCE_end = xr.open_dataset ( os.path.join ( TmpDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()391 d_ICE_beg = xr.open_dataset ( os.path.join ( TmpDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()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()393 394 ## 365 d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 366 d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 367 d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 368 d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 369 370 ## Write the full configuration 395 371 config_out = open (FullIniFile, 'w') 396 372 config.write (config_out ) … … 473 449 OCE_sum_e3tn_end = OCE_stock_int ( OCE_e3tn_end * OCE_msk3) 474 450 475 echo ( 'OCE_sum_ssh_beg = {:12.6e} m^3 - OCE_sum_ssh_end = {:12.6e} m^3'.format (OCE_sum_ssh_beg, OCE_sum_ssh_end))451 echo ( f'OCE_sum_ssh_beg = {OCE_sum_ssh_beg:12.6e} m^3 - OCE_sum_ssh_end = {OCE_sum_ssh_end:12.6e} m^3' ) 476 452 dOCE_ssh_vol = ( OCE_sum_ssh_end - OCE_sum_ssh_beg ) 477 453 dOCE_ssh_mas = dOCE_ssh_vol * OCE_rho_liq 478 454 479 455 if NEMO == 3.6 : 480 echo ( 'OCE_sum_e3tn_beg = {:12.6e} m^3 - OCE_sum_e3tn_end = {:12.6e} m^3'.format (OCE_sum_e3tn_beg, OCE_sum_e3tn_end))456 echo ( f'OCE_sum_e3tn_beg = {OCE_sum_e3tn_beg:12.6e} m^3 - OCE_sum_e3tn_end = {OCE_sum_e3tn_end:12.6e} m^3' ) 481 457 dOCE_e3tn_vol = ( OCE_sum_e3tn_end - OCE_sum_e3tn_beg ) 482 458 dOCE_e3tn_mas = dOCE_e3tn_vol * OCE_rho_liq … … 484 460 dOCE_vol_wat = dOCE_ssh_vol ; dOCE_mas_wat = dOCE_ssh_mas 485 461 486 echo ( 'dOCE vol = {:12.3e} m^3'.format ( dOCE_vol_wat))487 echo ( 'dOCE ssh = {:12.3e} m '.format ( dOCE_vol_wat/OCE_aire_tot))462 echo ( f'dOCE vol = {dOCE_vol_wat :12.3e} m^3' ) 463 echo ( f'dOCE ssh = {dOCE_vol_wat/OCE_aire_tot:12.3e} m ' ) 488 464 prtFlux ( 'dOCE mass ', dOCE_mas_wat, 'e' ) 489 465 490 466 if NEMO == 3.6 : 491 echo ( 'dOCE e3tn vol = {:12.3e} m^3'.format ( dOCE_e3tn_vol))467 echo ( f'dOCE e3tn vol = {dOCE_e3tn_vol:12.3e} m^3' ) 492 468 prtFlux ( 'dOCE e3tn mass', dOCE_e3tn_mas, 'e' ) 493 469 … … 547 523 dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat 548 524 549 echo ( 'ICE_vol_ice_beg = {:12.6e} m^3 | ICE_vol_ice_end = {:12.6e} m^3'.format (ICE_vol_ice_beg, ICE_vol_ice_end))550 echo ( 'ICE_vol_sno_beg = {:12.6e} m^3 | ICE_vol_sno_end = {:12.6e} m^3'.format (ICE_vol_sno_beg, ICE_vol_sno_end))551 echo ( 'ICE_vol_pnd_beg = {:12.6e} m^3 | ICE_vol_pnd_end = {:12.6e} m^3'.format (ICE_vol_pnd_beg, ICE_vol_pnd_end))552 echo ( 'ICE_vol_fzl_beg = {:12.6e} m^3 | ICE_vol_fzl_end = {:12.6e} m^3'.format (ICE_vol_fzl_beg, ICE_vol_fzl_end))553 554 echo ( 'dICE_vol_ice = {:12.3e} m^3'.format (dICE_vol_ice))555 echo ( 'dICE_vol_sno = {:12.3e} m^3'.format (dICE_vol_sno))556 echo ( 'dICE_vol_pnd = {:12.3e} m^3'.format (dICE_vol_pnd))557 echo ( 'dICE_mas_ice = {:12.3e} m^3'.format (dICE_mas_ice))558 echo ( 'dICE_mas_sno = {:12.3e} m^3'.format (dICE_mas_sno))559 echo ( 'dICE_mas_pnd = {:12.3e} m^3'.format (dICE_mas_pnd))560 echo ( 'dICE_mas_fzl = {:12.3e} m^3'.format (dICE_mas_fzl))525 echo ( f'ICE_vol_ice_beg = {ICE_vol_ice_beg:12.6e} m^3 | ICE_vol_ice_end = {ICE_vol_ice_end:12.6e} m^3' ) 526 echo ( f'ICE_vol_sno_beg = {ICE_vol_sno_beg:12.6e} m^3 | ICE_vol_sno_end = {ICE_vol_sno_end:12.6e} m^3' ) 527 echo ( f'ICE_vol_pnd_beg = {ICE_vol_pnd_beg:12.6e} m^3 | ICE_vol_pnd_end = {ICE_vol_pnd_end:12.6e} m^3' ) 528 echo ( f'ICE_vol_fzl_beg = {ICE_vol_fzl_beg:12.6e} m^3 | ICE_vol_fzl_end = {ICE_vol_fzl_end:12.6e} m^3' ) 529 530 echo ( f'dICE_vol_ice = {dICE_vol_ice:12.3e} m^3' ) 531 echo ( f'dICE_vol_sno = {dICE_vol_sno:12.3e} m^3' ) 532 echo ( f'dICE_vol_pnd = {dICE_vol_pnd:12.3e} m^3' ) 533 echo ( f'dICE_mas_ice = {dICE_mas_ice:12.3e} m^3' ) 534 echo ( f'dICE_mas_sno = {dICE_mas_sno:12.3e} m^3' ) 535 echo ( f'dICE_mas_pnd = {dICE_mas_pnd:12.3e} m^3' ) 536 echo ( f'dICE_mas_fzl = {dICE_mas_fzl:12.3e} m^3' ) 561 537 562 538 echo ( '\n------------------------------------------------------------') 563 echo ( ' Variation du contenu en eau ocean + glace' )539 echo ( 'Change in water content (ocean + ice) ' ) 564 540 prtFlux ( 'dMass (ocean)', dSEA_mas_wat, 'e', True ) 565 541 … … 584 560 585 561 echo ( '\n------------------------------------------------------------------------------------' ) 586 echo ( '-- Checks in NEMO - from budget_modipsl.sh (Cl ement Rousset)' )562 echo ( '-- Checks in NEMO - from budget_modipsl.sh (Clément Rousset)' ) 587 563 588 564 # Read variable and computes integral over space and time … … 658 634 prtFlux (' WFXSNW_PRE ', ICE_mas_wfxsnw_pre, 'e', True) 659 635 prtFlux (' WFXSUB_ERR ', ICE_mas_wfxsub_err, 'e', True) 660 prtFlux (' EVAP_OCE ', OCE_mas_evap_oce , 'e' )636 prtFlux (' EVAP_OCE ', OCE_mas_evap_oce , 'e' ) 661 637 prtFlux (' EVAP_ICE ', ICE_mas_evap_ice , 'e', True) 662 638 prtFlux (' SNOW_OCE ', OCE_mas_snow_oce , 'e', True) 663 639 prtFlux (' SNOW_ICE ', OCE_mas_snow_ice , 'e', True) 664 prtFlux (' RAIN ', OCE_mas_rain , 'e' )640 prtFlux (' RAIN ', OCE_mas_rain , 'e' ) 665 641 prtFlux (' FWB ', OCE_mas_fwb , 'e', True) 666 642 prtFlux (' SSR ', OCE_mas_ssr , 'e', True) … … 705 681 706 682 echo ( '\n------------------------------------------------------------------------------------' ) 707 echo ( '-- C alculs dans la note PDF de Clement')683 echo ( '-- Computations in the PDF note of Clément Rousset') 708 684 echo ( 'Freshwater budget of the ice-ocean system = emp_oce + emp_ice - runoffs - iceberg - iceshelf = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceberg - OCE_mas_iceshelf )) 709 685 echo ( 'Freshwater budget of the ice-ocean system = emp_oce + emp_ice - friver - iceberg - iceshelf = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_friver - OCE_mas_iceberg - OCE_mas_iceshelf )) -
TOOLS/WATER_BUDGET/lmdz.py
r6265 r6508 238 238 if math == xr : 239 239 p2D = xr.DataArray (p2D) 240 p2D.attrs = p1D.attrs 240 for attr in p1D.attrs : 241 p2D.attrs[attr] = p1D.attrs[attr] 241 242 p2D = p2D.rename ( {p2D.dims[0]:p1D.dims[0], p2D.dims[-1]:'x', p2D.dims[-2]:'y'} ) 242 243 243 244 return p2D 244 245 245 def geo2point ( p2D, cumulPoles=False, dim1D=' points_physiques') :246 def geo2point ( p2D, cumulPoles=False, dim1D='cell', debug=False ) : 246 247 ''' 247 248 From 2D (lat, lon) to 1D (points_phyiques) … … 259 260 260 261 jpn = jpi*(jpj-2) + 2 262 263 if debug : 264 print ( f'{jpi=} {jpj=} {jpn=} {jpt=}' ) 261 265 262 266 if jpt == -1 : … … 277 281 if math == xr : 278 282 p1D = xr.DataArray (p1D) 279 p1D.attrs = p2D.attrs 283 for attr in p2D.attrs : 284 p1D.attrs[attr] = p2D.attrs[attr] 280 285 p1D = p1D.rename ( {p1D.dims[0]:p2D.dims[0], p1D.dims[-1]:dim1D} ) 281 286 … … 286 291 return p1D 287 292 288 def geo3point ( p3D, cumulPoles=False, dim1D=' points_physiques' ) :293 def geo3point ( p3D, cumulPoles=False, dim1D='cell' ) : 289 294 ''' 290 295 From 3D (lev, lat, lon) to 2D (lev, points_phyiques) … … 314 319 if math == xr : 315 320 p2D = xr.DataArray (p2D) 316 p2D.attrs = p3D.attrs 321 for attr in p2D.attrs : 322 p2D.attrs[attr] = p3D.attrs[attr] 317 323 p2D = p2D.rename ( {p2D.dims[-1]:dim1D, p2D.dims[-2]:p3D.dims[-3]} ) 318 324
Note: See TracChangeset
for help on using the changeset viewer.