source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6676

Last change on this file since 6676 was 6676, checked in by omamce, 8 months ago

O.M. : WATER_BUDGET

Code restructuration
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 64.8 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
26import os
27import configparser
28
29## Import needed scientific modules
30import numpy as np
31import xarray as xr
32
33## Import local modules
34import WaterUtils as wu
35import lmdz
36
37## Read command line arguments
38## ---------------------------
39print ( "Name of Python script:", sys.argv[0] )
40IniFile = sys.argv[1]
41
42# Test existence of IniFile
43if not os.path.exists (IniFile ) :
44    raise FileExistsError ( f"File not found : {IniFile = }" )
45
46if 'full' in IniFile or 'ATM' in IniFile :
47    FullIniFile = IniFile
48else :
49    FullIniFile = 'ATM_' + IniFile
50
51print ("Output file : ", FullIniFile )
52
53## Experiment parameters
54## --------------------
55dpar = wu.ReadConfig ( IniFile )
56
57## Configure all needed parameter from existant parameters
58## -------------------------------------------------------
59dpar = wu.SetDatesAndFiles ( dpar )
60
61## Output file with water budget diagnostics
62## -----------------------------------------
63f_out = dpar['Files']['f_out']
64
65## Put dpar values in local namespace
66## ----------------------------------
67for Section in dpar.keys () : 
68    print ( f'\nReading [{Section}]' )
69    for VarName in dpar[Section].keys() :
70        globals()[VarName] = dpar[Section][VarName]
71        print ( f'    {VarName:21} set to : {globals()[VarName]}' )
72
73## Debuging and timer
74Timer = wu.functools.partial (wu.Timer, debug=Debug, timer=Timing)
75
76## Useful functions
77## ----------------
78if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float :
79    def rprec (ptab) :
80        '''This version does nothing
81
82        rprec may be used to reduce floating precision when reading history files
83        '''
84        return ptab
85else                 :
86    def rprec (ptab) :
87        '''Returns float with a different precision'''
88        return ptab.astype(readPrec).astype(float)
89
90def kg2Sv    (pval, rho=ATM_RHO) :
91    '''From kg to Sverdrup'''
92    return pval/dtime_sec*1.0e-6/rho
93
94def kg2myear (pval, rho=ATM_RHO) :
95    '''From kg to m/year'''
96    return pval/ATM_aire_sea_tot/rho/NbYear
97
98def var2prt (pvar, small=False, rho=ATM_RHO) :
99    '''Formats value for printing'''
100    if small :  return pvar, kg2Sv (pvar, rho=rho)*1000., kg2myear (pvar, rho=rho)*1000
101    else     :  return pvar, kg2Sv (pvar, rho=rho)      , kg2myear (pvar, rho=rho)
102
103def prtFlux (Desc, pvar, Form='F', small=False, rho=ATM_RHO, width=15) :
104    '''Pretty print of formatted value'''
105    if small :
106        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
107        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
108    else :
109        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
110        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
111    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt (pvar, small=small, rho=rho), width=width ) )
112
113def echo (string, end='\n') :
114    '''Function to print to stdout *and* output file'''
115    print ( str(string), end=end  )
116    sys.stdout.flush ()
117    f_out.write ( str(string) + end )
118    f_out.flush ()
119
120## Open history files
121## ------------------
122d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
123if SRF :
124    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
125    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
126    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
127
128echo ( f'{file_ATM_his = }' )
129if SRF :
130    echo ( f'{file_SRF_his = }' )
131   
132## Compute run length
133## ------------------
134dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
135echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' )
136dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
137
138##-- Compute length of each period
139dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
140echo ( '\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
141dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
142dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
143dtime_per_sec.attrs['unit'] = 's'
144
145##-- Number of years (approximative)
146NbYear = dtime_sec / YEAR_LENGTH
147
148## Define restart periods and file names
149## -------------------------------------
150dpar = wu.SetRestartNames ( dpar, f_out) 
151
152## Put dpar values in local namespace
153## ----------------------------------
154for Section in dpar.keys () : 
155    print ( f'\nReading [{Section}]' )
156    for VarName in dpar[Section].keys() :
157        locals()[VarName] = dpar[Section][VarName]
158        print ( f'    {VarName:21} set to : {locals()[VarName]}' )
159       
160## Extract restart files from tar
161## ----------------------------------
162
163liste_beg = [file_ATM_beg, ]
164liste_end = [file_ATM_end, ]
165if file_DYN_beg : liste_beg.append ( file_DYN_beg )
166if file_DYN_end : liste_end.append ( file_DYN_end )
167if file_SRF_beg : liste_beg.append ( file_SRF_beg )
168if file_SRF_end : liste_end.append ( file_SRF_end )
169
170if ICO :
171    if not file_DYN_aire :
172        file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ResolAtm+'_grid.nc' )
173    dpar['Files']['file_DYN_aire'] = file_DYN_aire
174
175if SRF and Routing == 'SIMPLE' :
176    liste_beg.append ( file_RUN_beg )
177    liste_end.append ( file_RUN_end )
178
179echo ( '\nExtract restart files from tar : ATM, ICO', end='')
180if SRF : echo ( ' and SRF')
181else   : echo (' ')
182
183## Write the full configuration
184## ----------------------------
185params_out = open (FullIniFile, 'w', encoding = 'utf-8')
186params = wu.dict2config ( dpar )
187params.write ( params_out )
188params_out.close ()
189   
190@Timer
191def extract ( file_name=file_ATM_beg, tar_restart=tar_restart_end, file_dir_comp=FileDir ) :
192    '''
193    Extract restart files from a tar
194    '''
195    echo ( f'----------')
196    error_count = 0
197    echo ( f'file to extract : {file_name = }' )
198    if os.path.exists ( os.path.join (FileDir, file_name) ) :
199        echo ( f'file found      : {file_name = }' )
200    else :
201        echo ( f'file not found  : {file_name = }' )
202        base_resFile = os.path.basename (file_name)
203        if os.path.exists ( tar_restart ) :
204            command =  f'cd {file_dir_comp} ; tar xf {tar_restart} {base_resFile}'
205            echo ( f'{command = }' )
206            err = os.system ( command )
207            if err != 0 :
208                if ContinueOnError :
209                    error_count += 1
210                    echo ( f'****** Command failed : {command}' )
211                    echo ( '****** Trying to continue' )
212                    echo ( ' ')
213                else :
214                    raise OSError ( f'**** command failed : {command} - Stopping' )
215            else :
216                echo ( f'tar done : {base_resFile}' )
217        else :
218            echo ( f'****** Tar restart file {tar_restart = } not found ' )
219            if ContinueOnError :
220                error_count += 1
221                echo ( f'****** Command failed : {command}' )
222                echo (  '****** Trying to continue' )
223                echo (  ' ')
224            else :
225                raise OSError ( f'****** tar file not found {tar_restart = } - Stopping' )
226    return error_count
227
228ErrorCount = 0
229
230ErrorCount += extract ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, file_dir_comp=FileDir )
231ErrorCount += extract ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, file_dir_comp=FileDir )
232
233ErrorCount += extract ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, file_dir_comp=FileDir )
234ErrorCount += extract ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, file_dir_comp=FileDir )
235
236if SRF :
237    ErrorCount += extract ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, file_dir_comp=FileDir )
238    ErrorCount += extract ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, file_dir_comp=FileDir )
239
240    if Routing == 'SIMPLE' :
241        ErrorCount += extract ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, file_dir_comp=FileDir )
242        ErrorCount += extract ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, file_dir_comp=FileDir )
243
244##-- Exit in case of error in the opening file phase
245if ErrorCount > 0 :
246    echo ( '  ' )
247    raise RuntimeError ( f'**** Some files missing - Stopping - {ErrorCount = }' )
248
249##
250echo ('\nOpening ATM SRF and ICO restart files')
251d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
252d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
253if SRF :
254    d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
255    d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
256d_DYN_beg = xr.open_dataset ( os.path.join (FileDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
257d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
258
259if SRF :
260    for var in d_SRF_beg.variables :
261        d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.)
262        d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.)
263
264    if Routing == 'SIMPLE' :
265        d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
266        d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
267
268@Timer
269def to_cell ( dd, newname='cell' ) :
270    '''Set space dimension to newname
271    '''
272    for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] :
273        if oldname in dd.dims and oldname != newname :
274            dd = dd.rename ( {oldname : newname} )
275    return dd
276
277d_ATM_beg = to_cell ( d_ATM_beg )
278d_ATM_end = to_cell ( d_ATM_end )
279if SRF :
280    d_SRF_beg = to_cell ( d_SRF_beg )
281    d_SRF_end = to_cell ( d_SRF_end )
282d_DYN_beg = to_cell ( d_DYN_beg )
283d_DYN_end = to_cell ( d_DYN_end )
284
285if SRF and Routing == 'SIMPLE' :
286    d_RUN_beg = to_cell ( d_RUN_beg )
287    d_RUN_end = to_cell ( d_RUN_end )
288
289d_ATM_his = to_cell ( d_ATM_his )
290if SRF : d_SRF_his = to_cell ( d_SRF_his )
291
292echo ( f'{file_ATM_beg = }' )
293echo ( f'{file_ATM_end = }' )
294echo ( f'{file_DYN_beg = }' )
295echo ( f'{file_DYN_end = }' )
296if SRF :
297    echo ( f'{file_SRF_beg = }' )
298    echo ( f'{file_SRF_end = }' )
299    if Routing == 'SIMPLE' :
300        echo ( f'{file_RUN_beg = }' )
301        echo ( f'{file_RUN_end = }' )
302
303# ATM grid with cell surfaces
304if LMDZ :
305    echo ('ATM grid with cell surfaces : LMDZ')
306    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
307    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
308    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumul_poles=True, dim1d='cell' )
309    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
310    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
311    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
312    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
313    if SRF :
314        SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
315        SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
316        SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1d='cell', cumul_poles=True )
317        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1d='cell', cumul_poles=True )
318        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1d='cell' )
319
320if ICO :
321    if ATM_HIS == 'latlon' :
322        echo ( 'ATM areas and fractions on LATLON grid' )
323        if 'lat_dom_out' in d_ATM_his.variables :
324            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
325            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
326        else :
327            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
328            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
329        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumul_poles=True, dim1d='cell' )
330        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
331        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
332        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
333        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
334
335    if ATM_HIS == 'ico' :
336        echo ( 'ATM areas and fractions on ICO grid' )
337        ATM_aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
338        ATM_lat  =  rprec (d_ATM_his ['lat']         )
339        ATM_lon  =  rprec (d_ATM_his ['lon']         )
340        ATM_fter =  rprec (d_ATM_his ['fract_ter'][0])
341        ATM_foce =  rprec (d_ATM_his ['fract_oce'][0])
342        ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0])
343        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0])
344
345    if SRF :
346        if SRF_HIS == 'latlon' :
347            echo ( 'SRF areas and fractions on LATLON grid' )
348            if 'lat_domain_landpoints_out' in d_SRF_his  :
349                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
350                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
351            else :
352                if 'lat_domain_landpoints_out' in d_SRF_his :
353                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
354                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
355                else :
356                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
357                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
358
359            SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1d='cell', cumul_poles=True )
360            SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1d='cell', cumul_poles=True )
361            SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1d='cell', cumul_poles=True )
362            SRF_aire = SRF_areafrac
363
364        if SRF_HIS == 'ico' :
365            echo ( 'SRF areas and fractions on ICO grid' )
366            SRF_lat       =  rprec (d_SRF_his ['lat']     )
367            SRF_lon       =  rprec (d_SRF_his ['lon']     )
368            SRF_areas     =  rprec (d_SRF_his ['Areas']   )
369            SRF_contfrac  =  rprec (d_SRF_his ['Contfrac'])
370            SRF_aire      =  SRF_areas * SRF_contfrac
371
372ATM_fsea      = ATM_foce + ATM_fsic
373ATM_flnd      = ATM_fter + ATM_flic
374ATM_aire_fter = ATM_aire * ATM_fter
375ATM_aire_flic = ATM_aire * ATM_flic
376ATM_aire_fsic = ATM_aire * ATM_fsic
377ATM_aire_foce = ATM_aire * ATM_foce
378ATM_aire_flnd = ATM_aire * ATM_flnd
379ATM_aire_fsea = ATM_aire * ATM_fsea
380
381#SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. )
382
383if ICO :
384    # Area on icosahedron grid
385    echo ( f'{file_DYN_aire = }' )
386    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
387   
388    if SortIco :
389            # Creation d'une clef de tri pour le fichier aire
390            DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) )
391    else :
392            DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) )
393           
394    DYN_lat = d_DYN_aire['lat']
395    DYN_lon = d_DYN_aire['lon']
396       
397    DYN_aire = d_DYN_aire['aire']
398    DYN_fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
399       
400    DYN_flnd = 1.0 - DYN_fsea
401    DYN_fter = d_ATM_beg['FTER']
402    DYN_flic = d_ATM_beg['FLIC']
403    DYN_foce = d_ATM_beg['FOCE']
404    DYN_aire_fter = DYN_aire * DYN_fter
405
406#if ATM_HIS == 'ico' :
407#    ATM_aire      = DYN_aire
408#    ATM_aire_fter = DYN_aire * ATM_fter
409#    ATM_aire_flic = DYN_aire * ATM_flic
410#    ATM_aire_fsic = DYN_aire * ATM_fsic
411#    ATM_aire_foce = DYN_aire * ATM_foce
412#    ATM_aire_flnd = DYN_aire * ATM_flnd
413#    ATM_aire_fsea = DYN_aire * ATM_fsea
414
415if ATM_HIS == 'ico' :
416    DYN_aire      = ATM_aire
417    DYN_foce      = ATM_foce
418    DYN_fsic      = ATM_fsic
419    DYN_flic      = ATM_flic
420    DYN_fter      = ATM_fter
421    DYN_fsea      = ATM_fsea
422    DYN_flnd      = ATM_flnd
423    DYN_aire_fter = ATM_aire_fter
424    DYN_aire_flic = ATM_aire_flic
425    DYN_aire_fsic = ATM_aire_fsic
426    DYN_aire_foce = ATM_aire_foce
427    DYN_aire_flnd = ATM_aire_flnd
428    DYN_aire_fsea = ATM_aire_fsea
429
430if LMDZ :
431    # Area on lon/lat grid
432    DYN_aire = ATM_aire
433    DYN_fsea = ATM_fsea
434    DYN_flnd = ATM_flnd
435    DYN_fter = rprec (d_ATM_beg['FTER'])
436    DYN_flic = rprec (d_ATM_beg['FLIC'])
437    DYN_aire_fter = DYN_aire * DYN_fter
438
439if ICO and ATM_HIS == 'ico' :
440    # Comparaison des aires sur ATM et DYN
441    aire_diff = ATM_aire - DYN_aire
442    echo ( 'f{Difference Aire hist file vs. grid file {aire_diff.mean()=} {aire_diff.min()=}  {aire_diff.max()=} ' )
443   
444
445# Functions computing integrals and sum
446@Timer
447def DYN_stock_int (stock) :
448    '''Integrate (* surface) stock on atmosphere grid'''
449    return wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() )
450
451@Timer
452def ATM_flux_int (flux) :
453    '''Integrate (* time * surface) flux on atmosphere grid'''
454    return wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() )
455
456@Timer
457def LIC_flux_int (flux) :
458    '''Integrate (* time * surface) flux on land ice grid'''
459    return wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() )
460
461if SRF :
462    @Timer
463    def SRF_stock_int (stock) :
464        '''Integrate (* surface) stock on land grid'''
465        return wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
466
467    @Timer
468    def SRF_flux_int (flux) :
469        '''Integrate (* time * surface) flux on land grid'''
470        return wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() )
471
472@Timer
473def ONE_stock_int (stock) :
474    '''Sum stock'''
475    return wu.Psum ( stock.to_masked_array().ravel() )
476
477@Timer
478def ONE_flux_int (flux) :
479    '''Integrate (* time) flux on area=1 grid'''
480    return wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() )
481
482@Timer
483def Sprec ( tlist ) :
484    '''Accurate sum of list of scalar elements'''
485    return wu.Psum ( np.array ( tlist) )
486
487ATM_aire_sea     = ATM_aire * ATM_fsea
488
489ATM_aire_tot     = ONE_stock_int (ATM_aire)
490ATM_aire_sea_tot = ONE_stock_int (ATM_aire_fsea)
491ATM_aire_ter_tot = ONE_stock_int (ATM_aire_fter)
492ATM_aire_lic_tot = ONE_stock_int (ATM_aire_flic)
493
494DYN_aire_tot     = ONE_stock_int ( DYN_aire )
495if SRF : SRF_aire_tot     = ONE_stock_int ( SRF_aire )
496
497echo ('')
498echo ( f'ATM DYN : Area of atmosphere        : {DYN_aire_tot:18.9e}' )
499echo ( f'ATM HIS : Area of atmosphere        : {ATM_aire_tot:18.9e}' )
500echo ( f'ATM HIS : Area of ter in atmosphere : {ATM_aire_ter_tot:18.9e}' )
501echo ( f'ATM HIS : Area of lic in atmosphere : {ATM_aire_lic_tot:18.9e}' )
502if SRF :
503    echo ( f'ATM SRF : Area of atmosphere        : {SRF_aire_tot:18.9e}' )
504echo ('')
505echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN_aire_tot/(RA*RA*4*np.pi) ) )
506echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM_aire_tot/(RA*RA*4*np.pi) ) )
507echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_ter_tot/(RA*RA*4*np.pi) ) )
508if SRF :
509    echo ( 'ATM SRF : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF_aire_tot/(RA*RA*4*np.pi) ) )
510    echo ('')
511    echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' )
512
513
514if np.abs (ATM_aire_tot/(RA*RA*4*np.pi) - 1.0) > 0.01 :
515    raise RuntimeError ('Error of atmosphere surface interpolated on lon/lat grid')
516
517echo ( '\n====================================================================================' )
518echo ( f'-- ATM changes in stores  -- {Title} ' )
519
520#-- Change in precipitable water from the atmosphere daily and monthly files
521#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
522
523echo ( 'ATM vertical grid' )
524ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
525ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
526
527echo ( 'Surface pressure' )
528if ICO :
529    DYN_psol_beg = d_DYN_beg['ps']
530    DYN_psol_end = d_DYN_end['ps']
531if LMDZ :
532    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' )
533    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' )
534
535echo ( '3D Pressure at the interface layers (not scalar points)' )
536DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg
537DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end
538
539echo ( 'Check computation of pressure levels' )
540
541ind = np.empty (8)
542
543ind[0] = (DYN_pres_beg[ 0]-DYN_psol_beg).min().item()
544ind[1] = (DYN_pres_beg[ 0]-DYN_psol_beg).max().item()
545ind[2] = (DYN_pres_beg[-1]).min().item()
546ind[3] = (DYN_pres_beg[-1]).max().item()
547ind[4] = (DYN_pres_end[ 0]-DYN_psol_end).min().item()
548ind[5] = (DYN_pres_end[ 0]-DYN_psol_end).max().item()
549ind[6] = (DYN_pres_end[-1]).min().item()
550ind[7] = (DYN_pres_end[-1]).max().item()
551
552if any ( ind != 0) :
553    echo ( 'All values should be zero' )
554    echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).min().item() = {ind[0]}' )
555    echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).max().item() = {ind[1]}' )
556    echo ( f'(DYN_pres_beg[-1]).min().item()              = {ind[2]}' )
557    echo ( f'(DYN_pres_beg[-1]).max().item()              = {ind[3]}' )
558    echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).min().item() = {ind[4]}' )
559    echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).max().item() = {ind[5]}' )
560    echo ( f'(DYN_pres_end[-1]).min().item()              = {ind[6]}' )
561    echo ( f'(DYN_pres_end[-1]).max().item()              = {ind[7]}' )
562    raise RuntimeError
563
564klevp1 = ATM_Bhyb.shape[-1]
565cell   = DYN_psol_beg.shape[-1]
566klev   = klevp1 - 1
567
568echo ( 'Layer thickness (pressure)' )
569DYN_mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
570DYN_mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
571
572for k in np.arange (klev) :
573    DYN_mass_beg[k,:] = ( DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] ) / GRAV
574    DYN_mass_end[k,:] = ( DYN_pres_end[k,:] - DYN_pres_end[k+1,:] ) / GRAV
575
576DYN_mass_beg_2D = DYN_mass_beg.sum (dim='sigs')
577DYN_mass_end_2D = DYN_mass_end.sum (dim='sigs')
578
579DYN_mas_air_beg = DYN_stock_int ( DYN_mass_beg_2D )
580DYN_mas_air_end = DYN_stock_int ( DYN_mass_end_2D )
581
582echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' )
583if LMDZ :
584    if 'H2Ov' in d_DYN_beg.variables :
585        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' )
586        DYN_wat_beg = lmdz.geo3point ( d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' )
587        DYN_wat_end = lmdz.geo3point ( d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' )
588    if 'H2Ov_g' in d_DYN_beg.variables :
589        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' )
590        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' )
591        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' )
592if ICO :
593    if 'H2Ov_g' in d_DYN_beg.variables :
594        echo ('reading ICO : H2Ov_g, H2Ov_l, H2Ov_s' )
595        DYN_wat_beg = d_DYN_beg['H2Ov_g'] + d_DYN_beg['H2Ov_l'] + d_DYN_beg['H2Ov_s']
596        DYN_wat_end = d_DYN_end['H2Ov_g'] + d_DYN_end['H2Ov_l'] + d_DYN_end['H2Ov_s']
597    elif 'H2O_g' in d_DYN_beg.variables :
598        echo ('reading ICO : H2O_g, H2O_l, H2O_s' )
599        DYN_wat_beg = d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l']   + d_DYN_beg['H2O_s']
600        DYN_wat_end = d_DYN_end['H2O_g'] + d_DYN_end['H2O_l']   + d_DYN_end['H2O_s']
601    elif 'q' in d_DYN_beg.variables :
602        echo ('reading ICO : q' )
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)
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)
605
606if 'lev' in DYN_wat_beg.dims :
607    DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} )
608    DYN_wat_end = DYN_wat_end.rename ( {'lev':'sigs'} )
609
610echo ( 'Compute water content : vertical and horizontal integral' )
611
612DYN_wat_beg_2D =  (DYN_mass_beg * DYN_wat_beg).sum (dim='sigs')
613DYN_wat_end_2D =  (DYN_mass_end * DYN_wat_end).sum (dim='sigs')
614
615DYN_mas_wat_beg = DYN_stock_int ( DYN_wat_beg_2D )
616DYN_mas_wat_end = DYN_stock_int ( DYN_wat_end_2D )
617
618echo ( 'Variation of water content' )
619dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg
620
621echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' )
622echo ( '------------------------------------------------------------------------------------' )
623echo ( 'DYN_mas_air_beg = {:12.6e} kg | DYN_mas_air_end = {:12.6e} kg'.format (DYN_mas_air_beg, DYN_mas_air_end) )
624echo ( 'DYN_mas_wat_beg = {:12.6e} kg | DYN_mas_wat_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
625prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True )
626
627ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \
628              d_ATM_beg['SNOW03']*d_ATM_beg['FOCE'] + d_ATM_beg['SNOW04']*d_ATM_beg['FSIC']
629ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \
630              d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC']
631
632ATM_qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \
633              d_ATM_beg['QS03']  *d_ATM_beg['FOCE'] + d_ATM_beg['QS04']  *d_ATM_beg['FSIC']
634ATM_qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \
635              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC']
636
637ATM_qsol_beg     = d_ATM_beg['QSOL']
638ATM_qsol_end     = d_ATM_end['QSOL']
639
640LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']
641LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC']
642
643LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0']
644LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0']
645
646ATM_qs01_beg     = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
647ATM_qs02_beg     = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
648ATM_qs03_beg     = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
649ATM_qs04_beg     = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
650
651ATM_qs01_end     = d_ATM_end['QS01'] * d_ATM_end['FTER']
652ATM_qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC']
653ATM_qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE']
654ATM_qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC']
655
656LIC_qs_beg = ATM_qs02_beg
657LIC_qs_end = ATM_qs02_end
658
659ATM_mas_sno_beg     = DYN_stock_int ( ATM_sno_beg  )
660ATM_mas_sno_end     = DYN_stock_int ( ATM_sno_end  )
661
662ATM_mas_qs_beg      = DYN_stock_int ( ATM_qs_beg   )
663ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   )
664ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg )
665ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg )
666ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg )
667ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg )
668ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg )
669ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end )
670ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end )
671ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end )
672ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end )
673ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end )
674
675LIC_mas_sno_beg     = DYN_stock_int ( LIC_sno_beg     )
676LIC_mas_sno_end     = DYN_stock_int ( LIC_sno_end     )
677LIC_mas_runlic0_beg = DYN_stock_int ( LIC_runlic0_beg )
678LIC_mas_runlic0_end = DYN_stock_int ( LIC_runlic0_end )
679
680LIC_mas_qs_beg  = ATM_mas_qs02_beg
681LIC_mas_qs_end  = ATM_mas_qs02_end
682LIC_mas_wat_beg = LIC_mas_qs_beg + LIC_mas_sno_beg
683LIC_mas_wat_end = LIC_mas_qs_end + LIC_mas_sno_end
684
685dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg
686dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg
687dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg
688
689dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg
690dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg
691dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg
692dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg
693
694dLIC_mas_qs      = LIC_mas_qs_end      - LIC_mas_qs_beg
695dLIC_mas_sno     = LIC_mas_sno_end     - LIC_mas_sno_beg
696dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg
697
698dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno # + dLIC_mas_runlic0
699
700echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' )
701echo ( '------------------------------------------------------------------------------------' )
702echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
703prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True  )
704
705echo ( f'\nChange of soil humidity content -- {Title} ' )
706echo ( '------------------------------------------------------------------------------------' )
707echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) )
708prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True )
709
710echo ( f'\nChange of atmosphere water+snow content -- {Title}  ' )
711echo ( '------------------------------------------------------------------------------------' )
712prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True)
713
714if SRF :
715    echo ( '\n====================================================================================' )
716    echo ( f'-- SRF changes  -- {Title} ' )
717
718    if Routing == 'SIMPLE' :
719        RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   )
720        RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   )
721        RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] )
722        RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
723        RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
724        RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
725
726        RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   )
727        RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   )
728        RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] )
729
730        RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
731        RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
732        RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
733
734    if Routing == 'SECHIBA' :
735        RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
736        RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
737        RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
738        RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
739        RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
740        RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
741
742        RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
743        RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
744        RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
745        RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
746        RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
747        RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
748
749    RUN_mas_wat_beg = Sprec ( [RUN_mas_wat_fast_beg , RUN_mas_wat_slow_beg, RUN_mas_wat_stream_beg,
750                            RUN_mas_wat_flood_beg, RUN_mas_wat_lake_beg, RUN_mas_wat_pond_beg] )
751
752    RUN_mas_wat_end = Sprec ( [RUN_mas_wat_fast_end  , RUN_mas_wat_slow_end , RUN_mas_wat_stream_end,
753                            RUN_mas_wat_flood_end , RUN_mas_wat_lake_end , RUN_mas_wat_pond_end] )
754
755    dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg
756    dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg
757    dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg
758    dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg
759    dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg
760    dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg
761
762    dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg
763
764    echo ( f'\nRunoff reservoirs -- {Title} ')
765    echo (  '------------------------------------------------------------------------------------' )
766    echo ( f'RUN_mas_wat_fast_beg   = {RUN_mas_wat_fast_beg  :12.6e} kg | RUN_mas_wat_fast_end   = {RUN_mas_wat_fast_end  :12.6e} kg ' )
767    echo ( f'RUN_mas_wat_slow_beg   = {RUN_mas_wat_slow_beg  :12.6e} kg | RUN_mas_wat_slow_end   = {RUN_mas_wat_slow_end  :12.6e} kg ' )
768    echo ( f'RUN_mas_wat_stream_beg = {RUN_mas_wat_stream_beg:12.6e} kg | RUN_mas_wat_stream_end = {RUN_mas_wat_stream_end:12.6e} kg ' )
769    echo ( f'RUN_mas_wat_flood_beg  = {RUN_mas_wat_flood_beg :12.6e} kg | RUN_mas_wat_flood_end  = {RUN_mas_wat_flood_end :12.6e} kg ' )
770    echo ( f'RUN_mas_wat_lake_beg   = {RUN_mas_wat_lake_beg  :12.6e} kg | RUN_mas_wat_lake_end   = {RUN_mas_wat_lake_end  :12.6e} kg ' )
771    echo ( f'RUN_mas_wat_pond_beg   = {RUN_mas_wat_pond_beg  :12.6e} kg | RUN_mas_wat_pond_end   = {RUN_mas_wat_pond_end  :12.6e} kg ' )
772    echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' )
773
774    echo ( '------------------------------------------------------------------------------------' )
775
776    prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True )
777    prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True )
778    prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True )
779    prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True )
780    prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True )
781    prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True )
782    prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True )
783
784    echo ( f'\nWater content in routing  -- {Title} ' )
785    echo ( '------------------------------------------------------------------------------------' )
786    echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' )
787    prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   )
788
789    echo ( '\n====================================================================================' )
790    print ( 'Reading SRF restart')
791    SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.)
792    SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.)
793    SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.)
794    SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.)
795
796    SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.)
797    SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.)
798    SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.)
799    SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.)
800
801    if LMDZ :
802        SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1d='cell')
803        SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1d='cell')
804        SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1d='cell')
805        SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1d='cell')
806        SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1d='cell')
807        SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1d='cell')
808        SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1d='cell')
809        SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1d='cell')
810
811
812    # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood
813
814    SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg
815    SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end
816
817    echo ( '\n====================================================================================' )
818    print ('Computing integrals')
819
820    print ( ' 1/8', end='' ) ; sys.stdout.flush ()
821    SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  )
822    print ( ' 2/8', end='' ) ; sys.stdout.flush ()
823    SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg )
824    print ( ' 3/8', end='' ) ; sys.stdout.flush ()
825    SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        )
826    print ( ' 4/8', end='' ) ; sys.stdout.flush ()
827    SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     )
828    print ( ' 5/8', end='' ) ; sys.stdout.flush ()
829
830    SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  )
831    print ( ' 6/8', end='' ) ; sys.stdout.flush ()
832    SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end )
833    print ( ' 7/8', end='' ) ; sys.stdout.flush ()
834    SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        )
835    print ( ' 8/8', end='' ) ; sys.stdout.flush ()
836    SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     )
837
838    print (' -- ') ; sys.stdout.flush ()
839
840    dSRF_mas_watveg   = Sprec ( [SRF_mas_watveg_end , -SRF_mas_watveg_beg] )
841    dSRF_mas_watsoil  = Sprec ( [SRF_mas_watsoil_end, -SRF_mas_watsoil_beg] )
842    dSRF_mas_snow     = Sprec ( [SRF_mas_snow_end   , -SRF_mas_snow_beg] )
843    dSRF_mas_lake     = Sprec ( [SRF_mas_lake_end   , -SRF_mas_lake_beg] )
844
845    echo ( '------------------------------------------------------------------------------------' )
846    echo ( f'\nSurface reservoirs  -- {Title} ')
847    echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' )
848    echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' )
849    echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' )
850    echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' )
851
852    prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True )
853    prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True )
854    prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True )
855    prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True )
856
857    SRF_mas_wat_beg = Sprec ( [SRF_mas_watveg_beg , SRF_mas_watsoil_beg, SRF_mas_snow_beg] )
858    SRF_mas_wat_end = Sprec ( [SRF_mas_watveg_end , SRF_mas_watsoil_end, SRF_mas_snow_end] )
859
860    dSRF_mas_wat    = Sprec ( [+SRF_mas_watveg_end , +SRF_mas_watsoil_end, +SRF_mas_snow_end,
861                               -SRF_mas_watveg_beg , -SRF_mas_watsoil_beg, -SRF_mas_snow_beg]  )
862
863    echo ( '------------------------------------------------------------------------------------' )
864    echo ( f'Water content in surface  -- {Title} ' )
865    echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ')
866    prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True )
867
868    echo ( '------------------------------------------------------------------------------------' )
869    echo ( 'Water content in  ATM + SRF + RUN + LAKE' )
870    echo ( '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 ) )
873    prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True)
874
875echo ( '\n====================================================================================' )
876echo ( f'-- ATM Fluxes  -- {Title} ' )
877
878if ATM_HIS == 'latlon' :
879    echo ( ' latlon case' )
880    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1d='cell' )
881    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1d='cell' )
882    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1d='cell' )
883    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1d='cell' )
884    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' )
885    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1d='cell' )
886    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1d='cell' )
887    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1d='cell' )
888    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1d='cell' )
889    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1d='cell' )
890    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1d='cell' )
891    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1d='cell' )
892    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1d='cell' )
893    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1d='cell' )
894    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1d='cell' )
895    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1d='cell' )
896    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1d='cell' )
897    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1d='cell' )
898    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1d='cell' )
899    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1d='cell' )
900    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1d='cell' )
901    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1d='cell' )
902    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' )
903    echo ( 'End of LATLON case')
904
905if ATM_HIS == 'ico' :
906    echo (' ico case')
907    ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
908    ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
909    ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
910    ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
911    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
912    ATM_fqcalving   = rprec (d_ATM_his ['fqcalving'])
913    ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  )
914    ATM_precip      = rprec (d_ATM_his ['precip']   )
915    ATM_snowf       = rprec (d_ATM_his ['snow']     )
916    ATM_evap        = rprec (d_ATM_his ['evap']     )
917    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
918    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
919    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
920    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
921    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
922    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
923    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
924    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
925    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
926    ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
927    ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
928    ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
929    ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
930    ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
931    ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
932    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
933    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
934    echo ( 'End of ico case ')
935
936echo ( 'ATM wprecip_oce' )
937ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
938ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
939ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
940ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
941
942ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic
943ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic
944ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic
945ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic
946ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic
947ATM_wemp        = ATM_wevap - ATM_wprecip
948ATM_emp         = ATM_evap  - ATM_precip
949
950ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
951ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
952ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
953ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
954ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
955
956ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter
957ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce
958ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic
959ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic
960ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce
961
962if SRF :
963    if RUN_HIS == 'latlon' :
964        echo ( 'RUN costalflow Grille LATLON' )
965        if TestInterp :
966            echo ( 'RUN runoff TestInterp' )
967            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1d='cell' )
968            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1d='cell' )
969        else :
970            echo ( 'RUN runoff' )
971            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1d='cell' )
972            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1d='cell' )
973
974        RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1d='cell' )
975        RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1d='cell' )
976        RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1d='cell' )
977        RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1d='cell' )
978        RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1d='cell' )
979
980    if RUN_HIS == 'ico' :
981        echo ( 'RUN costalflow Grille ICO' )
982        RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])
983        RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )
984        RUN_runoff      =  rprec (d_RUN_his ['runoff']     )
985        RUN_drainage    =  rprec (d_RUN_his ['drainage']   )
986        RUN_riversret   =  rprec (d_RUN_his ['riversret']  )
987
988        RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
989        RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
990
991    Step = 0
992
993    if SRF_HIS == 'latlon' :
994        if TestInterp :
995            echo ( 'SRF rain TestInterp' )
996            SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1d='cell')
997            SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1d='cell')
998            SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1d='cell')
999            SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1d='cell')
1000            SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1d='cell' )
1001            #SRF_rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs )
1002            #SRF_evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs )
1003            #SRF_snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs )
1004            #SRF_subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs )
1005            #SRF_transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs )
1006        else :
1007            echo ( 'SRF rain' )
1008            SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1d='cell')
1009            SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1d='cell')
1010            SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1d='cell')
1011            SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1d='cell')
1012            SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1d='cell' )
1013
1014    if SRF_HIS == 'ico' :
1015        echo ( 'SRF rain')
1016        SRF_rain     = rprec (d_SRF_his ['rain'] )
1017        SRF_evap     = rprec (d_SRF_his ['evap'] )
1018        SRF_snowf    = rprec (d_SRF_his ['snowf'])
1019        SRF_subli    = rprec (d_SRF_his ['subli'])
1020        SRF_transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget')
1021
1022    echo ( 'SRF emp' )
1023    SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
1024    SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units']
1025
1026    ## Correcting units of SECHIBA variables
1027    def mmd2si ( pvar ) :
1028        '''Change unit from mm/d or m^3/s to kg/s if needed'''
1029        if 'units' in pvar.attrs :
1030            if pvar.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
1031                pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg/s'
1032            if pvar.attrs['units'] == 'mm/d' :
1033                pvar.values = pvar.values  * ATM_RHO * (1e-3/lmdz.RDAY) ;  pvar.attrs['units'] = 'kg/s'
1034            if pvar.attrs['units'] in ['m^3', 'm3', ] :
1035                pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg'
1036
1037    for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
1038        zvar = locals()['RUN_' + var]
1039        mmd2si (zvar)
1040
1041    for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
1042        zvar = locals()['SRF_' + var]
1043        mmd2si (zvar)
1044
1045    echo ( 'RUN input' )
1046    RUN_input  = RUN_runoff      + RUN_drainage
1047    RUN_output = RUN_coastalflow + RUN_riverflow
1048
1049echo ( 'ATM flw_wbilo' )
1050ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
1051ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      )
1052ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    )
1053ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      )
1054ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      )
1055ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       )
1056
1057ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
1058ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
1059ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
1060ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
1061ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
1062# Type d'integration a verifier
1063ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
1064ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
1065
1066LIC_flx_calving     = LIC_flux_int ( ATM_fqcalving  )
1067LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    )
1068
1069echo ( 'ATM flx precip' )
1070ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
1071ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
1072ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
1073ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
1074
1075LIC_flx_precip      = LIC_flux_int ( ATM_precip     )
1076LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      )
1077LIC_flx_evap        = LIC_flux_int ( ATM_evap       )
1078LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  )
1079
1080echo ( 'ATM flx_wrain_ter' )
1081ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter )
1082ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce )
1083ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic )
1084ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic )
1085ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea )
1086
1087ATM_flx_wsnow_ter    = ATM_flux_int ( ATM_wsnow_ter )
1088ATM_flx_wsnow_oce    = ATM_flux_int ( ATM_wsnow_oce )
1089ATM_flx_wsnow_lic    = ATM_flux_int ( ATM_wsnow_lic )
1090ATM_flx_wsnow_sic    = ATM_flux_int ( ATM_wsnow_sic )
1091ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea )
1092
1093echo ( 'ATM flx_evap_ter' )
1094ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter )
1095ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce )
1096ATM_flx_wevap_lic    = ATM_flux_int ( ATM_wevap_lic )
1097ATM_flx_wevap_sic    = ATM_flux_int ( ATM_wevap_sic )
1098ATM_flx_wevap_sea    = ATM_flux_int ( ATM_wevap_sea )
1099ATM_flx_wprecip_lic  = ATM_flux_int ( ATM_wprecip_lic )
1100ATM_flx_wprecip_oce  = ATM_flux_int ( ATM_wprecip_oce )
1101ATM_flx_wprecip_sic  = ATM_flux_int ( ATM_wprecip_sic )
1102ATM_flx_wprecip_ter  = ATM_flux_int ( ATM_wprecip_ter )
1103ATM_flx_wprecip_sea  = ATM_flux_int ( ATM_wprecip_sea )
1104ATM_flx_wemp_lic     = ATM_flux_int ( ATM_wemp_lic )
1105ATM_flx_wemp_oce     = ATM_flux_int ( ATM_wemp_oce )
1106ATM_flx_wemp_sic     = ATM_flux_int ( ATM_wemp_sic )
1107ATM_flx_wemp_ter     = ATM_flux_int ( ATM_wemp_ter )
1108ATM_flx_wemp_sea     = ATM_flux_int ( ATM_wemp_sea )
1109
1110ATM_flx_emp          = ATM_flux_int ( ATM_emp )
1111
1112if SRF :
1113    echo ( 'RUN flx_coastal' )
1114    RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
1115    echo ( 'RUN flx_river' )
1116    RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
1117    echo ( 'RUN flx_coastal_cpl' )
1118    RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
1119    echo ( 'RUN flx_river_cpl' )
1120    RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
1121    echo ( 'RUN flx_drainage' )
1122    RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
1123    echo ( 'RUN flx_riversset' )
1124    RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
1125    echo ( 'RUN flx_runoff' )
1126    RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
1127    echo ( 'RUN flx_input' )
1128    RUN_flx_input       = SRF_flux_int ( RUN_input      )
1129    echo ( 'RUN flx_output' )
1130    RUN_flx_output      = ONE_flux_int ( RUN_output     )
1131
1132    echo ( 'RUN flx_bil' ) ; Step += 1
1133    #RUN_flx_bil    = RUN_flx_input   - RUN_flx_output
1134    #RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river
1135
1136    RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output)
1137    RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow)
1138
1139prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )
1140prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )
1141prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )
1142prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )
1143prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )
1144prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True)
1145prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True)
1146prtFlux ('calving              ', ATM_flx_calving       , 'f' )
1147prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )
1148prtFlux ('precip               ', ATM_flx_precip        , 'f' )
1149prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )
1150prtFlux ('evap                 ', ATM_flx_evap          , 'f' )
1151prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' )
1152
1153prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' )
1154prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' )
1155prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' )
1156prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' )
1157prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True )
1158
1159echo ( 'Errors <field> vs. wbil_<field>' )
1160prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True )
1161prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True )
1162prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True )
1163prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True )
1164
1165if SRF :
1166    echo ( '\n====================================================================================' )
1167    echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
1168    prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' )
1169    prtFlux ('riverflow     ', RUN_flx_river      , 'f' )
1170    prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' )
1171    prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )
1172    prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )
1173    prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )
1174    prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )
1175    prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )
1176    prtFlux ('river in      ', RUN_flx_input      , 'f' )
1177    prtFlux ('river out     ', RUN_flx_output     , 'f' )
1178    prtFlux ('river bil     ', RUN_flx_bil        , 'f' )
1179
1180ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+# ATM_flx_fqfonte #+ RUN_flx_river
1181
1182
1183echo ('')
1184#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 ))
1185
1186#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 ))
1187
1188ATM_flx_toSRF = -ATM_flx_wbilo_ter
1189
1190echo (' ')
1191echo ( '\n====================================================================================' )
1192echo ( f'--  Atmosphere  -- {Title} ' )
1193echo ( f'Mass begin = {DYN_mas_wat_beg:12.6e} kg | Mass end = {DYN_mas_wat_end:12.6e} kg' )
1194prtFlux ( 'dmass (atm)  = ', dDYN_mas_wat , 'e', True )
1195prtFlux ( 'Sum wbilo_*  = ', ATM_flx_wbilo, 'e', True )
1196prtFlux ( 'E-P          = ', ATM_flx_emp  , 'e', True )
1197echo ( ' ' )
1198prtFlux ( 'Water loss atm from wbil_*', ATM_flx_wbilo - dDYN_mas_wat, 'f', True )
1199echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) )
1200
1201echo (' ')
1202prtFlux ( 'Water loss atm from E-P', ATM_flx_emp  - dDYN_mas_wat , 'f', True )
1203echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) )
1204echo (' ')
1205
1206ATM_error =  ATM_flx_emp  - dDYN_mas_wat
1207
1208
1209echo (' ')
1210echo ( '\n====================================================================================' )
1211
1212LIC_flx_budget1 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_calving , -LIC_flx_fqfonte] )
1213LIC_flx_budget2 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_calving , -LIC_flx_fqfonte] )
1214LIC_flx_budget3 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_runlic] )
1215LIC_flx_budget4 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_runlic] )
1216
1217echo ( f'--  LIC  -- {Title} ' )
1218echo ( f'Mass total   begin = {LIC_mas_wat_beg    :12.6e} kg | Mass end = {LIC_mas_wat_end    :12.6e} kg' )
1219echo ( f'Mass snow    begin = {LIC_mas_sno_beg    :12.6e} kg | Mass end = {LIC_mas_sno_end    :12.6e} kg' )
1220echo ( f'Mass qs      begin = {LIC_mas_qs_beg     :12.6e} kg | Mass end = {LIC_mas_qs_end     :12.6e} kg' )
1221echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' )
1222prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=45 )
1223prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=45 )
1224prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=45 )
1225prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=45 )
1226prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=45 )
1227prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=45 )
1228prtFlux ( 'LIC LIC_flx_fqfonte   ', LIC_flx_fqfonte       , 'f', True, width=45 )
1229prtFlux ( 'LIC LIC_flx_calving   ', LIC_flx_calving       , 'f', True, width=45 )
1230prtFlux ( 'LIC LIC_flx_runofflic ', LIC_flx_runlic        , 'f', True, width=45 )
1231prtFlux ( 'LIC fqfonte + calving ', LIC_flx_calving+LIC_flx_fqfonte , 'f', True, width=45 )
1232prtFlux ( 'LIC fluxes 1 ( wemp_lic   - fqcalving - fqfonte)) ', LIC_flx_budget1              , 'f', True, width=45 )
1233prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)   ', LIC_flx_budget2              , 'f', True, width=45 )
1234prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)    ', LIC_flx_budget3              , 'f', True, width=45 )
1235prtFlux ( 'LIC fluxes 3 ( wemp_lic  - runofflic*frac_lic)    ', LIC_flx_budget4              , 'f', True, width=45 )
1236prtFlux ( 'LIC error 1                                       ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=45 )
1237prtFlux ( 'LIC error 2                                       ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=45 )
1238prtFlux ( 'LIC error 3                                       ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=45 )
1239echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) )
1240echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) )
1241echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) )
1242
1243if SRF :
1244    echo ( '\n====================================================================================' )
1245    echo ( f'-- SECHIBA fluxes  -- {Title} ' )
1246
1247    SRF_flx_rain     = SRF_flux_int ( SRF_rain     )
1248    SRF_flx_evap     = SRF_flux_int ( SRF_evap     )
1249    SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    )
1250    SRF_flx_subli    = SRF_flux_int ( SRF_subli    )
1251    SRF_flx_transpir = SRF_flux_int ( SRF_transpir )
1252    SRF_flx_emp      = SRF_flux_int ( SRF_emp      )
1253
1254    RUN_flx_torouting   = SRF_flux_int  ( RUN_runoff       + RUN_drainage)
1255    RUN_flx_fromrouting = ONE_flux_int  ( RUN_coastalflow + RUN_riverflow )
1256
1257    SRF_flx_all =  SRF_flux_int  ( SRF_rain + SRF_snowf - SRF_evap - RUN_runoff - RUN_drainage )
1258
1259    prtFlux ('rain         ', SRF_flx_rain       , 'f' )
1260    prtFlux ('evap         ', SRF_flx_evap       , 'f' )
1261    prtFlux ('snowf        ', SRF_flx_snowf      , 'f' )
1262    prtFlux ('E-P          ', SRF_flx_emp        , 'f' )
1263    prtFlux ('subli        ', SRF_flx_subli      , 'f' )
1264    prtFlux ('transpir     ', SRF_flx_transpir   , 'f' )
1265    prtFlux ('to routing   ', RUN_flx_torouting  , 'f' )
1266    prtFlux ('budget       ', SRF_flx_all        , 'f', small=True )
1267
1268    echo ( '\n------------------------------------------------------------------------------------' )
1269    echo ( 'Water content in surface ' )
1270    echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' )
1271    prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True)
1272    prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True )
1273    echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) )
1274
1275    echo ( '\n====================================================================================' )
1276    echo ( f'-- Check ATM vs. SRF -- {Title} ' )
1277    prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' )
1278    prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' )
1279    prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' )
1280    prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True)
1281    echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) )
1282
1283    echo ('')
1284    echo ( '\n====================================================================================' )
1285    echo ( f'-- RUNOFF fluxes  -- {Title} ' )
1286    RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal
1287    prtFlux ('runoff    ', RUN_flx_runoff      , 'f' )
1288    prtFlux ('drainage  ', RUN_flx_drainage    , 'f' )
1289    prtFlux ('run+drain ', RUN_flx_torouting   , 'f' )
1290    prtFlux ('river     ', RUN_flx_river       , 'f' )
1291    prtFlux ('coastal   ', RUN_flx_coastal     , 'f' )
1292    prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' )
1293    prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True)
1294
1295    echo ( '\n------------------------------------------------------------------------------------' )
1296    echo ( f'Water content in routing+lake  -- {Title} ' )
1297    echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1298    prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True)
1299    prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True )
1300    echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) )
1301
1302    echo ( '\n------------------------------------------------------------------------------------' )
1303    echo ( f'Water content in routing  -- {Title} ' )
1304    echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1305    prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  )
1306    prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True)
1307    echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) )
1308
1309echo ( ' ' )
1310echo ( f'{Title = }' )
1311
1312echo ( 'SVN Information' )
1313for kk in SVN.keys():
1314    print ( SVN[kk] )
1315
1316## Write the full configuration
1317params_out = open (FullIniFile, 'w', encoding = 'utf-8')
1318params = wu.dict2config ( dpar )
1319params.write ( params_out )
1320params_out.close ()
Note: See TracBrowser for help on using the repository browser.