source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6518

Last change on this file since 6518 was 6518, checked in by omamce, 13 months ago

O.M. - WATER_BUDGET
Evolutions pour Sechiba
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 67.9 KB
Line 
1#!/usr/bin/env python3
2###
3### Script to check water conservation in the IPSL coupled model
4###
5##  Warning, to install, configure, run, use any of included software or
6##  to read the associated documentation you'll need at least one (1)
7##  brain in a reasonably working order. Lack of this implement will
8##  void any warranties (either express or implied).  Authors assumes
9##  no responsability for errors, omissions, data loss, or any other
10##  consequences caused directly or indirectly by the usage of his
11##  software by incorrectly or partially configured personal
12##
13##
14# SVN Information
15SVN = {
16    'Author'  : "$Author$",
17    'Date'    : "$Date$",
18    'Revision': "$Revision$",
19    'Id'      : "$Id: ATM_waterbudget.py 6508 2023-06-13 10:58:38Z omamce $",
20    'HeadURL' : "$HeadUrl: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/WATER_BUDGET/ATM_waterbudget.py $"
21    }
22
23###
24## Import system modules
25import sys, os, shutil, subprocess, platform
26import numpy as np
27import configparser, re
28
29# Check python version
30if 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
35import WaterUtils as wu
36
37## Creates parser for reading .ini input file
38config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
39config.optionxform = str # To keep capitals
40
41## Experiment parameters
42ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
43OCE_icb=False ; Coupled=False ; Routing=None
44TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
45
46##
47ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
48TmpDir=None ; RunDir=None ; FileOut=None
49dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
50FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
51file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
52tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
53file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
54file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
55file_ICE_beg=None ; file_OCE_beg=None
56file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
57tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
58tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
59tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
60tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
61ContinueOnError=False ; ErrorCount=0 ; SortIco = False
62
63# Read command line arguments
64print ( "Name of Python script:", sys.argv[0] )
65IniFile = sys.argv[1]
66# Text existence of IniFile
67if not os.path.exists (IniFile ) :
68    raise FileExistsError ( f"File not found : {IniFile = }" )
69
70if 'full' in IniFile : FullIniFile = IniFile
71else                 : FullIniFile = 'full_' + IniFile
72   
73print ("Input file : ", IniFile )
74config.read (IniFile)
75FullIniFile = 'full_' + IniFile
76
77## Reading config file
78for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
79   if Section in config.keys () : 
80        print ( f'\nReading [{Section}]' )
81        for VarName in config[Section].keys() :
82            locals()[VarName] = config[Section][VarName]
83            exec  ( f'{VarName} = wu.setBool ({VarName})' )
84            exec  ( f'{VarName} = wu.setNum  ({VarName})' )
85            exec  ( f'{VarName} = wu.setNone ({VarName})' )
86            exec  ( f'wu.{VarName} = {VarName}' )
87            print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
88            #exec ( f'del {VarName}' )
89
90print ( f'\nConfig file readed : {IniFile} ' )
91           
92##-- Some physical constants
93if wu.unDefined ( 'Ra' )           : Ra          = 6366197.7236758135 #-- Earth Radius (m)
94if wu.unDefined ( 'Grav' )         : Grav        = 9.81               #-- Gravity (m^2/s)
95if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = 917.0              #-- Ice volumic mass (kg/m3) in LIM3 or SI3
96if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3 or SI3
97if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = 1026.0             #-- Ocean water volumic mass (kg/m3) in NEMO
98if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = 1000.0             #-- Water volumic mass in atmosphere (kg/m^3)
99if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = 1000.0             #-- Water volumic mass in surface reservoir (kg/m^3)
100if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = 1000.0             #-- Water volumic mass of rivers (kg/m^3)
101if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = 1000.              #-- Water volumic mass in ice ponds (kg/m^3)
102if wu.unDefined ( 'YearLength' )   : YearLength  = 365.25 * 86400.    #-- Year length (s) - Use only to compute drif in approximate unit
103           
104##-- Set libIGCM and machine dependant values
105if not 'Files' in config.keys () : config['Files'] = {}
106
107config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno,
108                      'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
109
110config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco }
111
112## --------------------------
113ICO  = ( 'ICO' in wu.ATM )
114LMDZ = ( 'LMD' in wu.ATM )
115   
116with open ('SetLibIGCM.py') as f: exec ( f.read() )
117config['Files']['TmpDir'] = TmpDir
118
119if libIGCM :
120    config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
121
122##-- Import specific module
123import nemo, lmdz
124##-- Now import needed scientific modules
125import xarray as xr
126
127##- Output file with water budget diagnostics
128if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
129config['Files']['FileOut'] = FileOut
130
131f_out = open ( FileOut, mode = 'w' )
132
133##- Useful functions
134def kg2Sv    (val, rho=ATM_rho) :
135    '''From kg to Sverdrup'''
136    return val/dtime_sec*1.0e-6/rho
137
138def kg2myear (val, rho=ATM_rho) :
139    '''From kg to m/year'''
140    return val/ATM_aire_sea_tot/rho/NbYear
141
142def var2prt (var, small=False, rho=ATM_rho) :
143    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
144    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
145
146def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
147    if small :
148        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
149        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
150    else :
151        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
152        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
153    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) )
154    return None
155
156def echo (string, end='\n') :
157    '''Function to print to stdout *and* output file'''
158    print ( str(string), end=end  )
159    sys.stdout.flush ()
160    f_out.write ( str(string) + end )
161    f_out.flush ()
162    return None
163   
164##- Set libIGCM directories
165if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
166if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
167if wu.unDefined ('L_EXP'      ) : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName)
168if wu.unDefined ('R_SAVE'     ) : R_SAVE      = os.path.join ( R_OUT, L_EXP )
169if wu.unDefined ('R_BUFR'     ) : R_BUFR      = os.path.join ( R_BUF, L_EXP )
170if wu.unDefined ('POST_DIR'   ) : POST_DIR    = os.path.join ( R_BUFR, 'Out' )
171if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' )
172if wu.unDefined ('R_BUF_KSH'  ) : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
173if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
174
175config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
176                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR }
177
178##- Set directory to extract files
179if wu.unDefined ( 'RunDir' ) : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
180config['Files']['RunDir'] = RunDir
181
182if not os.path.isdir ( RunDir ) : os.makedirs ( RunDir )
183
184##- Set directories to rebuild ocean and ice restart files
185if wu.unDefined ( 'RunDirOCE' ) : RunDirOCE = os.path.join ( RunDir, 'OCE' )
186if wu.unDefined ( 'RunDirICE' ) : RunDirICE = os.path.join ( RunDir, 'ICE' )
187if not os.path.exists ( RunDirOCE ) : os.mkdir ( RunDirOCE )
188if not os.path.exists ( RunDirICE ) : os.mkdir ( RunDirICE )
189
190echo (' ')
191echo ( f'JobName   : {JobName}'   )
192echo ( f'Comment   : {Comment}'   )
193echo ( f'TmpDir    : {TmpDir}'    )
194echo ( f'RunDir    : {RunDir}'    )
195echo ( f'RunDirOCE : {RunDirOCE}' )
196echo ( f'RunDirICE : {RunDirICE}' )
197
198echo ( f'\nDealing with {L_EXP}'  )
199
200##-- Creates model output directory names
201if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' )
202if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' )
203if wu.unDefined ('dir_ATM_his' ) :
204    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
205    config['Files']['dir_ATM_his'] = dir_ATM_his
206if wu.unDefined ( 'dir_SRF_his' ) : 
207    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
208    config['Files']['dir_SRF_his'] = dir_SRF_his
209   
210echo ( f'The analysis relies on files from the following model output directories : ' )
211echo ( f'{dir_ATM_his = }' )
212echo ( f'{dir_SRF_his = }' )
213
214##-- Creates files names
215if wu.unDefined ( 'Period' ) :
216    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
217    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
218    config['Files']['Period'] = Period
219if wu.unDefined ( 'FileCommon' ) :
220    FileCommon = f'{JobName}_{Period}'
221    config['Files']['FileCommon'] = FileCommon
222
223if wu.unDefined ( 'Title' ) :
224    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
225    config['Files']['Title'] = Title
226     
227echo ('\nOpen history files' )
228if wu.unDefined ( 'file_ATM_his' ) :
229    if ATM_HIS == 'latlon' :
230        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' )
231    if ATM_HIS == 'ico' :
232        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' )
233    config['Files']['file_ATM_his'] = file_ATM_his
234if wu.unDefined ( 'file_SRF_his' ) :
235    if ATM_HIS == 'latlon' :
236        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
237    if ATM_HIS == 'ico' :
238        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
239    config['Files']['file_SRF_his'] = file_SRF_his
240
241if Routing == 'SIMPLE' :
242    if file_RUN_his == None :
243        if ATM_HIS == 'latlon' :
244            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
245        if ATM_HIS == 'ico' :
246            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
247        config['Files']['file_RUN_his'] = file_RUN_his
248
249d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
250d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
251if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
252if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
253
254echo ( f'{file_ATM_his = }' )
255echo ( f'{file_SRF_his = }' )
256if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
257
258##-- Compute run length
259dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
260echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
261dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
262
263##-- Compute length of each period
264dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
265echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
266dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
267dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
268dtime_per_sec.attrs['unit'] = 's'
269
270##-- Number of years (approximative)
271NbYear = dtime_sec / YearLength
272
273##-- Extract restart files from tar
274YearRes = YearBegin - 1              # Year of the restart of beginning of simulation
275YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
276
277config['Files']['YearPre'] = f'{YearBegin}'
278
279if wu.unDefined ( 'TarRestartPeriod_beg' ) : 
280    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
281    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231'
282    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
283
284if wu.unDefined ( 'TarRestartPeriod_end ' ) : 
285    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
286    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
287    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231'
288    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
289
290if wu.unDefined ( 'tar_restart_beg' ) :
291    tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' )
292    config['Files']['tar_restart_beg'] = tar_restart_beg
293if wu.unDefined ( 'tar_restart_end' ) :
294    tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' )
295    config['Files']['tar_restart_end'] = tar_restart_end
296   
297echo ( f'{tar_restart_beg = }' )
298echo ( f'{tar_restart_end = }' )
299
300##-- Names of tar files with restarts
301if wu.unDefined ( 'SRF_HIS' ) : SRF_HIS = ATM_HIS
302
303if wu.unDefined ( 'tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg
304if wu.unDefined ( 'tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg
305if wu.unDefined ( 'tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg
306if wu.unDefined ( 'tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg
307if wu.unDefined ( 'tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg
308if wu.unDefined ( 'tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg
309
310if wu.unDefined ( 'tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end
311if wu.unDefined ( 'tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end
312if wu.unDefined ( 'tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end
313if wu.unDefined ( 'tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end
314if wu.unDefined ( 'tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end
315if wu.unDefined ( 'tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end
316
317if wu.unDefined ( 'file_ATM_beg' ) : 
318    file_ATM_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc'
319    config['Files']['file_ATM_beg'] = file_ATM_beg
320if wu.unDefined ( 'file_ATM_end' ) :
321    file_ATM_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc'
322    config['Files']['file_ATM_end'] = file_ATM_end
323
324liste_beg = [file_ATM_beg, ]
325liste_end = [file_ATM_end, ]
326
327if wu.unDefined ( 'file_DYN_beg' ) :
328    if LMDZ : file_DYN_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restart.nc'
329    if ICO  : file_DYN_beg = f'{RunDir}/ICO_{JobName}_{YearRes}1231_restart.nc'
330    liste_beg.append (file_DYN_beg)
331    config['Files']['file_DYN_beg'] = file_DYN_beg
332   
333if wu.unDefined ( 'file_DYN_end' ) : 
334    if LMDZ : file_DYN_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restart.nc'
335    if ICO  : file_DYN_end = f'{RunDir}/ICO_{JobName}_{YearEnd}1231_restart.nc'
336    liste_end.append ( file_DYN_end )
337    config['Files']['file_DYN_end'] = file_DYN_end
338
339if wu.unDefined ( 'file_SRF_beg' ) :
340    file_SRF_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'
341    config['Files']['file_SRF_beg'] = file_SRF_beg
342if wu.unDefined ( 'file_SRF_end' ) :
343    file_SRF_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc'
344    config['Files']['file_SRF_end'] = file_SRF_end
345   
346liste_beg.append ( file_SRF_beg )
347liste_end.append ( file_SRF_end )
348
349echo ( f'{file_ATM_beg = }')
350echo ( f'{file_ATM_end = }')
351echo ( f'{file_DYN_beg = }')
352echo ( f'{file_DYN_end = }')
353echo ( f'{file_SRF_beg = }')
354echo ( f'{file_SRF_end = }')
355
356if ICO :
357    if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
358    config['Files'] ['file_DYN_aire'] = file_DYN_aire
359
360if Routing == 'SIMPLE' :
361    if wu.unDefined ( 'file_RUN_beg' ) : 
362        file_RUN_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc'
363        config['Files']['file_RUN_beg'] = file_RUN_beg
364    if wu.unDefined ( 'file_RUN_end' ) : 
365        file_RUN_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc'
366        config['Files']['file_RUN_end'] = file_RUN_end
367       
368    liste_beg.append ( file_RUN_beg )
369    liste_end.append ( file_RUN_end )
370    echo ( f'{file_RUN_beg = }' )
371    echo ( f'{file_RUN_end = }' )
372
373echo ('\nExtract restart files from tar : ATM, ICO and SRF')
374
375def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, RunDirComp=RunDir, ErrorCount=ErrorCount ) :
376    echo ( f'----------')
377    echo ( f'file to extract : {file_name = }' )
378    if os.path.exists ( os.path.join (RunDir, file_name) ) :
379        echo ( f'file found : {file_name = }' )
380    else :
381        echo ( f'file not found : {file_name = }' )
382        base_resFile = os.path.basename (file_name)
383        if os.path.exists ( tar_restart ) :
384            command =  f'cd {RunDir} ; tar xf {tar_restart} {base_resFile}'
385            echo ( f'{command = }' )
386            try : 
387                os.system ( command )
388            except :
389                if ContinueOnError :
390                    ErrorCount += 1
391                    echo ( f'****** Command failed : {command}' )
392                    echo ( f'****** Trying to continue' )
393                    echo ( ' ')
394                else :
395                    raise Exception ( f'**** command failed : {command} - Stopping' )
396            else :
397                echo ( f'tar done : {base_resFile}' )
398        else :
399            echo ( f'****** Tar restart file {tar_restart = } not found ' )
400            if ContinueOnError :
401                ErrorCount += 1
402            else : 
403                raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' )
404    return ErrorCount
405 
406ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, RunDirComp=RunDir )
407ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, RunDirComp=RunDir )
408ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, RunDirComp=RunDir )
409
410ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, RunDirComp=RunDir )
411ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, RunDirComp=RunDir )
412ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, RunDirComp=RunDir )
413
414if Routing == 'SIMPLE' :
415    ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, RunDirComp=RunDir )
416    ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, RunDirComp=RunDir )
417
418##-- Exit in case of error in the opening file phase
419if ErrorCount > 0 :
420    echo ( '  ' )
421    raise Exception ( f'**** Some files missing - Stopping - {ErrorCount = }' )
422
423##
424echo ('\nOpening ATM SRF and ICO restart files')
425d_ATM_beg = xr.open_dataset ( os.path.join (RunDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
426d_ATM_end = xr.open_dataset ( os.path.join (RunDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
427d_SRF_beg = xr.open_dataset ( os.path.join (RunDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
428d_SRF_end = xr.open_dataset ( os.path.join (RunDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
429d_DYN_beg = xr.open_dataset ( os.path.join (RunDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
430d_DYN_end = xr.open_dataset ( os.path.join (RunDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
431
432for var in d_SRF_beg.variables :
433    d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.)
434    d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.)
435
436if Routing == 'SIMPLE' :
437    d_RUN_beg = xr.open_dataset ( os.path.join (RunDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
438    d_RUN_end = xr.open_dataset ( os.path.join (RunDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
439
440def to_cell ( dd, newname='cell' ) :
441    '''Set space dimension to 'cell' '''
442    for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] :
443        try    : dd = dd.rename ( {oldname : newname} )
444        except : pass
445    return dd
446
447d_ATM_beg = to_cell ( d_ATM_beg )
448d_ATM_end = to_cell ( d_ATM_end )
449d_SRF_beg = to_cell ( d_SRF_beg )
450d_SRF_end = to_cell ( d_SRF_end )
451d_DYN_beg = to_cell ( d_DYN_beg )
452d_DYN_end = to_cell ( d_DYN_end )
453
454if Routing == 'SIMPLE' :
455    d_RUN_beg = to_cell ( d_RUN_beg )
456    d_RUN_end = to_cell ( d_RUN_end )
457
458d_ATM_his = to_cell ( d_ATM_his )
459d_SRF_his = to_cell ( d_SRF_his )
460   
461echo ( f'{file_ATM_beg = }' )
462echo ( f'{file_ATM_end = }' )
463echo ( f'{file_DYN_beg = }' )
464echo ( f'{file_DYN_end = }' )
465echo ( f'{file_SRF_beg = }' )
466echo ( f'{file_SRF_end = }' )
467if Routing == 'SIMPLE' :
468    echo ( f'{file_RUN_beg = }' )
469    echo ( f'{file_RUN_end = }' )
470
471## Write the full configuration
472#config_out = open (FullIniFile, 'w')
473#config.write ( config_out )
474#config_out.close ()
475
476if ICO :
477    if SortIco : 
478        # Creation d'une clef de tri pour les restarts (pour chaque restart !)
479        DYN_beg_keysort = np.lexsort ( (d_DYN_beg['lat_mesh'], d_DYN_beg['lon_mesh'] ) )
480        ATM_beg_keysort = np.lexsort ( (d_ATM_beg['latitude'], d_ATM_beg['longitude']) )
481        SRF_beg_keysort = np.lexsort ( (d_SRF_beg['nav_lat' ], d_SRF_beg['nav_lon']  ) )
482        DYN_end_keysort = np.lexsort ( (d_DYN_end['lat_mesh'], d_DYN_end['lon_mesh'] ) )
483        ATM_end_keysort = np.lexsort ( (d_ATM_end['latitude'], d_ATM_end['longitude']) )
484        SRF_end_keysort = np.lexsort ( (d_SRF_end['nav_lat' ], d_SRF_end['nav_lon']  ) )
485       
486        if ATM_HIS == 'ico' :
487            ATM_his_keysort = np.lexsort ( (d_ATM_his['lat'], d_ATM_his['lon'] ) )
488        if SRF_HIS == 'ico' :
489            SRF_his_keysort = np.lexsort ( (d_SRF_his['lat'], d_SRF_his['lon'] ) )
490            RUN_his_keysort = SRF_his_keysort
491    else : 
492        DYN_beg_keysort = np.arange ( len ( d_DYN_beg['lat_mesh'] ) )
493        ATM_beg_keysort = DYN_beg_keysort
494        SRF_beg_keysort = DYN_beg_keysort
495
496        DYN_end_keysort = DYN_beg_keysort
497        ATM_end_keysort = DYN_beg_keysort
498        SRF_end_keysort = DYN_beg_keysort
499   
500        ATM_his_keysort = DYN_beg_keysort
501        SRF_his_keysort = DYN_beg_keysort
502        RUN_his_keysort = SRF_his_keysort
503       
504# ATM grid with cell surfaces
505if LMDZ :
506    ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' )
507    ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' )
508    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire']     [0], cumulPoles=True, dim1D='cell_latlon' )
509    ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' )
510    ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' )
511    ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' )
512    ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' )
513    SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' )
514    SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' )
515    SRF_aire      = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell_latlonl', cumulPoles=True )
516    SRF_areas     = lmdz.geo2point ( d_SRF_his ['Areas'],  dim1D='cell_latlonl', cumulPoles=True )
517    SRF_contfrac  = lmdz.geo2point ( d_SRF_his ['Contfrac'], dim1D='cell_latlonl' )
518
519    #SRF_res_aire = SRF_aire
520   
521if ICO :
522    if ATM_HIS == 'latlon' :
523        jpja, jpia = d_ATM_his['aire'][0].shape
524        if wu.unDefined ( 'file_ATM_aire' ) : 
525            file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' )
526            config['Files']['file_ATM_aire'] = file_ATM_aire
527        echo ( f'Aire sur grille reguliere : {file_ATM_aire = }' )
528        d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze()
529        try :
530            ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' )
531            ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' )
532        except :
533            ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat_dom_out']+0*d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' )
534            ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat_dom_out']+  d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' )
535        ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True, dim1D='cell_latlon' )
536        ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' )
537        ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' )
538        ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' )
539        ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' )
540
541    if SRF_HIS == 'latlon' :
542        try :
543            SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' )
544            SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' )
545        except :
546            try : 
547                SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_dom_out']+0*d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' )
548                SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_dom_out']+  d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' )
549            except : 
550                SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_domain_landpoints_out']+0*d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' )
551                SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_domain_landpoints_out']+  d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' )
552        SRF_areas    = d_SRF_his ['Areas'].values
553        SRF_areafrac = d_SRF_his ['AreaFrac'].values
554        SRF_areas    = xr.DataArray ( SRF_areas   , coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
555        SRF_areafrac = xr.DataArray ( SRF_areafrac, coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
556        SRF_contfrac = SRF_areafrac / SRF_areas
557        #SRF_aire      = lmdz.geo2point (SRF_areas * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
558        #SRF_aire      = lmdz.geo2point (SRF_areas, dim1D='cell_latlon', cumulPoles=True )
559        #SRF_aire      = lmdz.geo2point (d_ATM_aire ['aire'].squeeze().values * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
560        SRF_areas     = lmdz.geo2point ( SRF_areas   , dim1D='cell_latlon', cumulPoles=True )
561        SRF_areafrac  = lmdz.geo2point ( SRF_areafrac, dim1D='cell_latlon', cumulPoles=True )
562        SRF_contfrac  = lmdz.geo2point ( SRF_contfrac, dim1D='cell_latlon' )
563        SRF_aire      = SRF_areafrac
564
565    if ATM_HIS == 'ico' :
566        echo ( f'Grille icosaedre' )
567        ATM_aire =  d_ATM_his ['aire']     [0][ATM_his_keysort].squeeze()
568        ATM_lat  =  d_ATM_his ['lat']         [ATM_his_keysort]
569        ATM_lon  =  d_ATM_his ['lon']         [ATM_his_keysort]
570        ATM_fter =  d_ATM_his ['fract_ter'][0][ATM_his_keysort]
571        ATM_foce =  d_ATM_his ['fract_oce'][0][ATM_his_keysort]
572        ATM_fsic =  d_ATM_his ['fract_sic'][0][ATM_his_keysort]
573        ATM_flic =  d_ATM_his ['fract_lic'][0][ATM_his_keysort]
574       
575    if SRF_HIS == 'ico' : 
576        SRF_lat       =  d_SRF_his ['lat'][SRF_his_keysort]
577        SRF_lon       =  d_SRF_his ['lon'][SRF_his_keysort]
578        SRF_aire      =  d_SRF_his ['Areas'][SRF_his_keysort] * d_SRF_his ['Contfrac'][SRF_his_keysort]
579        SRF_areas     =  d_SRF_his ['Areas'][SRF_his_keysort] 
580        SRF_contfrac  =  d_SRF_his ['Contfrac'][SRF_his_keysort]
581
582ATM_fsea      = ATM_foce + ATM_fsic
583ATM_flnd      = ATM_fter + ATM_flic
584ATM_aire_fter = ATM_aire * ATM_fter
585ATM_aire_flic = ATM_aire * ATM_flic
586ATM_aire_fsic = ATM_aire * ATM_fsic
587ATM_aire_foce = ATM_aire * ATM_foce
588ATM_aire_flnd = ATM_aire * ATM_flnd
589ATM_aire_fsea = ATM_aire * ATM_fsea
590
591SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0.)
592
593## Write the full configuration
594config_out = open (FullIniFile, 'w')
595config.write ( config_out )
596config_out.close ()
597           
598if ICO :
599    # Area on icosahedron grid
600    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
601
602    if SortIco : 
603        # Creation d'une clef de tri pour le fichier aire
604        DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) )
605    else : 
606        DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) )
607           
608    DYN_lat = d_DYN_aire['lat'][DYN_aire_keysort]
609    DYN_lon = d_DYN_aire['lon'][DYN_aire_keysort]
610
611    DYN_aire = d_DYN_aire['aire'][DYN_aire_keysort]
612    DYN_fsea = d_DYN_aire['fract_oce'][DYN_aire_keysort] + d_DYN_aire['fract_sic'][DYN_aire_keysort]
613    DYN_flnd = 1.0 - DYN_fsea
614    DYN_fter = d_ATM_beg['FTER'][ATM_beg_keysort]
615    DYN_flic = d_ATM_beg['FLIC'][ATM_beg_keysort]
616    DYN_aire_fter = DYN_aire * DYN_fter
617   
618if LMDZ :
619    # Area on lon/lat grid
620    DYN_aire = ATM_aire
621    DYN_fsea = ATM_fsea
622    DYN_flnd = ATM_flnd
623    DYN_fter = d_ATM_beg['FTER']
624    DYN_flic = d_ATM_beg['FLIC']
625    DYN_aire_fter = DYN_aire * DYN_fter
626
627# Functions computing integrals and sum
628def DYN_stock_int (stock) :
629    '''Integrate (* surface) stock on atmosphere grid'''
630    DYN_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
631    return DYN_stock_int
632
633def ATM_flux_int (flux) :
634    '''Integrate (* time * surface) flux on atmosphere grid'''
635    ATM_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() )
636    return ATM_flux_int
637
638def SRF_stock_int (stock) :
639    '''Integrate (* surface) stock on land grid'''
640    #SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) # Marche avec ICO et lon/lat na
641    SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
642    return SRF_stock_int
643
644def SRF_flux_int (flux) :
645    '''Integrate (* time * surface) flux on land grid'''
646    SRF_flux_int  = wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() )
647    return SRF_flux_int
648
649def ONE_stock_int (stock) :
650    '''Sum stock '''
651    ONE_stock_int  = wu.Psum ( stock.to_masked_array().ravel() )
652    return ONE_stock_int
653
654def ONE_flux_int (flux) :
655    '''Integrate (* time) flux on area=1 grid '''
656    ONE_flux_int  = wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() )
657    return ONE_flux_int
658
659ATM_aire_sea     = ATM_aire * ATM_fsea
660ATM_aire_tot     = ONE_stock_int (ATM_aire)
661ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea)
662
663DYN_aire_tot     =  ONE_stock_int ( DYN_aire )
664SRF_aire_tot     =  ONE_stock_int ( SRF_aire )
665
666echo ('')
667echo ( f'ATM DYN : Area of atmosphere : {DYN_aire_tot:18.9e}' )
668echo ( f'ATM HIS : Area of atmosphere : {ATM_aire_tot:18.9e}' )
669echo ( f'ATM SRF : Area of atmosphere : {SRF_aire_tot:18.9e}' )
670echo ('')
671echo ( 'ATM DYN : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) )
672echo ( 'ATM HIS : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )
673echo ( 'ATM SRF : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) )
674echo ('')
675echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' )
676
677
678if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) :
679    raise Exception ('Error of atmosphere surface interpolated on lon/lat grid')
680
681echo ( '\n====================================================================================' )
682echo ( f'-- ATM changes in stores  -- {Title} ' )
683
684#-- Change in precipitable water from the atmosphere daily and monthly files
685#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
686
687echo ( 'ATM vertical grid' )
688ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
689ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
690
691echo ( 'Surface pressure' )
692if ICO :
693    DYN_psol_beg = d_DYN_beg['ps'][DYN_beg_keysort]
694    DYN_psol_end = d_DYN_end['ps'][DYN_end_keysort]
695if LMDZ : 
696    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' )
697    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' )
698   
699echo ( '3D Pressure at the interface layers (not scalar points)' )
700DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg
701DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end
702
703klevp1 = ATM_Bhyb.shape[-1]
704cell   = DYN_psol_beg.shape[-1]
705klev   = klevp1 - 1
706
707echo ( 'Layer thickness (pressure)' )
708DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
709DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
710   
711for k in np.arange (klevp1-1) :
712    DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:]
713    DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:]
714
715echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' )
716if LMDZ :
717    try : 
718        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']  + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' )
719        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']  + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' )
720    except :
721        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' )
722        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' )
723if ICO :
724    try :
725        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'])[..., DYN_beg_keysort]
726        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'])[..., DYN_beg_keysort]
727       
728    except : 
729        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]
730        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]
731       
732    DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} )
733    DYN_wat_beg = DYN_wat_end.rename ( {'lev':'sigs'} )
734   
735echo ( 'Compute water content : vertical and horizontal integral' )
736DYN_mas_wat_beg = DYN_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav
737DYN_mas_wat_end = DYN_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav
738
739echo ( 'Variation of water content' )
740dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg
741
742# \([a-z,A-Z,_]*\)/dtime_sec\*1e-9  kg2Sv(\1)
743# \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear  kg2myear(\1)
744
745echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' )
746echo ( '------------------------------------------------------------------------------------' )
747echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
748prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True )
749
750ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \
751              d_ATM_beg['SNOW03']*d_ATM_beg['FOCE'] + d_ATM_beg['SNOW04']*d_ATM_beg['FSIC']
752ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \
753              d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC']
754
755ATM_qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \
756              d_ATM_beg['QS03']  *d_ATM_beg['FOCE'] + d_ATM_beg['QS04']  *d_ATM_beg['FSIC']
757ATM_qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \
758              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC']
759
760ATM_qsol_beg = d_ATM_beg['QSOL']
761ATM_qsol_end = d_ATM_end['QSOL']
762
763LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']
764LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC']
765
766LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0']
767LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0']
768
769ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
770ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER']
771ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
772ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC']
773ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
774ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE']
775ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
776ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
777
778if ICO :
779     ATM_sno_beg       = ATM_sno_beg    [ATM_beg_keysort]
780     ATM_sno_end       = ATM_sno_end    [ATM_end_keysort]
781     ATM_qs_beg        = ATM_qs_beg     [ATM_beg_keysort]
782     ATM_qs_end        = ATM_qs_end     [ATM_end_keysort]
783     ATM_qsol_beg      = ATM_qsol_beg   [ATM_beg_keysort]
784     ATM_qsol_end      = ATM_qsol_end   [ATM_end_keysort]
785     ATM_qs01_beg      = ATM_qs01_beg   [ATM_beg_keysort]
786     ATM_qs01_end      = ATM_qs01_end   [ATM_end_keysort]
787     ATM_qs02_beg      = ATM_qs02_beg   [ATM_beg_keysort]
788     ATM_qs02_end      = ATM_qs02_end   [ATM_end_keysort]
789     ATM_qs03_beg      = ATM_qs03_beg   [ATM_beg_keysort]
790     ATM_qs03_end      = ATM_qs03_end   [ATM_end_keysort]
791     ATM_qs04_beg      = ATM_qs04_beg   [ATM_beg_keysort]
792     ATM_qs04_end      = ATM_qs04_end   [ATM_end_keysort]
793     LIC_sno_beg       = LIC_sno_beg    [ATM_beg_keysort]
794     LIC_sno_end       = LIC_sno_end    [ATM_end_keysort]
795     LIC_runlic0_beg   = LIC_runlic0_beg[ATM_beg_keysort]
796     LIC_runlic0_end   = LIC_runlic0_end[ATM_end_keysort]
797   
798LIC_qs_beg = ATM_qs02_beg
799LIC_qs_end = ATM_qs02_end
800
801ATM_mas_sno_beg     = DYN_stock_int ( ATM_sno_beg  )
802ATM_mas_sno_end     = DYN_stock_int ( ATM_sno_end  )
803
804ATM_mas_qs_beg      = DYN_stock_int ( ATM_qs_beg   )
805ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   )
806ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg )
807ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end )
808ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg )
809ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end )
810ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg )
811ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end )
812ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg )
813ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end )
814ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg )
815ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end )
816
817LIC_mas_sno_beg     = DYN_stock_int ( LIC_sno_beg     )
818LIC_mas_sno_end     = DYN_stock_int ( LIC_sno_end     )
819LIC_mas_runlic0_beg = DYN_stock_int ( LIC_runlic0_beg )
820LIC_mas_runlic0_end = DYN_stock_int ( LIC_runlic0_end )
821
822LIC_mas_qs_beg  = ATM_mas_qs02_beg
823LIC_mas_qs_end  = ATM_mas_qs02_end
824LIC_mas_wat_beg = LIC_mas_qs_beg + LIC_mas_sno_beg
825LIC_mas_wat_end = LIC_mas_qs_end + LIC_mas_sno_end
826
827dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg
828dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg
829dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg
830
831dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg
832dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg
833dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg
834dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg
835
836dLIC_mas_qs      = LIC_mas_qs_end      - LIC_mas_qs_beg
837dLIC_mas_sno     = LIC_mas_sno_end     - LIC_mas_sno_beg
838dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg
839dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno
840
841echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' )
842echo ( '------------------------------------------------------------------------------------' )
843echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
844prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True  )
845
846echo ( f'\nChange of soil humidity content -- {Title} ' )
847echo ( '------------------------------------------------------------------------------------' )
848echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) )
849prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True )
850
851echo ( f'\nChange of atmosphere water+snow content -- {Title}  ' )
852echo ( '------------------------------------------------------------------------------------' )
853prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True)
854
855echo ( '\n====================================================================================' )
856echo ( f'-- SRF changes  -- {Title} ' )
857
858if Routing == 'SIMPLE' :
859    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   )
860    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   )
861    RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] )
862   
863    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
864    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
865    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
866   
867    RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   )
868    RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   )
869    RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] )
870   
871    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
872    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
873    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
874
875if Routing == 'SECHIBA' :
876    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
877    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
878    RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
879    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
880    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
881    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
882   
883    RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
884    RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
885    RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
886    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
887    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
888    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
889
890RUN_mas_wat_beg = RUN_mas_wat_fast_beg  + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + \
891                  RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg
892RUN_mas_wat_end = RUN_mas_wat_fast_end  + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + \
893                  RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end
894
895dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg
896dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg
897dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg
898dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg
899dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg
900dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg
901
902dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg
903
904echo ( f'\nRunoff reservoirs -- {Title} ')
905echo ( f'------------------------------------------------------------------------------------' )
906echo ( 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 ' )
907echo ( 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 ' )
908echo ( 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 ' )
909echo ( 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 ' )
910echo ( 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 ' )
911echo ( 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 ' )
912echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' )
913
914echo ( '------------------------------------------------------------------------------------' )
915
916prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True )
917prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True )
918prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True )
919prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True )
920prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True )
921prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True )
922prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True )
923
924echo ( f'\nWater content in routing  -- {Title} ' )
925echo ( '------------------------------------------------------------------------------------' )
926echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' )
927prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   )
928
929echo ( '\n====================================================================================' )
930print (f'Reading SRF restart')
931SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.)
932SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.)
933SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.)
934SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.)
935
936SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.)
937SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.)
938SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.)
939SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.)
940
941if LMDZ :
942    SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell_latlon')
943    SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell_latlon')
944    SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell_latlon')
945    SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell_latlon')
946    SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell_latlon')
947    SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell_latlon')
948    SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell_latlon')
949    SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell_latlon')
950 
951if ICO :
952    SRF_tot_watveg_beg  = SRF_tot_watveg_beg  [SRF_beg_keysort]
953    SRF_tot_watsoil_beg = SRF_tot_watsoil_beg [SRF_beg_keysort]
954    SRF_snow_beg        = SRF_snow_beg        [SRF_beg_keysort]
955    SRF_lakeres_beg     = SRF_lakeres_beg     [SRF_beg_keysort]
956    SRF_tot_watveg_end  = SRF_tot_watveg_end  [SRF_end_keysort]
957    SRF_tot_watsoil_end = SRF_tot_watsoil_end [SRF_end_keysort]
958    SRF_snow_end        = SRF_snow_end        [SRF_end_keysort]
959    SRF_lakeres_end     = SRF_lakeres_end     [SRF_end_keysort]
960
961# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood
962
963SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg
964SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end
965
966echo ( '\n====================================================================================' )
967print ('Computing integrals')
968
969print ( ' 1/8', end='' ) ; sys.stdout.flush ()
970SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  )
971print ( ' 2/8', end='' ) ; sys.stdout.flush ()
972SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg )
973print ( ' 3/8', end='' ) ; sys.stdout.flush ()
974SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        )
975print ( ' 4/8', end='' ) ; sys.stdout.flush ()
976SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     )
977print ( ' 5/8', end='' ) ; sys.stdout.flush ()
978
979SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  )
980print ( ' 6/8', end='' ) ; sys.stdout.flush ()
981SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end )
982print ( ' 7/8', end='' ) ; sys.stdout.flush ()
983SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        )
984print ( ' 8/8', end='' ) ; sys.stdout.flush ()
985SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     )
986
987print (' -- ') ; sys.stdout.flush ()
988
989dSRF_mas_watveg   = SRF_mas_watveg_end  - SRF_mas_watveg_beg
990dSRF_mas_watsoil  = SRF_mas_watsoil_end - SRF_mas_watsoil_beg
991dSRF_mas_snow     = SRF_mas_snow_end    - SRF_mas_snow_beg
992dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg
993
994echo ( '------------------------------------------------------------------------------------' )
995echo ( f'\nSurface reservoirs  -- {Title} ')
996echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' )
997echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' )
998echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' )
999echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' )
1000
1001prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True )
1002prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True )
1003prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True )
1004prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True )
1005
1006SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg
1007SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end
1008dSRF_mas_wat    = SRF_mas_wat_end - SRF_mas_wat_beg
1009
1010echo ( '------------------------------------------------------------------------------------' )
1011echo ( f'Water content in surface  -- {Title} ' )
1012echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ')
1013prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True )
1014
1015echo ( '------------------------------------------------------------------------------------' )
1016echo ( 'Water content in  ATM + SRF + RUN + LAKE' )
1017echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '.
1018           format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg ,
1019                   DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) )
1020prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True)
1021
1022echo ( '\n====================================================================================' )
1023echo ( f'-- ATM Fluxes  -- {Title} ' )
1024
1025print ( 'Step 1', end='' )
1026if ATM_HIS == 'latlon' :
1027    ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell_latlon' )
1028    ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic'], dim1D='cell_latlon' )
1029    ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter'], dim1D='cell_latlon' )
1030    ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic'], dim1D='cell_latlon' )
1031    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' )
1032    ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving'], dim1D='cell_latlon' )
1033    ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']  , dim1D='cell_latlon' )
1034    ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']   , dim1D='cell_latlon' )
1035    ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']     , dim1D='cell_latlon' )
1036    ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']     , dim1D='cell_latlon' )
1037    ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell_latlon' )
1038    ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell_latlon' )
1039    ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' )
1040    ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' )
1041    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' )
1042   
1043print ( ' - Step 2', end='' )
1044if ATM_HIS == 'ico' :
1045    ATM_wbilo_oce   = d_ATM_his ['wbilo_oce'][..., ATM_his_keysort]
1046    ATM_wbilo_sic   = d_ATM_his ['wbilo_sic'][..., ATM_his_keysort]
1047    ATM_wbilo_ter   = d_ATM_his ['wbilo_ter'][..., ATM_his_keysort]
1048    ATM_wbilo_lic   = d_ATM_his ['wbilo_lic'][..., ATM_his_keysort]
1049    ATM_runofflic   = d_ATM_his ['runofflic'][..., ATM_his_keysort]
1050    ATM_fqcalving   = d_ATM_his ['fqcalving'][..., ATM_his_keysort]
1051    ATM_fqfonte     = d_ATM_his ['fqfonte']  [..., ATM_his_keysort]
1052    ATM_precip      = d_ATM_his ['precip']   [..., ATM_his_keysort]
1053    ATM_snowf       = d_ATM_his ['snow']     [..., ATM_his_keysort]
1054    ATM_evap        = d_ATM_his ['evap']     [..., ATM_his_keysort]
1055    ATM_wevap_ter   = d_ATM_his ['wevap_ter'][..., ATM_his_keysort]
1056    ATM_wevap_oce   = d_ATM_his ['wevap_oce'][..., ATM_his_keysort]
1057    ATM_wevap_lic   = d_ATM_his ['wevap_lic'][..., ATM_his_keysort]
1058    ATM_wevap_sic   = d_ATM_his ['wevap_sic'][..., ATM_his_keysort]
1059    ATM_runofflic   = d_ATM_his ['runofflic'][..., ATM_his_keysort]
1060
1061print ( ' - Step 3', end='' )
1062ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic
1063ATM_emp         = ATM_evap - ATM_precip
1064ATM_wbilo_sea   = ATM_wbilo_oce + ATM_wbilo_sic
1065ATM_wevap_sea   = ATM_wevap_sic + ATM_wevap_oce
1066
1067ATM_wemp_ter = ATM_wevap_ter - ATM_precip * ATM_fter
1068ATM_wemp_oce = ATM_wevap_oce - ATM_precip * ATM_foce
1069ATM_wemp_sic = ATM_wevap_sic - ATM_precip * ATM_fsic
1070ATM_wemp_lic = ATM_wevap_lic - ATM_precip * ATM_flic
1071ATM_wemp_sea = ATM_wevap_sic + ATM_wevap_oce - ATM_precip * ATM_fsea
1072
1073print ( ' - Step 4', end='' )
1074if RUN_HIS == 'latlon' : 
1075    RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell_latlon' ) 
1076    RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']  , dim1D='cell_latlon' )
1077    RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']     , dim1D='cell_latlon' )
1078    RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']   , dim1D='cell_latlon' )
1079    RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']  , dim1D='cell_latlon' )
1080   
1081    RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'], dim1D='cell_latlon' ) 
1082    RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']  , dim1D='cell_latlon' )
1083
1084print ( ' - Step 5', end='' )
1085if RUN_HIS == 'ico' :
1086    RUN_coastalflow =  d_RUN_his ['coastalflow'][..., RUN_his_keysort]
1087    RUN_riverflow   =  d_RUN_his ['riverflow']  [..., RUN_his_keysort]
1088    RUN_runoff      =  d_RUN_his ['runoff']     [..., RUN_his_keysort]
1089    RUN_drainage    =  d_RUN_his ['drainage']   [..., RUN_his_keysort]
1090    RUN_riversret   =  d_RUN_his ['riversret']  [..., RUN_his_keysort]
1091   
1092    RUN_coastalflow_cpl = d_RUN_his ['coastalflow_cpl'][..., RUN_his_keysort]
1093    RUN_riverflow_cpl   = d_RUN_his ['riverflow_cpl']  [..., RUN_his_keysort]
1094
1095print ( ' - Step 6', end='' )
1096if SRF_HIS == 'latlon' :
1097    SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  , dim1D='cell_latlon')
1098    SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  , dim1D='cell_latlon')
1099    SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] , dim1D='cell_latlon')
1100    SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] , dim1D='cell_latlon')
1101    SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell_latlon' )
1102
1103print ( '- Step 7', end='' )
1104if SRF_HIS == 'ico' :
1105    SRF_rain     = d_SRF_his ['rain'] [..., SRF_his_keysort]
1106    SRF_evap     = d_SRF_his ['evap'] [..., SRF_his_keysort]
1107    SRF_snowf    = d_SRF_his ['snowf'][..., SRF_his_keysort]
1108    SRF_subli    = d_SRF_his ['subli'][..., SRF_his_keysort]
1109    SRF_transpir = np.sum (d_SRF_his ['transpir'], axis=1)[..., SRF_his_keysort]
1110
1111print ( '- Step 8', end='' )
1112SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
1113SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units']
1114   
1115## Correcting units of SECHIBA variables
1116def mmd2SI ( Var ) :
1117    '''Change unit from mm/d or m^3/s to kg/s if needed'''
1118    if 'units' in VarT.attrs : 
1119        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
1120            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s'
1121        if VarT.attrs['units'] == 'mm/d' :
1122            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
1123        if VarT.attrs['units'] in ['m^3', 'm3', ] :
1124            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg'
1125           
1126for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
1127    VarT = locals()['RUN_' + var]
1128    mmd2SI (VarT)
1129
1130for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
1131    VarT = locals()['SRF_' + var]
1132    mmd2SI (VarT)
1133
1134print ( '- Step 9', end='' )
1135RUN_input  = RUN_runoff      + RUN_drainage
1136RUN_output = RUN_coastalflow + RUN_riverflow
1137
1138print ( '- Step 10', end='' )
1139ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
1140ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
1141ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
1142ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
1143ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
1144ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
1145ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
1146ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
1147
1148print ( '- Step 11', end='' )
1149ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
1150ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
1151ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
1152ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
1153
1154print ( '- Step 12', end='' )
1155ATM_flx_wevap_ter   = ATM_flux_int ( ATM_wevap_ter )
1156ATM_flx_wevap_sea   = ATM_flux_int ( ATM_wevap_sea )
1157ATM_flx_precip_ter  = ATM_flux_int ( ATM_precip * ATM_fter )
1158ATM_flx_precip_sea  = ATM_flux_int ( ATM_precip * ATM_fsea )
1159ATM_flx_emp         = ATM_flux_int ( ATM_emp)
1160ATM_flx_wemp_ter    = ATM_flux_int ( ATM_wemp_ter )
1161ATM_flx_wemp_oce    = ATM_flux_int ( ATM_wemp_oce )
1162ATM_flx_wemp_lic    = ATM_flux_int ( ATM_wemp_lic )
1163ATM_flx_wemp_sea    = ATM_flux_int ( ATM_wemp_sea )
1164
1165print ( '- Step 13', end='' )
1166RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
1167RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
1168RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
1169RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
1170RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
1171RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
1172RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
1173RUN_flx_input       = SRF_flux_int ( RUN_input      )
1174RUN_flx_output      = ONE_flux_int ( RUN_output     )
1175
1176print ( '- Step 14' )
1177RUN_flx_bil    = RUN_flx_input   - RUN_flx_output
1178RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river
1179
1180prtFlux ('wbilo_oce     ', ATM_flx_wbilo_oce  , 'f' )         
1181prtFlux ('wbilo_sic     ', ATM_flx_wbilo_sic  , 'f' )         
1182prtFlux ('wbilo_sic+oce ', ATM_flx_wbilo_sea  , 'f' )         
1183prtFlux ('wbilo_ter     ', ATM_flx_wbilo_ter  , 'f' )         
1184prtFlux ('wbilo_lic     ', ATM_flx_wbilo_lic  , 'f' )         
1185prtFlux ('Sum wbilo_*   ', ATM_flx_wbilo      , 'f', True) 
1186prtFlux ('E-P           ', ATM_flx_emp        , 'f', True) 
1187prtFlux ('calving       ', ATM_flx_calving    , 'f' ) 
1188prtFlux ('fqfonte       ', ATM_flx_fqfonte    , 'f' )       
1189prtFlux ('precip        ', ATM_flx_precip     , 'f' )       
1190prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )       
1191prtFlux ('evap          ', ATM_flx_evap       , 'f' )
1192prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )   
1193
1194echo ( '\n====================================================================================' )
1195echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
1196prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
1197prtFlux ('riverflow     ', RUN_flx_river      , 'f' )       
1198prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
1199prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )   
1200prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
1201prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
1202prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
1203prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
1204prtFlux ('river in      ', RUN_flx_input      , 'f' )   
1205prtFlux ('river out     ', RUN_flx_output     , 'f' )   
1206prtFlux ('river bil     ', RUN_flx_bil        , 'f' )   
1207
1208ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river
1209echo ('')
1210#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 ))
1211
1212#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 ))
1213
1214ATM_flx_toSRF = -ATM_flx_wbilo_ter
1215
1216echo (' ')
1217echo ( '\n====================================================================================' )
1218echo ( f'--  Atmosphere  -- {Title} ' )
1219echo ( f'Mass begin = {DYN_mas_wat_beg:12.6e} kg | Mass end = {DYN_mas_wat_end:12.6e} kg' )
1220prtFlux ( 'dmass (atm)  = ', dDYN_mas_wat , 'e', True )
1221prtFlux ( 'Sum wbilo_*  = ', ATM_flx_wbilo, 'e', True )
1222prtFlux ( 'E-P          = ', ATM_flx_emp  , 'e', True )
1223echo ( ' ' )
1224prtFlux ( 'Water loss atm', ATM_flx_wbilo - dDYN_mas_wat, 'f', True )
1225echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) )
1226
1227echo (' ')
1228prtFlux ( 'Water loss atm', ATM_flx_emp  - dDYN_mas_wat , 'f', True )
1229echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) )
1230echo (' ')
1231
1232echo (' ')
1233echo ( '\n====================================================================================' )
1234
1235ATM_flx_runofflic_lic = ATM_flux_int ( ATM_runofflic * ATM_flic )
1236
1237LIC_flx_budget1 = -ATM_flx_wemp_lic  - ATM_flx_calving - ATM_flx_fqfonte
1238LIC_flx_budget2 = -ATM_flx_wbilo_lic - ATM_flx_calving - ATM_flx_fqfonte
1239LIC_flx_budget3 = -ATM_flx_wbilo_lic - ATM_flx_runofflic_lic
1240
1241echo ( f'--  LIC  -- {Title} ' )
1242echo ( f'Mass total   begin = {LIC_mas_wat_beg    :12.6e} kg | Mass end = {LIC_mas_wat_end    :12.6e} kg' )
1243echo ( f'Mass snow    begin = {LIC_mas_sno_beg    :12.6e} kg | Mass end = {LIC_mas_sno_end    :12.6e} kg' )
1244echo ( f'Mass qs      begin = {LIC_mas_qs_beg     :12.6e} kg | Mass end = {LIC_mas_qs_end     :12.6e} kg' )
1245echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' )
1246prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=65 )
1247prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=65 )
1248prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=65 )
1249prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=65 )
1250prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=65 )
1251prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=65 )
1252prtFlux ( 'LIC ATM_flx_fqfonte   ', ATM_flx_fqfonte       , 'f', True, width=65 )
1253prtFlux ( 'LIC ATM_flx_calving   ', ATM_flx_calving       , 'f', True, width=65 )
1254prtFlux ( 'LIC ATM_flx_runofflic ', ATM_flx_runofflic_lic , 'f', True, width=65 )
1255prtFlux ( 'LIC fqfonte + calving ', ATM_flx_calving+ATM_flx_fqfonte , 'f', True, width=65 )
1256prtFlux ( 'LIC fluxes 1 (wevap - precip*frac_lic  - fqcalving - fqfonte) ', LIC_flx_budget1              , 'f', True, width=65 )
1257prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)               ', LIC_flx_budget2              , 'f', True, width=65 )
1258prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)                ', LIC_flx_budget3              , 'f', True, width=65 )
1259prtFlux ( 'LIC error 1                                                   ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=65 )
1260prtFlux ( 'LIC error 2                                                   ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=65 )
1261prtFlux ( 'LIC error 3                                                   ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=65 )
1262echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) )
1263echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) )
1264echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) )
1265
1266echo ( '\n====================================================================================' )
1267echo ( f'-- SECHIBA fluxes  -- {Title} ' )
1268
1269SRF_flx_rain     = SRF_flux_int ( SRF_rain     )
1270SRF_flx_evap     = SRF_flux_int ( SRF_evap     )
1271SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    )
1272SRF_flx_subli    = SRF_flux_int ( SRF_subli    )
1273SRF_flx_transpir = SRF_flux_int ( SRF_transpir )
1274SRF_flx_emp      = SRF_flux_int ( SRF_emp      )
1275
1276RUN_flx_torouting  =  RUN_flx_runoff  + RUN_flx_drainage
1277RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river
1278
1279SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage
1280
1281prtFlux ('rain         ', SRF_flx_rain       , 'f' )
1282prtFlux ('evap         ', SRF_flx_evap       , 'f' )
1283prtFlux ('snowf        ', SRF_flx_snowf      , 'f' )
1284prtFlux ('E-P          ', SRF_flx_emp        , 'f' )
1285prtFlux ('subli        ', SRF_flx_subli      , 'f' )
1286prtFlux ('transpir     ', SRF_flx_transpir   , 'f' )
1287prtFlux ('to routing   ', RUN_flx_torouting  , 'f' )
1288prtFlux ('budget       ', SRF_flx_all        , 'f', small=True )
1289
1290echo ( '\n------------------------------------------------------------------------------------' )
1291echo ( 'Water content in surface ' )
1292echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' )
1293prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True)
1294prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True )
1295echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) )
1296
1297echo ( '\n====================================================================================' )
1298echo ( f'-- Check ATM vs. SRF -- {Title} ' )
1299prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' )
1300prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' )
1301prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' )
1302prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True)
1303echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) )
1304
1305echo ('')
1306echo ( '\n====================================================================================' )
1307echo ( f'-- RUNOFF fluxes  -- {Title} ' )
1308RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal
1309prtFlux ('runoff    ', RUN_flx_runoff      , 'f' ) 
1310prtFlux ('drainage  ', RUN_flx_drainage    , 'f' ) 
1311prtFlux ('run+drain ', RUN_flx_torouting   , 'f' ) 
1312prtFlux ('river     ', RUN_flx_river       , 'f' ) 
1313prtFlux ('coastal   ', RUN_flx_coastal     , 'f' ) 
1314prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' ) 
1315prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True)
1316
1317echo ( '\n------------------------------------------------------------------------------------' )
1318echo ( f'Water content in routing+lake  -- {Title} ' )
1319echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1320prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True)
1321prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True )
1322echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) )
1323
1324echo ( '\n------------------------------------------------------------------------------------' )
1325echo ( f'Water content in routing  -- {Title} ' )
1326echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1327prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  )
1328prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True)
1329echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) )
1330
1331echo ( ' ' )
1332echo ( f'{Title = }' )
1333
1334echo ( 'SVN Information' )
1335for clef in SVN.keys () : echo ( SVN[clef] )
Note: See TracBrowser for help on using the repository browser.