source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6446

Last change on this file since 6446 was 6277, checked in by omamce, 18 months ago

O.M. : WATER_BUDGET

Switch to .ini file to read configuratio
Correct SECHIBA budget
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 50.1 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
15#  $Author$
16#  $Date$
17#  $Revision$
18#  $Id$
19#  $HeadURL$
20
21###
22## Import system modules
23import sys, os, shutil, subprocess, platform
24import numpy as np
25import configparser, re
26
27## Creates parser
28config = configparser.ConfigParser()
29config.optionxform = str # To keep capitals
30
31## Read experiment parameters
32ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ;  TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
33
34##
35ARCHIVE =None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
36TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None
37file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
38file_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
39file_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
41d_ATM_his=None ; d_SRF_his=None ; d_RUN_his=None ; d_OCE_his=None ;  d_ICE_his=None ; d_OCE_sca=None
42d_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
43d_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
46print ( "Name of Python script:", sys.argv[0] )
47IniFile = sys.argv[1]
48print ("Input file : ", IniFile )
49config.read (IniFile)
50FullIniFile = 'full_' + IniFile
51
52def 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
59def 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
72def 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
81for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] :
82    if Section in config.keys() : 
83        print ( f'[{Section}]' )
84        for VarName in config[Section].keys() :
85            locals()[VarName] = config[Section][VarName]
86            locals()[VarName] = setBool (locals()[VarName])
87            locals()[VarName] = setNum  (locals()[VarName])
88            print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) )
89
90if not 'Files' in config.keys() : config['Files'] = {}
91
92def unDefined (char) :
93    if char in globals () :
94        if char == None : return True
95        else : return False
96    else : return True
97
98##-- Some physical constants
99#-- Earth Radius
100if not 'Ra'            in locals () : Ra = 6366197.7236758135
101#-- Gravity
102if not 'Grav'          in locals () : Grav = 9.81
103#-- Ice volumic mass (kg/m3) in LIM3
104if not 'ICE_rho_ice'   in locals () : ICE_rho_ice = 917.0
105#-- Snow volumic mass (kg/m3) in LIM3
106if not 'ICE_rho_sno'   in locals () : ICE_rho_sno = 330.0
107#-- Ocean water volumic mass (kg/m3) in NEMO
108if not 'OCE_rho_liq'   in locals () : OCE_rho_liq = 1026.
109#-- Water volumic mass in atmosphere
110if not 'ATM_rho'       in locals () : ATM_rho = 1000.
111#-- Water volumic mass in surface reservoirs
112if not 'SRF_rho'       in locals () : SRF_rho = 1000.
113#-- Water volumic mass of rivers
114if not 'RUN_rho'       in locals () : RUN_rho = 1000.
115#-- Year length
116if not 'YearLength'    in locals () : YearLength = 365.25 * 24. * 60. * 60.
117
118config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
119
120## --
121ICO  = ( 'ICO' in ATM )
122LMDZ = ( 'LMD' in ATM )
123   
124
125# Where do we run ?
126SysName, NodeName, Release, Version, Machine = os.uname ()
127TGCC  = ( 'irene'   in NodeName )
128IDRIS = ( 'jeanzay' in NodeName )
129
130## Set site specific libIGCM directories, and other specific stuff
131if 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       
153if IDRIS :
154    raise Exception ("Pour IDRIS : repertoires et chemins a definir") 
155
156config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU,
157                    'Program'  : "Generated by : " + str(sys.argv), 
158                    'HOSTNAME' : platform.node (), 'LOGNAME'  : os.getlogin (),
159                    'Python'   : f'{platform.python_version ()}',
160                    'OS'       : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}',
161                    'SVN_Author'   : "$Author$", 
162                    'SVN_Date'     : "$Date$",
163                    'SVN_Revision' : "$Revision$",
164                    'SVN_Id'       : "$Id$",
165                    'SVN_HeadURL'  : "$HeadURL$"}
166
167if libIGCM :
168    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
169
170## Import specific module
171import nemo, lmdz
172## Now import needed scientific modules
173import xarray as xr
174
175# Output file
176if FileOut == None : 
177    FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
178    config['Files']['FileOut'] = FileOut
179
180f_out = open ( FileOut, mode = 'w' )
181
182# Function to print to stdout *and* output file
183def echo (string, end='\n') :
184    print ( str(string), end=end  )
185    sys.stdout.flush ()
186    f_out.write ( str(string) + end )
187    f_out.flush ()
188    return None
189   
190## Set libIGCM directories
191R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT')
192R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT')
193
194L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName)
195R_SAVE      = os.path.join ( R_OUT, L_EXP )
196R_BUFR      = os.path.join ( R_BUF, L_EXP )
197POST_DIR    = os.path.join ( R_BUF, L_EXP, 'Out' )
198REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' )
199R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
200R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
201
202#if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir )
203if not os.path.isdir (TmpDir) : os.mkdir (TmpDir)
204TmpDirOCE = os.path.join (TmpDir, 'OCE')
205TmpDirICE = os.path.join (TmpDir, 'ICE')
206if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE )
207if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE )
208
209echo (' ')
210echo ( f'JobName : {JobName}' )
211echo (Comment)
212echo ( f'Working in TMPDIR : {TmpDir}' )
213
214echo ( f'\nDealing with {L_EXP}' )
215
216#-- Model output directories
217if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' )
218if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' )
219if dir_ATM_his == None :
220    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
221    config['Files']['dir_ATM_his'] = dir_ATM_his
222if dir_SRF_his == None : 
223    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
224    config['Files']['dir_SRF_his'] = dir_SRF_his
225   
226echo ( f'The analysis relies on files from the following model output directories : ' )
227echo ( f'{dir_ATM_his}' )
228echo ( f'{dir_SRF_his}' )
229
230#-- Files Names
231if Period == None :
232    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
233    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
234    config['Files']['Period'] = Period
235if FileCommon == None :
236    FileCommon = f'{JobName}_{Period}'
237    config['Files']['FileCommon'] = FileCommon
238
239if Title == None :
240    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
241    config['Files']['Title'] = Title
242     
243echo ('\nOpen history files' )
244if file_ATM_his == None : 
245    file_ATM_his = os.path.join (  dir_ATM_his, f'{FileCommon}_histmth.nc' )
246    config['Files']['file_ATM_his'] = file_ATM_his
247if file_SRF_his == None :
248    file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
249    config['Files']['file_SRF_his'] = file_SRF_his
250#if Routing == 'SECHIBA'  :
251#     file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
252if Routing == 'SIMPLE' :
253    if file_RUN_his == None : 
254        file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
255        config['Files']['file_RUN_his'] = file_RUN_his
256
257d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
258d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
259if Routing == 'SECHIBA' :
260    d_RUN_his = d_SRF_his
261if Routing == 'SIMPLE' : 
262    d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
263
264echo ( file_ATM_his )
265echo ( file_SRF_his )
266if Routing == 'SIMPLE' : 
267    echo ( file_RUN_his )
268
269## Compute run length
270dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
271echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
272dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
273
274## Compute length of each period
275dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
276echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
277dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
278dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
279
280## Number of years
281NbYear = dtime_sec / YearLength
282
283#-- Open restart files
284
285YearRes = YearBegin - 1              # Year of the restart of beginning of simulation
286YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
287
288if 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
293if 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
299if 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
302if 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   
306echo ( f'{file_restart_beg}' )
307echo ( f'{file_restart_end}' )
308
309if file_ATM_beg == None : 
310    file_ATM_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc'
311    config['Files']['file_ATM_beg'] = file_ATM_beg
312if file_ATM_end == None :
313    file_ATM_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc'
314    config['Files']['file_ATM_end'] = file_ATM_end
315
316liste_beg = [file_ATM_beg, ]
317liste_end = [file_ATM_end, ]
318   
319
320if file_DYN_beg == None :
321    if LMDZ :
322        file_DYN_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restart.nc'
323    if ICO  :
324        file_DYN_beg = f'{TmpDir}/ICO_{JobName}_{YearRes}1231_restart.nc'
325    liste_beg.append (file_DYN_beg)
326    config['Files']['file_DYN_beg'] = file_DYN_beg
327   
328if file_DYN_end == None : 
329    if LMDZ :
330        file_DYN_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restart.nc'
331    if ICO  :
332        file_DYN_end = f'{TmpDir}/ICO_{JobName}_{YearEnd}1231_restart.nc'
333    liste_end.append ( file_DYN_end )
334    config['Files']['file_DYN_end'] = file_DYN_end
335
336if file_SRF_beg == None : 
337    file_SRF_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'
338    config['Files']['file_SRF_beg'] = file_SRF_beg
339if file_SRF_end == None : 
340    file_SRF_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc'
341    config['Files']['file_SRF_end'] = file_SRF_end
342   
343liste_beg.append ( file_SRF_beg )
344liste_end.append ( file_SRF_end )
345
346echo ( f'{file_ATM_beg}')
347echo ( f'{file_ATM_end}')
348echo ( f'{file_DYN_beg}')
349echo ( f'{file_DYN_end}')
350echo ( f'{file_SRF_beg}')
351echo ( f'{file_SRF_end}')
352
353if Routing == 'SIMPLE' :
354    if file_RUN_beg == None : 
355        file_RUN_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc'
356        config['Files']['file_RUN_beg'] = file_RUN_beg
357    if file_RUN_end == None : 
358        file_RUN_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc'
359        config['Files']['file_RUN_end'] = file_RUN_end
360       
361    liste_beg.append ( file_RUN_beg )
362    liste_end.append ( file_RUN_end )
363    echo ( f'{file_RUN_beg}')
364    echo ( f'{file_RUN_end}')
365
366echo ('\nExtract restart files from tar : ATM, ICO and SRF')
367for resFile in liste_beg :
368    if os.path.exists ( os.path.join (TmpDir, resFile) ) :
369        echo ( f'file found : {resFile}' )
370    else :
371        base_resFile = os.path.basename (resFile)
372        command =  f'cd {TmpDir} ; tar xf {file_restart_beg} {base_resFile}'
373        echo ( command )
374        os.system ( command )
375       
376for 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       
385echo ('\nOpening ATM SRF and ICO restart files')
386d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze()
387d_ATM_end = xr.open_dataset ( os.path.join (TmpDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze()
388d_SRF_beg = xr.open_dataset ( os.path.join (TmpDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze()
389d_SRF_end = xr.open_dataset ( os.path.join (TmpDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze()
390d_DYN_beg = xr.open_dataset ( os.path.join (TmpDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze()
391d_DYN_end = xr.open_dataset ( os.path.join (TmpDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze()
392
393for var in d_SRF_beg.variables :
394    d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.)
395    d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.)
396
397if 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   
401echo ( file_ATM_beg )
402echo ( file_ATM_end )
403echo ( file_DYN_beg )
404echo ( file_DYN_end )
405echo ( file_SRF_beg )
406echo ( file_SRF_end )
407if Routing == 'SIMPLE' :
408    echo ( file_RUN_beg )
409    echo ( file_RUN_end )
410
411##
412config_out = open (FullIniFile, 'w')
413config.write (config_out )
414config_out.close ()
415
416def 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
432def 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
439def 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
447Psum = KSsum
448
449def 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
454def 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
459def 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
464def 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
469def ONE_stock_int (stock) :
470    '''Sum stock '''
471    ONE_stock_int  = Psum ( stock.to_masked_array().ravel() )
472    return ONE_stock_int
473
474def 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
479def kg2Sv  (val, rho=ATM_rho) :
480    '''From kg to Sverdrup'''
481    return val/dtime_sec*1.0e-6/rho
482
483def kg2myear (val, rho=ATM_rho) :
484    '''From kg to m/year'''
485    return val/ATM_aire_sea_tot/rho/NbYear
486   
487# ATM grid with cell surfaces
488if 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   
505if 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
512SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.)
513
514if ICO :
515    # Area on icosahedron grid
516    file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
517    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']
522    DYN_flnd = 1.0 - DYN_fsea
523    DYN_fter = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'})
524    DYN_flic = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'})
525   
526if LMDZ :
527    DYN_aire = ATM_aire
528    DYN_fsea = ATM_fsea
529    DYN_flnd = ATM_flnd
530    DYN_fter = d_ATM_beg['FTER']
531    DYN_flic = d_ATM_beg['FTER']
532
533   
534DYN_aire_fter = DYN_aire * DYN_fter
535
536#if LMDZ :
537#    d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} )
538
539ATM_aire_sea     = ATM_aire * ATM_fsea
540ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea )
541
542ATM_aire_tot     = ONE_stock_int (ATM_aire)
543ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea)
544
545
546echo ( 'Aire atmosphere/4pi R^2 : {:12.5f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )
547
548if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) :
549    raise Exception ('Erreur surface interpolee sur grille reguliere')
550
551echo ( '\n====================================================================================' )
552echo ( f'-- ATM changes in stores  -- {Title} ' )
553
554#-- Change in precipitable water from the atmosphere daily and monthly files
555#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
556
557# ATM vertical grid
558ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
559ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
560
561# Surface pressure
562if ICO : 
563    DYN_psol_beg = d_DYN_beg['ps']
564    DYN_psol_end = d_DYN_end['ps']
565if 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)
570DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg
571DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end
572
573klevp1 = ATM_Bhyb.shape[-1]
574if ICO  : cell_mesh        = DYN_psol_beg.shape[-1]
575if LMDZ : points_physiques = DYN_psol_beg.shape[-1]
576klev = klevp1 - 1
577
578# Layer thickness (pressure)
579if 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)        ) )
582if 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
586for k in np.arange (klevp1-1) :
587    DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:]
588    DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:]
589
590##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases
591if LMDZ :
592    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) ) )
595    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) ) )
598if ICO :
599    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'} )
602    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
607DYN_mas_wat_beg = ATM_stock_int (DYN_dpres_beg * DYN_wat_beg) / Grav
608DYN_mas_wat_end = ATM_stock_int (DYN_dpres_end * DYN_wat_end) / Grav
609
610# Variation of water content
611dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg
612
613# \([a-z,A-Z,_]*\)/dtime_sec\*1e-9  kg2Sv(\1)
614# \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear  kg2myear(\1)
615
616def var2prt (var, small=False) :
617    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000.
618    else     :  return var , kg2Sv(var)      , kg2myear(var)
619
620def 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
630echo ( f'\nVariation du contenu en eau atmosphere (dynamique)  -- {Title} ' )
631echo ( '------------------------------------------------------------------------------------' )
632echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
633prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True )
634
635ATM_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']
636ATM_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
638ATM_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']
639ATM_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
641ATM_qsol_beg = d_ATM_beg['QSOL']
642ATM_qsol_end = d_ATM_end['QSOL']
643
644ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
645ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER']
646ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
647ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC']
648ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
649ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE']
650ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
651ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
652
653if 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
669ATM_mas_sno_beg   = ATM_stock_int ( ATM_sno_beg  )
670ATM_mas_sno_end   = ATM_stock_int ( ATM_sno_end  )
671ATM_mas_qs_beg    = ATM_stock_int ( ATM_qs_beg   )
672ATM_mas_qs_end    = ATM_stock_int ( ATM_qs_end   )
673ATM_mas_qsol_beg  = ATM_stock_int ( ATM_qsol_beg )
674ATM_mas_qsol_end  = ATM_stock_int ( ATM_qsol_end )
675ATM_mas_qs01_beg  = ATM_stock_int ( ATM_qs01_beg )
676ATM_mas_qs01_end  = ATM_stock_int ( ATM_qs01_end )
677ATM_mas_qs02_beg  = ATM_stock_int ( ATM_qs02_beg )
678ATM_mas_qs02_end  = ATM_stock_int ( ATM_qs02_end )
679ATM_mas_qs03_beg  = ATM_stock_int ( ATM_qs03_beg )
680ATM_mas_qs03_end  = ATM_stock_int ( ATM_qs03_end )
681ATM_mas_qs04_beg  = ATM_stock_int ( ATM_qs04_beg )
682ATM_mas_qs04_end  = ATM_stock_int ( ATM_qs04_end )
683
684dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg
685dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg
686dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg
687
688dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg
689dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg
690dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg
691dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg
692
693dATM_mas_sno_2 = ATM_stock_int ( ATM_sno_end - ATM_sno_beg )
694
695echo ( f'\nVariation du contenu en neige atmosphere (calottes)  -- {Title} ' )
696echo ( '------------------------------------------------------------------------------------' )
697echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
698prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True  )
699prtFlux ( 'dMass (neige atm) ', dATM_mas_sno_2, 'e', True  )
700
701echo ( f'\nVariation du contenu humidite du sol  -- {Title} ' )
702echo ( '------------------------------------------------------------------------------------' )
703echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) )
704prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True )
705
706echo ( f'\nVariation du contenu en eau+neige atmosphere  -- {Title}  ' )
707echo ( '------------------------------------------------------------------------------------' )
708prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True)
709
710echo ( '\n====================================================================================' )
711echo ( f'-- SRF changes  -- {Title} ' )
712
713if Routing == 'SIMPLE' :
714    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   )
715    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   )
716    RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] )
717    RUN_mas_wat_flood_beg  = 0.0
718    RUN_mas_wat_lake_beg   = 0.0
719    RUN_mas_wat_pond_beg   = 0.0
720   
721    RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   )
722    RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   )
723    RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] )
724    RUN_mas_wat_flood_end  = 0.0
725    RUN_mas_wat_lake_end   = 0.0
726    RUN_mas_wat_pond_end   = 0.0
727
728if Routing == 'SECHIBA' :
729    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
730    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
731    RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
732    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
733    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
734    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
735
736   
737    RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
738    RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
739    RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
740    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
741    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
742    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
743
744RUN_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
745RUN_mas_wat_end = RUN_mas_wat_fast_end + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end
746
747dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg
748dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg
749dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg
750dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg
751dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg
752dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg
753
754dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg
755
756echo ( f'\nLes differents reservoirs  -- {Title} ')
757echo ( '------------------------------------------------------------------------------------' )
758echo ( '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   ) )
759echo ( '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   ) )
760echo ( '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 ) )
761echo ( '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  ) )
762echo ( '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   ) )
763echo ( 'RUN_mas_wat_pond_beg   = {:12.6e} kg | RUN_mas_wat_pond_end   = {:12.6e} kg '.format (RUN_mas_wat_pond_beg  , RUN_mas_wat_pond_end   ) )
764
765echo ( '------------------------------------------------------------------------------------' )
766
767prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True )
768prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True )
769prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True )
770prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True )
771prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True )
772prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True )
773prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True )
774
775echo ( f'\nWater content in routing  -- {Title} ' )
776echo ( '------------------------------------------------------------------------------------' )
777echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) )
778prtFlux ( 'dMass (routing) ', dRUN_mas_wat       )
779
780echo ( '\n====================================================================================' )
781print (f'Reading SRF restart')
782SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg < 1E15, 0.)
783SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg< 1E15, 0.)
784SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg       < 1E15, 0.)
785SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg    < 1E15, 0.)
786
787SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end < 1E15, 0.)
788SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end< 1E15, 0.)
789SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end       < 1E15, 0.)
790SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end    < 1E15, 0.)
791
792if 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    )
801 
802if 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'} )
811
812# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood
813
814SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg
815SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end
816
817echo ( '\n====================================================================================' )
818print ('Computing integrals')
819
820print ( ' 1/8', end='' ) ; sys.stdout.flush ()
821SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  )
822print ( ' 2/8', end='' ) ; sys.stdout.flush ()
823SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg )
824print ( ' 3/8', end='' ) ; sys.stdout.flush ()
825SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        )
826print ( ' 4/8', end='' ) ; sys.stdout.flush ()
827SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     )
828print ( ' 5/8', end='' ) ; sys.stdout.flush ()
829
830SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  )
831print ( ' 6/8', end='' ) ; sys.stdout.flush ()
832SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end )
833print ( ' 7/8', end='' ) ; sys.stdout.flush ()
834SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        )
835print ( ' 8/8', end='' ) ; sys.stdout.flush ()
836SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     )
837
838print (' -- ') ; sys.stdout.flush ()
839
840dSRF_mas_watveg   = SRF_mas_watveg_end  - SRF_mas_watveg_beg
841dSRF_mas_watsoil  = SRF_mas_watsoil_end - SRF_mas_watsoil_beg
842dSRF_mas_snow     = SRF_mas_snow_end    - SRF_mas_snow_beg
843dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg
844
845
846
847echo ( '------------------------------------------------------------------------------------' )
848echo ( f'\nLes differents reservoirs  -- {Title} ')
849echo ( 'SRF_mas_watveg_beg   = {:12.6e} kg | SRF_mas_watveg_end   = {:12.6e} kg '.format (SRF_mas_watveg_beg  , SRF_mas_watveg_end   ) )
850echo ( 'SRF_mas_watsoil_beg  = {:12.6e} kg | SRF_mas_watsoil_end  = {:12.6e} kg '.format (SRF_mas_watsoil_beg , SRF_mas_watsoil_end  ) )
851echo ( 'SRF_mas_snow_beg     = {:12.6e} kg | SRF_mas_snow_end     = {:12.6e} kg '.format (SRF_mas_snow_beg    , SRF_mas_snow_end     ) )
852echo ( 'SRF_mas_lake_beg     = {:12.6e} kg | SRF_mas_lake_end     = {:12.6e} kg '.format (SRF_mas_lake_beg    , SRF_mas_lake_end     ) )
853
854prtFlux ('dMass (watveg) ', dSRF_mas_watveg , 'e' , True)
855prtFlux ('dMass (watsoil)', dSRF_mas_watsoil, 'e' , True)
856prtFlux ('dMass (snow)   ', dSRF_mas_snow   , 'e' , True)
857prtFlux ('dMass (lake)   ', dSRF_mas_lake   , 'e' , True)
858
859SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg
860SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end
861dSRF_mas_wat    = SRF_mas_wat_end - SRF_mas_wat_beg
862
863echo ( '------------------------------------------------------------------------------------' )
864echo ( f'Water content in surface  -- {Title} ' )
865echo ( 'SRF_mas_wat_beg   = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) )
866prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True )
867
868echo ( '------------------------------------------------------------------------------------' )
869echo ( 'Water content in  ATM + SRF + RUN + LAKE' )
870echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '.
871           format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg ,
872                   DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) )
873prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True)
874
875echo ( '\n====================================================================================' )
876echo ( f'-- ATM Fluxes  -- {Title} ' )
877
878ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce']   )
879ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic']   )
880ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter']   )
881ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic']   )
882ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic']   )
883ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving']   )
884ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']     )
885ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']      )
886ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']        )
887ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']        )
888ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic
889ATM_emp         = ATM_evap - ATM_precip
890
891RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'] ) 
892RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']   )
893RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']      )
894RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']    )
895RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']   )
896
897RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'] ) 
898RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']   )
899
900SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  )
901SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  )
902SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] )
903SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] )
904SRF_transpir = lmdz.geo2point ( np.sum(d_SRF_his ['transpir'], axis=1) ) ; SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
905SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units']
906
907## Correcting units of SECHIBA variables
908def mmd2SI ( Var ) :
909    '''Change unit from mm/d or m^3/s to kg/s if needed'''
910    if 'units' in VarT.attrs : 
911        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-3'] :
912            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s'
913        if VarT.attrs['units'] == 'mm/d' :
914            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
915
916for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
917    VarT = locals()['RUN_' + var]
918    mmd2SI (VarT)
919
920for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
921    VarT = locals()['SRF_' + var]
922    mmd2SI (VarT)
923
924
925RUN_input  = RUN_runoff      + RUN_drainage
926RUN_output = RUN_coastalflow + RUN_riverflow
927
928ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic
929
930ATM_flx_oce         = ATM_flux_int ( ATM_wbilo_oce  )
931ATM_flx_sic         = ATM_flux_int ( ATM_wbilo_sic  )
932ATM_flx_sea         = ATM_flux_int ( ATM_wbilo_sea  )
933ATM_flx_ter         = ATM_flux_int ( ATM_wbilo_ter  )
934ATM_flx_lic         = ATM_flux_int ( ATM_wbilo_lic  )
935ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
936ATM_flx_qfonte      = ATM_flux_int ( ATM_fqfonte    )
937ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
938ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
939ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
940ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
941
942RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
943RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
944RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
945RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
946RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
947RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
948RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
949RUN_flx_input       = SRF_flux_int ( RUN_input      )
950RUN_flx_output      = ONE_flux_int ( RUN_output     )
951
952ATM_flx_emp   = ATM_flux_int (ATM_emp)
953ATM_flx_wbilo = ATM_flux_int (ATM_wbilo)
954
955RUN_flx_bil   = RUN_flx_input - RUN_flx_output
956RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river
957
958prtFlux ('wbilo_oce     ', ATM_flx_oce, 'f' )
959prtFlux ('wbilo_oce     ', ATM_flx_oce, 'f' )         
960prtFlux ('wbilo_sic     ', ATM_flx_sic, 'f' )         
961prtFlux ('wbilo_sic+oce ', ATM_flx_sea, 'f' )         
962prtFlux ('wbilo_ter     ', ATM_flx_ter, 'f' )         
963prtFlux ('wbilo_lic     ', ATM_flx_lic, 'f' )         
964prtFlux ('Sum wbilo_*   ', ATM_flx_wbilo      , 'f', True) 
965prtFlux ('E-P           ', ATM_flx_emp        , 'f', True) 
966prtFlux ('calving       ', ATM_flx_calving    , 'f' ) 
967prtFlux ('fqfonte       ', ATM_flx_qfonte     , 'f' )       
968prtFlux ('precip        ', ATM_flx_precip     , 'f' )       
969prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )       
970prtFlux ('evap          ', ATM_flx_evap       , 'f' )         
971prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
972prtFlux ('riverflow     ', RUN_flx_river      , 'f' )       
973prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
974prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )   
975prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
976prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
977prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
978prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
979prtFlux ('river in      ', RUN_flx_input      , 'f' )   
980prtFlux ('river out     ', RUN_flx_output     , 'f' )   
981prtFlux ('river bil     ', RUN_flx_bil        , 'f' )   
982prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )   
983
984ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river
985echo ('')
986echo ('  Global        {:12.3e} kg | {:12.4f} Sv | {:12.4f} m '.format ( ATM_flx_budget , ATM_flx_budget / dtime_sec*1E-9, ATM_flx_budget /ATM_aire_sea_tot/ATM_rho ))
987
988#echo ('  E-P-R         {:12.3e} kg | {:12.4e} Sv | {:12.4f} m '.format ( ATM_flx_emp      , ATM_flx_emp      / dtime_sec*1E-6/ATM_rho, ATM_flx_emp      /ATM_aire_sea_tot/ATM_rho ))
989
990echo (' ')
991echo ( '\n====================================================================================' )
992echo ( f'--  Atmosphere  -- {Title} ' )
993echo ( 'Mass begin = {:12.6e} kg | Mass end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
994echo ( '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. ))
995echo ( '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. ))
996echo ( '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. ))
997echo ( ' ' )
998prtFlux ( 'Perte eau atm ', ATM_flx_wbilo - dDYN_mas_wat, 'e', True )
999echo ( 'Perte eau atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) )
1000
1001echo (' ')
1002prtFlux ( 'Perte eau atm', ATM_flx_emp  - dDYN_mas_wat , 'e', True )
1003echo ( 'Perte eau atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) )
1004echo (' ')
1005
1006echo ( '\n====================================================================================' )
1007echo ( f'-- SECHIBA fluxes  -- {Title} ' )
1008
1009SRF_flx_rain     = SRF_flux_int ( SRF_rain     )
1010SRF_flx_evap     = SRF_flux_int ( SRF_evap     )
1011SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    )
1012SRF_flx_subli    = SRF_flux_int ( SRF_subli    )
1013SRF_flx_transpir = SRF_flux_int ( SRF_transpir )
1014SRF_flx_emp      = SRF_flux_int ( SRF_emp      )
1015
1016RUN_flx_torouting  =  RUN_flx_runoff  + RUN_flx_drainage
1017RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river
1018
1019SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage
1020
1021prtFlux ('rain       ', SRF_flx_rain     , 'f' )
1022prtFlux ('evap       ', SRF_flx_evap     , 'f' )
1023prtFlux ('snowf      ', SRF_flx_snowf    , 'f' )
1024prtFlux ('E-P        ', SRF_flx_emp      , 'f' )
1025prtFlux ('subli      ', SRF_flx_subli    , 'f' )
1026prtFlux ('transpir   ', SRF_flx_transpir , 'f' )
1027prtFlux ('to routing ', RUN_flx_torouting, 'f' )
1028prtFlux ('budget     ', SRF_flx_all      , 'f', True )
1029
1030echo ( '\n------------------------------------------------------------------------------------' )
1031echo ( 'Water content in surface ' )
1032echo ( 'SRF_mas_wat_beg  = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) )
1033prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True)
1034echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) )
1035
1036echo ('')
1037echo ( '\n====================================================================================' )
1038echo ( f'-- RUNOFF fluxes  -- {Title} ' )
1039RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal
1040prtFlux ('runoff    ', RUN_flx_runoff      , 'f' ) 
1041prtFlux ('drainage  ', RUN_flx_drainage    , 'f' ) 
1042prtFlux ('run+drain ', RUN_flx_torouting   , 'f' ) 
1043prtFlux ('river     ', RUN_flx_river       , 'f' ) 
1044prtFlux ('coastal   ', RUN_flx_coastal     , 'f' ) 
1045prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' ) 
1046prtFlux ('budget    ', RUN_flx_all         , 'f' , True)
1047
1048echo ( '\n------------------------------------------------------------------------------------' )
1049echo ( f'Water content in routing  -- {Title} ' )
1050echo ( 'RUN_mas_wat_beg  = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) )
1051prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', True)
1052
1053echo ( '\n------------------------------------------------------------------------------------' )
1054echo ( f'Water content in routing  -- {Title} ' )
1055echo ( 'RUN_mas_wat_beg  = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) )
1056prtFlux ( 'dMass (routing)  ', dRUN_mas_wat, 'f', True  )
1057
1058print ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) )
1059
1060# PRINT *,"routing water Balance ;  before : ", water_balance_before," ; after : ",water_balance_after,  &a
1061# 1150              " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%"
Note: See TracBrowser for help on using the repository browser.