source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6688

Last change on this file since 6688 was 6688, checked in by omamce, 7 months ago

O.M. :

WATER_BUDGET : put function common to ATM and OCE in WAterUtils library

use of dictionary for all input parameters
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 61.4 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 types
28import configparser
29
30## Import needed scientific modules
31import numpy as np
32import xarray as xr
33
34## Import local modules
35import WaterUtils as wu
36import lmdz
37
38## Read command line arguments
39## ---------------------------
40print ( "Name of Python script:", sys.argv[0] )
41IniFile = sys.argv[1]
42
43# Test existence of IniFile
44if not os.path.exists (IniFile ) :
45    raise FileExistsError ( f"File not found : {IniFile = }" )
46
47if 'full' in IniFile or 'ATM' in IniFile :
48    FullIniFile = IniFile
49else :
50    FullIniFile = 'ATM_' + IniFile
51
52print ("Output file : ", FullIniFile )
53
54## Experiment parameters : read in .ini file
55## -----------------------------------------
56dpar = wu.ReadConfig ( IniFile )
57
58## Configure all needed parameter from existant parameters
59## -------------------------------------------------------
60dpar = wu.SetDatesAndFiles ( dpar )
61
62## Output file with water budget diagnostics
63## -----------------------------------------
64f_out = dpar['Files']['f_out']
65
66## Put dpar values in local namespace
67## ----------------------------------
68for Section in dpar.keys () : 
69    print ( f'\nReading [{Section}]' )
70    for VarName in dpar[Section].keys() :
71        globals()[VarName] = dpar[Section][VarName]
72        print ( f'    {VarName:21} set to : {globals()[VarName]}' )
73
74## Debuging and timer
75Timer = wu.functools.partial (wu.Timer, debug=Debug, timer=Timing)
76
77## Useful functions
78## ----------------
79if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float :
80    def rprec (ptab) :
81        '''This version does nothing
82
83        rprec may be used to reduce floating precision when reading history files
84        '''
85        return ptab
86else                 :
87    def rprec (ptab) :
88        '''Returns float with a different precision'''
89        return ptab.astype(readPrec).astype(float)
90
91def kg2Sv    (pval, rho=ATM_RHO) :
92    '''From kg to Sverdrup'''
93    return pval/dtime_sec*1.0e-6/rho
94
95def kg2myear (pval, rho=ATM_RHO) :
96    '''From kg to m/year'''
97    return pval/ATM.aire_sea_tot/rho/NbYear
98
99def var2prt (pvar, small=False, rho=ATM_RHO) :
100    '''Formats value for printing'''
101    if small :  return pvar, kg2Sv (pvar, rho=rho)*1000., kg2myear (pvar, rho=rho)*1000
102    else     :  return pvar, kg2Sv (pvar, rho=rho)      , kg2myear (pvar, rho=rho)
103
104def prtFlux (Desc, pvar, Form='F', small=False, rho=ATM_RHO, width=15) :
105    '''Pretty print of formatted value'''
106    if small :
107        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
108        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
109    else :
110        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
111        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
112    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt (pvar, small=small, rho=rho), width=width ) )
113
114def echo (string, end='\n') :
115    '''Function to print to stdout *and* output file'''
116    print ( str(string), end=end  )
117    sys.stdout.flush ()
118    f_out.write ( str(string) + end )
119    f_out.flush ()
120
121## Open history files
122## ------------------
123d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
124if SECHIBA :
125    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
126    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
127    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
128
129echo ( f'{file_ATM_his = }' )
130if SECHIBA :
131    echo ( f'{file_SRF_his = }' )
132   
133## Compute run length
134## ------------------
135dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
136echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' )
137dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
138dpar['Experiment']['dtime_sec'] = dtime_sec
139
140##-- Compute length of each period
141dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
142echo ( '\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
143dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
144dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
145dtime_per_sec.attrs['unit'] = 's'
146dpar['Experiment']['dtime_per_sec'] = dtime_per_sec
147
148##-- Number of years (approximative)
149NbYear = dtime_sec / YEAR_LENGTH
150dpar['Experiment']['NbYear'] = NbYear
151
152## Define restart periods and file names
153## -------------------------------------
154dpar = wu.SetRestartNames ( dpar, f_out ) 
155
156## Put dpar values in local namespace
157## ----------------------------------
158for Section in dpar.keys () : 
159    print ( f'\nReading [{Section}]' )
160    for VarName in dpar[Section].keys() :
161        locals()[VarName] = dpar[Section][VarName]
162        print ( f'    {VarName:21} set to : {locals()[VarName]}' )
163       
164## Extract restart files from tar
165## ------------------------------
166liste_beg = [file_ATM_beg, ]
167liste_end = [file_ATM_end, ]
168if file_DYN_beg : liste_beg.append ( file_DYN_beg )
169if file_DYN_end : liste_end.append ( file_DYN_end )
170if file_SRF_beg : liste_beg.append ( file_SRF_beg )
171if file_SRF_end : liste_end.append ( file_SRF_end )
172
173if ICO :
174    if not file_DYN_aire :
175        file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ResolAtm+'_grid.nc' )
176    dpar['Files']['file_DYN_aire'] = file_DYN_aire
177
178if SECHIBA and Routing == 'SIMPLE' :
179    liste_beg.append ( file_RUN_beg )
180    liste_end.append ( file_RUN_end )
181
182echo ( '\nExtract restart files from tar : ATM, ICO', end='')
183if SECHIBA : echo ( ' and SECHIBA')
184else   : echo (' ')
185
186## Write the full configuration
187## ----------------------------
188params_out = open (FullIniFile, 'w', encoding = 'utf-8')
189params = wu.dict2config ( dpar )
190params.write ( params_out )
191params_out.close ()
192   
193@Timer
194def extract ( file_name=file_ATM_beg, tar_restart=tar_restart_end, file_dir_comp=FileDir ) :
195    '''
196    Extract restart files from a tar
197    '''
198    echo ( f'----------')
199    error_count = 0
200    echo ( f'file to extract : {file_name = }' )
201    if os.path.exists ( os.path.join (FileDir, file_name) ) :
202        echo ( f'file found      : {file_name = }' )
203    else :
204        echo ( f'file not found  : {file_name = }' )
205        base_resFile = os.path.basename (file_name)
206        if os.path.exists ( tar_restart ) :
207            command =  f'cd {file_dir_comp} ; tar xf {tar_restart} {base_resFile}'
208            echo ( f'{command = }' )
209            err = os.system ( command )
210            if err != 0 :
211                if ContinueOnError :
212                    error_count += 1
213                    echo ( f'****** Command failed : {command}' )
214                    echo ( '****** Trying to continue' )
215                    echo ( ' ')
216                else :
217                    raise OSError ( f'**** command failed : {command} - Stopping' )
218            else :
219                echo ( f'tar done : {base_resFile}' )
220        else :
221            echo ( f'****** Tar restart file {tar_restart = } not found ' )
222            if ContinueOnError :
223                error_count += 1
224                echo ( f'****** Command failed : {command}' )
225                echo (  '****** Trying to continue' )
226                echo (  ' ')
227            else :
228                raise OSError ( f'****** tar file not found {tar_restart = } - Stopping' )
229    return error_count
230
231ErrorCount = 0
232
233echo ( f'Extract {file_ATM_beg = }' )
234ErrorCount += extract ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, file_dir_comp=FileDir )
235echo ( f'Extract {file_DYN_beg = }' )
236ErrorCount += extract ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, file_dir_comp=FileDir )
237
238echo ( f'Extract {file_ATM_end = }' )
239ErrorCount += extract ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, file_dir_comp=FileDir )
240echo ( f'Extract {file_DYN_end = }' )
241ErrorCount += extract ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, file_dir_comp=FileDir )
242
243if SECHIBA :
244    echo ( f'Extract {file_SRF_beg = }' )
245    ErrorCount += extract ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, file_dir_comp=FileDir )
246    echo ( f'Extract {file_SRF_end = }' )
247    ErrorCount += extract ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, file_dir_comp=FileDir )
248
249    if Routing == 'SIMPLE' :
250        echo ( f'Extract {file_RUN_beg = }' )
251        ErrorCount += extract ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, file_dir_comp=FileDir )
252        echo ( f'Extract {file_RUN_end = }' )
253        ErrorCount += extract ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, file_dir_comp=FileDir )
254
255##-- Exit in case of error in the extracting file phase
256if ErrorCount > 0 :
257    echo ( '  ' )
258    raise RuntimeError ( f'**** Some files missing - Stopping - {ErrorCount = }' )
259
260##
261echo ('\nOpening ATM SECHIBA and ICO restart files')
262d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
263d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
264if SECHIBA :
265    d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
266    d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
267d_DYN_beg = xr.open_dataset ( os.path.join (FileDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
268d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
269
270if SECHIBA :
271    for var in d_SRF_beg.variables :
272        d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.)
273        d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.)
274
275    if Routing == 'SIMPLE' :
276        d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
277        d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
278
279@Timer
280def to_cell ( dd, newname='cell' ) :
281    '''Set space dimension to newname
282    '''
283    for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] :
284        if oldname in dd.dims and oldname != newname :
285            dd = dd.rename ( {oldname : newname} )
286    return dd
287
288d_ATM_beg = to_cell ( d_ATM_beg )
289d_ATM_end = to_cell ( d_ATM_end )
290if SECHIBA :
291    d_SRF_beg = to_cell ( d_SRF_beg )
292    d_SRF_end = to_cell ( d_SRF_end )
293d_DYN_beg = to_cell ( d_DYN_beg )
294d_DYN_end = to_cell ( d_DYN_end )
295
296if SECHIBA and Routing == 'SIMPLE' :
297    d_RUN_beg = to_cell ( d_RUN_beg )
298    d_RUN_end = to_cell ( d_RUN_end )
299
300d_ATM_his = to_cell ( d_ATM_his )
301if SECHIBA : d_SRF_his = to_cell ( d_SRF_his )
302
303echo ( f'{file_ATM_beg = }' )
304echo ( f'{file_ATM_end = }' )
305echo ( f'{file_DYN_beg = }' )
306echo ( f'{file_DYN_end = }' )
307if SECHIBA :
308    echo ( f'{file_SRF_beg = }' )
309    echo ( f'{file_SRF_end = }' )
310    if Routing == 'SIMPLE' :
311        echo ( f'{file_RUN_beg = }' )
312        echo ( f'{file_RUN_end = }' )
313
314## Compute aire, fractions, etc ...
315## --------------------------------
316if ICO :
317    # Area on icosahedron grid
318    #echo ( f'{file_DYN_aire = }' )
319    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
320else :
321    d_DYN_aire = None
322       
323dpar, ATM = wu.ComputeGridATM ( dpar, d_ATM_his ) 
324dpar, SRF = wu.ComputeGridSRF ( dpar, d_SRF_his )
325dpar, DYN = wu.ComputeGridDYN ( dpar, ATM, d_DYN_aire, d_ATM_beg )
326
327if ICO and ATM_HIS == 'ico' :
328    # Comparaison des aires sur ATM et DYN
329    aire_diff = ATM.aire - DYN.aire
330    echo ( f'Difference Aire hist file vs. grid file mean:{aire_diff.mean().values} min:{aire_diff.min().values} max:{aire_diff.max().values} ' )
331
332# Functions computing integrals and sum
333@Timer
334def DYN_stock_int (stock) :
335    '''Integrate (* surface) stock on atmosphere grid'''
336    return wu.P1sum ( stock * DYN.aire )
337
338@Timer
339def ATM_flux_int (flux) :
340    '''Integrate (* time * surface) flux on atmosphere grid'''
341    return wu.P1sum ( flux * dtime_per_sec * ATM.aire )
342
343@Timer
344def LIC_flux_int (flux) :
345    '''Integrate (* time * surface) flux on land ice grid'''
346    return wu.P1sum ( flux * dtime_per_sec * ATM.aire_flic )
347
348if SECHIBA :
349    @Timer
350    def SRF_stock_int (stock) :
351        '''Integrate (* surface) stock on land grid'''
352        return wu.P1sum ( stock * DYN.aire_fter )
353
354    @Timer
355    def SRF_flux_int (flux) :
356        '''Integrate (* time * surface) flux on land grid'''
357        return wu.P1sum ( flux * dtime_per_sec * SRF.aire )
358
359@Timer
360def ONE_stock_int (stock) :
361    '''Sum stock'''
362    return wu.P1sum ( stock ) 
363
364@Timer
365def ONE_flux_int (flux) :
366    '''Integrate (* time) flux on area=1 grid'''
367    return wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() )
368
369@Timer
370def Sprec ( tlist ) :
371    '''Accurate sum of list of scalar elements'''
372    return wu.Psum ( np.array ( tlist) )
373
374echo ('')
375echo ( f'ATM DYN : Area of atmosphere        : {DYN.aire_tot:18.9e}' )
376echo ( f'ATM HIS : Area of atmosphere        : {ATM.aire_tot:18.9e}' )
377echo ( f'ATM HIS : Area of ter in atmosphere : {ATM.aire_ter_tot:18.9e}' )
378echo ( f'ATM HIS : Area of lic in atmosphere : {ATM.aire_lic_tot:18.9e}' )
379if SECHIBA :
380    echo ( f'ATM SECHIBA : Area of atmosphere        : {SRF.aire_tot:18.9e}' )
381echo ('')
382echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN.aire_tot/(RA*RA*4*np.pi) ) )
383echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM.aire_tot/(RA*RA*4*np.pi) ) )
384echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM.aire_ter_tot/(RA*RA*4*np.pi) ) )
385if SECHIBA :
386    echo ( 'ATM SECHIBA : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF.aire_tot/(RA*RA*4*np.pi) ) )
387    echo ('')
388    echo ( f'ATM SECHIBA : Area of atmosphere (no contfrac): {ONE_stock_int (SRF.areas):18.9e}' )
389
390
391if np.abs (ATM.aire_tot/(RA*RA*4*np.pi) - 1.0) > 0.01 :
392    raise RuntimeError ('Error of atmosphere surface interpolated on lon/lat grid')
393
394echo ( '\n====================================================================================' )
395echo ( f'-- ATM changes in stores  -- {Title} ' )
396
397#-- Change in precipitable water from the atmosphere daily and monthly files
398#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
399
400echo ( 'ATM vertical grid' )
401ATM.Ahyb = d_ATM_his['Ahyb'].squeeze()
402ATM.Bhyb = d_ATM_his['Bhyb'].squeeze()
403
404echo ( 'Surface pressure' )
405if ICO :
406    DYN.psol_beg = d_DYN_beg['ps']
407    DYN.psol_end = d_DYN_end['ps']
408if LMDZ :
409    DYN.psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' )
410    DYN.psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' )
411
412echo ( '3D Pressure at the interface layers (not at scalar points)' )
413DYN.pres_beg = ATM.Ahyb + ATM.Bhyb * DYN.psol_beg
414DYN.pres_end = ATM.Ahyb + ATM.Bhyb * DYN.psol_end
415
416echo ( 'Check computation of pressure levels' )
417
418ind = np.empty (8)
419
420ind[0] = (DYN.pres_beg[ 0]-DYN.psol_beg).min().item()
421ind[1] = (DYN.pres_beg[ 0]-DYN.psol_beg).max().item()
422ind[2] = (DYN.pres_beg[-1]).min().item()
423ind[3] = (DYN.pres_beg[-1]).max().item()
424ind[4] = (DYN.pres_end[ 0]-DYN.psol_end).min().item()
425ind[5] = (DYN.pres_end[ 0]-DYN.psol_end).max().item()
426ind[6] = (DYN.pres_end[-1]).min().item()
427ind[7] = (DYN.pres_end[-1]).max().item()
428
429if any ( ind != 0) :
430    echo ( 'All values should be zero' )
431    echo ( f'(DYN.pres_beg[ 0]-DYN.psol_beg).min().item() = {ind[0]}' )
432    echo ( f'(DYN.pres_beg[ 0]-DYN.psol_beg).max().item() = {ind[1]}' )
433    echo ( f'(DYN.pres_beg[-1]).min().item()              = {ind[2]}' )
434    echo ( f'(DYN.pres_beg[-1]).max().item()              = {ind[3]}' )
435    echo ( f'(DYN.pres_end[ 0]-DYN.psol_end).min().item() = {ind[4]}' )
436    echo ( f'(DYN.pres_end[ 0]-DYN.psol_end).max().item() = {ind[5]}' )
437    echo ( f'(DYN.pres_end[-1]).min().item()              = {ind[6]}' )
438    echo ( f'(DYN.pres_end[-1]).max().item()              = {ind[7]}' )
439    raise RuntimeError
440
441klevp1 = ATM.Bhyb.shape[-1]
442cell   = DYN.psol_beg.shape[-1]
443klev   = klevp1 - 1
444
445echo ( 'Layer thickness (in pressure)' )
446DYN.mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
447DYN.mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
448
449echo ( 'Layer mass (kg/m2)' )
450for k in np.arange (klev) :
451    DYN.mass_beg[k,:] = ( DYN.pres_beg[k,:] - DYN.pres_beg[k+1,:] ) / GRAV
452    DYN.mass_end[k,:] = ( DYN.pres_end[k,:] - DYN.pres_end[k+1,:] ) / GRAV
453
454echo ( 'Column mass (kg/m2)' )
455DYN.mass_beg_2D = DYN.mass_beg.sum (dim='sigs')
456DYN.mass_end_2D = DYN.mass_end.sum (dim='sigs')
457
458echo ( 'Atmosphere mass (kg)' )
459DYN.mass_air_beg = DYN_stock_int ( DYN.mass_beg_2D )
460DYN.mass_air_end = DYN_stock_int ( DYN.mass_end_2D )
461
462echo ( 'Atmosphere mass change (kg)' )
463DYN.diff_mass_air = DYN.mass_air_end - DYN.mass_air_beg
464
465echo ( f'\nChange of atmosphere mass (dynamics)  -- {Title} ' )
466echo ( '------------------------------------------------------------------------------------' )
467echo ( 'DYN.mass_air_beg = {:12.6e} kg | DYN.mass_air_end = {:12.6e} kg'.format (DYN.mass_air_beg, DYN.mass_air_end) )
468echo ( f'dMass (atm) : {DYN.diff_mass_air:9.4e} kg' )
469echo ( f'dMass (atm) : {(DYN.diff_mass_air/DYN.mass_air_beg):9.4e} (-)' )
470
471echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' )
472if LMDZ :
473    if 'H2Ov' in d_DYN_beg.variables :
474        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' )
475        DYN.wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']   + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' )
476        DYN.wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']   + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' )
477    if 'H2O_g' in d_DYN_beg.variables :
478        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' )
479        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' )
480        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' )
481if ICO :
482    if   'H2Ov_g' in d_DYN_beg.variables :
483        echo ('reading ICO : H2Ov_g, H2Ov_l, H2Ov_s' )
484        DYN.wat_beg = d_DYN_beg['H2Ov_g'] + d_DYN_beg['H2Ov_l'] + d_DYN_beg['H2Ov_s']
485        DYN.wat_end = d_DYN_end['H2Ov_g'] + d_DYN_end['H2Ov_l'] + d_DYN_end['H2Ov_s']
486    if 'H2O_g' in d_DYN_beg.variables :
487        echo ('reading ICO : H2O_g, H2O_l, H2O_s' )
488        DYN.wat_beg = d_DYN_beg['H2O_g']  + d_DYN_beg['H2O_l']  + d_DYN_beg['H2O_s']
489        DYN.wat_end = d_DYN_end['H2O_g']  + d_DYN_end['H2O_l']  + d_DYN_end['H2O_s']
490    if 'q' in d_DYN_beg.variables :
491        echo ('reading ICO : q' )
492        DYN.wat_beg = d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2)
493        DYN.wat_end = d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2)
494
495if 'lev' in DYN.wat_beg.dims :
496    DYN.wat_beg = DYN.wat_beg.rename ( {'lev':'sigs'} )
497    DYN.wat_end = DYN.wat_end.rename ( {'lev':'sigs'} )
498
499echo ( 'Compute water content : vertical and horizontal integral' )
500
501DYN.wat_beg_2D =  (DYN.mass_beg * DYN.wat_beg).sum (dim='sigs')
502DYN.wat_end_2D =  (DYN.mass_end * DYN.wat_end).sum (dim='sigs')
503
504DYN.mass_wat_beg = DYN_stock_int ( DYN.wat_beg_2D )
505DYN.mass_wat_end = DYN_stock_int ( DYN.wat_end_2D )
506
507echo ( 'Variation of water content' )
508DYN.diff_mass_wat = DYN.mass_wat_end - DYN.mass_wat_beg
509
510echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' )
511echo ( '------------------------------------------------------------------------------------' )
512echo ( 'DYN.mass_wat_beg = {:12.6e} kg | DYN.mass_wat_end = {:12.6e} kg'.format (DYN.mass_wat_beg, DYN.mass_wat_end) )
513prtFlux ( 'dMass water (atm)', DYN.diff_mass_wat, 'e', True )
514echo (   f'dMass water (atm) : {(DYN.diff_mass_wat/DYN.mass_wat_beg):9.4e} (-)' )
515
516
517ATM.sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \
518              d_ATM_beg['SNOW03']*d_ATM_beg['FOCE'] + d_ATM_beg['SNOW04']*d_ATM_beg['FSIC']
519ATM.sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \
520              d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC']
521
522ATM.qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \
523              d_ATM_beg['QS03']  *d_ATM_beg['FOCE'] + d_ATM_beg['QS04']  *d_ATM_beg['FSIC']
524ATM.qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \
525              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC']
526
527ATM.qsol_beg     = d_ATM_beg['QSOL']
528ATM.qsol_end     = d_ATM_end['QSOL']
529
530ATM.LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']
531ATM.LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC']
532
533ATM.LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0']
534ATM.LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0']
535
536ATM.qs01_beg     = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
537ATM.qs02_beg     = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
538ATM.qs03_beg     = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
539ATM.qs04_beg     = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
540
541ATM.qs01_end     = d_ATM_end['QS01'] * d_ATM_end['FTER']
542ATM.qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC']
543ATM.qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE']
544ATM.qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC']
545
546ATM.LIC_qs_beg = ATM.qs02_beg
547ATM.LIC_qs_end = ATM.qs02_end
548
549ATM.mass_sno_beg     = DYN_stock_int ( ATM.sno_beg  )
550ATM.mass_sno_end     = DYN_stock_int ( ATM.sno_end  )
551
552ATM.mass_qs_beg      = DYN_stock_int ( ATM.qs_beg   )
553ATM.mass_qs_end      = DYN_stock_int ( ATM.qs_end   )
554ATM.mass_qsol_beg    = DYN_stock_int ( ATM.qsol_beg )
555ATM.mass_qs01_beg    = DYN_stock_int ( ATM.qs01_beg )
556ATM.mass_qs02_beg    = DYN_stock_int ( ATM.qs02_beg )
557ATM.mass_qs03_beg    = DYN_stock_int ( ATM.qs03_beg )
558ATM.mass_qs04_beg    = DYN_stock_int ( ATM.qs04_beg )
559ATM.mass_qsol_end    = DYN_stock_int ( ATM.qsol_end )
560ATM.mass_qs01_end    = DYN_stock_int ( ATM.qs01_end )
561ATM.mass_qs02_end    = DYN_stock_int ( ATM.qs02_end )
562ATM.mass_qs03_end    = DYN_stock_int ( ATM.qs03_end )
563ATM.mass_qs04_end    = DYN_stock_int ( ATM.qs04_end )
564
565ATM.LIC_mass_sno_beg     = DYN_stock_int ( ATM.LIC_sno_beg     )
566ATM.LIC_mass_sno_end     = DYN_stock_int ( ATM.LIC_sno_end     )
567ATM.LIC_mass_runlic0_beg = DYN_stock_int ( ATM.LIC_runlic0_beg )
568ATM.LIC_mass_runlic0_end = DYN_stock_int ( ATM.LIC_runlic0_end )
569
570ATM.LIC_mass_qs_beg  = ATM.mass_qs02_beg
571ATM.LIC_mass_qs_end  = ATM.mass_qs02_end
572ATM.LIC_mass_wat_beg = ATM.LIC_mass_qs_beg + ATM.LIC_mass_sno_beg
573ATM.LIC_mass_wat_end = ATM.LIC_mass_qs_end + ATM.LIC_mass_sno_end
574
575ATM.diff_mass_sno  = ATM.mass_sno_end  - ATM.mass_sno_beg
576ATM.diff_mass_qs   = ATM.mass_qs_end   - ATM.mass_qs_beg
577ATM.diff_mass_qsol = ATM.mass_qsol_end - ATM.mass_qsol_beg
578
579ATM.diff_mass_qs01 = ATM.mass_qs01_end - ATM.mass_qs01_beg
580ATM.diff_mass_qs02 = ATM.mass_qs02_end - ATM.mass_qs02_beg
581ATM.diff_mass_qs03 = ATM.mass_qs03_end - ATM.mass_qs03_beg
582ATM.diff_mass_qs04 = ATM.mass_qs04_end - ATM.mass_qs04_beg
583
584ATM.LIC_diff_mass_qs      = ATM.LIC_mass_qs_end      - ATM.LIC_mass_qs_beg
585ATM.LIC_diff_mass_sno     = ATM.LIC_mass_sno_end     - ATM.LIC_mass_sno_beg
586ATM.LIC_diff_mass_runlic0 = ATM.LIC_mass_runlic0_end - ATM.LIC_mass_runlic0_beg
587
588ATM.LIC_diff_mass_wat     = ATM.LIC_diff_mass_qs + ATM.LIC_diff_mass_sno # + ATM.LIC_diff_mass_runlic0
589
590echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' )
591echo ( '------------------------------------------------------------------------------------' )
592echo ( 'ATM.mass_sno_beg  = {:12.6e} kg | ATM.mass_sno_end  = {:12.6e} kg'.format (ATM.mass_sno_beg, ATM.mass_sno_end) )
593prtFlux ( 'dMass (neige atm) ', ATM.diff_mass_sno  , 'e', True  )
594
595echo ( f'\nChange of soil humidity content -- {Title} ' )
596echo ( '------------------------------------------------------------------------------------' )
597echo ( 'ATM.mass_qs_beg  = {:12.6e} kg | ATM.mass_qs_end  = {:12.6e} kg'.format (ATM.mass_qs_beg, ATM.mass_qs_end) )
598prtFlux ( 'dMass (neige atm) ', ATM.diff_mass_qs, 'e', True )
599
600echo ( f'\nChange of atmosphere water+snow content -- {Title}  ' )
601echo ( '------------------------------------------------------------------------------------' )
602prtFlux ( 'dMass (eau + neige atm) ', DYN.diff_mass_wat + ATM.diff_mass_sno , 'e', True)
603
604if SECHIBA :
605    echo ( '\n====================================================================================' )
606    echo ( f'-- SECHIBA changes  -- {Title} ' )
607
608    if Routing == 'SIMPLE' :
609        SRF.RUN_mass_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   )
610        SRF.RUN_mass_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   )
611        SRF.RUN_mass_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] )
612        SRF.RUN_mass_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
613        SRF.RUN_mass_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
614        SRF.RUN_mass_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
615
616        SRF.RUN_mass_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   )
617        SRF.RUN_mass_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   )
618        SRF.RUN_mass_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] )
619
620        SRF.RUN_mass_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
621        SRF.RUN_mass_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
622        SRF.RUN_mass_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
623
624    if Routing == 'SECHIBA' :
625        SRF.RUN_mass_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
626        SRF.RUN_mass_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
627        SRF.RUN_mass_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
628        SRF.RUN_mass_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
629        SRF.RUN_mass_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
630        SRF.RUN_mass_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
631
632        SRF.RUN_mass_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
633        SRF.RUN_mass_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
634        SRF.RUN_mass_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
635        SRF.RUN_mass_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
636        SRF.RUN_mass_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
637        SRF.RUN_mass_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
638
639    SRF.RUN_mass_wat_beg = Sprec ( [SRF.RUN_mass_wat_fast_beg , SRF.RUN_mass_wat_slow_beg, SRF.RUN_mass_wat_stream_beg,
640                                   SRF.RUN_mass_wat_flood_beg, SRF.RUN_mass_wat_lake_beg, SRF.RUN_mass_wat_pond_beg] )
641
642    SRF.RUN_mass_wat_end = Sprec ( [SRF.RUN_mass_wat_fast_end  , SRF.RUN_mass_wat_slow_end , SRF.RUN_mass_wat_stream_end,
643                                   SRF.RUN_mass_wat_flood_end , SRF.RUN_mass_wat_lake_end , SRF.RUN_mass_wat_pond_end] )
644
645    SRF.RUN_diff_mass_wat_fast   = SRF.RUN_mass_wat_fast_end   - SRF.RUN_mass_wat_fast_beg
646    SRF.RUN_diff_mass_wat_slow   = SRF.RUN_mass_wat_slow_end   - SRF.RUN_mass_wat_slow_beg
647    SRF.RUN_diff_mass_wat_stream = SRF.RUN_mass_wat_stream_end - SRF.RUN_mass_wat_stream_beg
648    SRF.RUN_diff_mass_wat_flood  = SRF.RUN_mass_wat_flood_end  - SRF.RUN_mass_wat_flood_beg
649    SRF.RUN_diff_mass_wat_lake   = SRF.RUN_mass_wat_lake_end   - SRF.RUN_mass_wat_lake_beg
650    SRF.RUN_diff_mass_wat_pond   = SRF.RUN_mass_wat_pond_end   - SRF.RUN_mass_wat_pond_beg
651
652    SRF.RUN_diff_mass_wat        = SRF.RUN_mass_wat_end        - SRF.RUN_mass_wat_beg
653
654    echo ( f'\nRunoff reservoirs -- {Title} ')
655    echo (  '------------------------------------------------------------------------------------' )
656    echo ( f'SRF.RUN_mass_wat_fast_beg   = {SRF.RUN_mass_wat_fast_beg  :12.6e} kg | SRF.RUN_mass_wat_fast_end   = {SRF.RUN_mass_wat_fast_end  :12.6e} kg ' )
657    echo ( f'SRF.RUN_mass_wat_slow_beg   = {SRF.RUN_mass_wat_slow_beg  :12.6e} kg | SRF.RUN_mass_wat_slow_end   = {SRF.RUN_mass_wat_slow_end  :12.6e} kg ' )
658    echo ( f'SRF.RUN_mass_wat_stream_beg = {SRF.RUN_mass_wat_stream_beg:12.6e} kg | SRF.RUN_mass_wat_stream_end = {SRF.RUN_mass_wat_stream_end:12.6e} kg ' )
659    echo ( f'SRF.RUN_mass_wat_flood_beg  = {SRF.RUN_mass_wat_flood_beg :12.6e} kg | SRF.RUN_mass_wat_flood_end  = {SRF.RUN_mass_wat_flood_end :12.6e} kg ' )
660    echo ( f'SRF.RUN_mass_wat_lake_beg   = {SRF.RUN_mass_wat_lake_beg  :12.6e} kg | SRF.RUN_mass_wat_lake_end   = {SRF.RUN_mass_wat_lake_end  :12.6e} kg ' )
661    echo ( f'SRF.RUN_mass_wat_pond_beg   = {SRF.RUN_mass_wat_pond_beg  :12.6e} kg | SRF.RUN_mass_wat_pond_end   = {SRF.RUN_mass_wat_pond_end  :12.6e} kg ' )
662    echo ( f'SRF.RUN_mass_wat_beg        = {SRF.RUN_mass_wat_beg       :12.6e} kg | SRF.RUN_mass_wat_end        = {SRF.RUN_mass_wat_end       :12.6e} kg ' )
663
664    echo ( '------------------------------------------------------------------------------------' )
665
666    prtFlux ( 'dMass (fast)   ', SRF.RUN_diff_mass_wat_fast  , 'e', True )
667    prtFlux ( 'dMass (slow)   ', SRF.RUN_diff_mass_wat_slow  , 'e', True )
668    prtFlux ( 'dMass (stream) ', SRF.RUN_diff_mass_wat_stream, 'e', True )
669    prtFlux ( 'dMass (flood)  ', SRF.RUN_diff_mass_wat_flood , 'e', True )
670    prtFlux ( 'dMass (lake)   ', SRF.RUN_diff_mass_wat_lake  , 'e', True )
671    prtFlux ( 'dMass (pond)   ', SRF.RUN_diff_mass_wat_pond  , 'e', True )
672    prtFlux ( 'dMass (all)    ', SRF.RUN_diff_mass_wat       , 'e', True )
673
674    echo ( f'\nWater content in routing  -- {Title} ' )
675    echo ( '------------------------------------------------------------------------------------' )
676    echo ( f'SRF.RUN_mass_wat_beg = {SRF.RUN_mass_wat_end:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg' )
677    prtFlux ( 'dMass (routing) ', SRF.RUN_diff_mass_wat , 'e', True   )
678
679    echo ( '\n====================================================================================' )
680    print ( 'Reading SECHIBA restart')
681    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.)
682    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.)
683    SRF.snow_beg        = d_SRF_beg['snow_beg']        ; SRF.snow_beg        = SRF.snow_beg       .where (SRF.snow_beg        < 1E15, 0.)
684    SRF.lakeres_beg     = d_SRF_beg['lakeres']         ; SRF.lakeres_beg     = SRF.lakeres_beg    .where (SRF.lakeres_beg     < 1E15, 0.)
685
686    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.)
687    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.)
688    SRF.snow_end        = d_SRF_end['snow_beg']        ; SRF.snow_end        = SRF.snow_end       .where (SRF.snow_end        < 1E15, 0.)
689    SRF.lakeres_end     = d_SRF_end['lakeres']         ; SRF.lakeres_end     = SRF.lakeres_end    .where (SRF.lakeres_end     < 1E15, 0.)
690
691    if LMDZ :
692        SRF.tot_watveg_beg  = lmdz.geo2point (SRF.tot_watveg_beg , dim1d='cell')
693        SRF.tot_watsoil_beg = lmdz.geo2point (SRF.tot_watsoil_beg, dim1d='cell')
694        SRF.snow_beg        = lmdz.geo2point (SRF.snow_beg       , dim1d='cell')
695        SRF.lakeres_beg     = lmdz.geo2point (SRF.lakeres_beg    , dim1d='cell')
696        SRF.tot_watveg_end  = lmdz.geo2point (SRF.tot_watveg_end , dim1d='cell')
697        SRF.tot_watsoil_end = lmdz.geo2point (SRF.tot_watsoil_end, dim1d='cell')
698        SRF.snow_end        = lmdz.geo2point (SRF.snow_end       , dim1d='cell')
699        SRF.lakeres_end     = lmdz.geo2point (SRF.lakeres_end    , dim1d='cell')
700
701
702    # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood
703
704    SRF.wat_beg = SRF.tot_watveg_beg + SRF.tot_watsoil_beg + SRF.snow_beg
705    SRF.wat_end = SRF.tot_watveg_end + SRF.tot_watsoil_end + SRF.snow_end
706
707    echo ( '\n====================================================================================' )
708    print ('Computing integrals')
709
710    print ( ' 1/8', end='' ) ; sys.stdout.flush ()
711    SRF.mass_watveg_beg   = SRF_stock_int ( SRF.tot_watveg_beg  )
712    print ( ' 2/8', end='' ) ; sys.stdout.flush ()
713    SRF.mass_watsoil_beg  = SRF_stock_int ( SRF.tot_watsoil_beg )
714    print ( ' 3/8', end='' ) ; sys.stdout.flush ()
715    SRF.mass_snow_beg     = SRF_stock_int ( SRF.snow_beg        )
716    print ( ' 4/8', end='' ) ; sys.stdout.flush ()
717    SRF.mass_lake_beg     = ONE_stock_int ( SRF.lakeres_beg     )
718    print ( ' 5/8', end='' ) ; sys.stdout.flush ()
719
720    SRF.mass_watveg_end   = SRF_stock_int ( SRF.tot_watveg_end  )
721    print ( ' 6/8', end='' ) ; sys.stdout.flush ()
722    SRF.mass_watsoil_end  = SRF_stock_int ( SRF.tot_watsoil_end )
723    print ( ' 7/8', end='' ) ; sys.stdout.flush ()
724    SRF.mass_snow_end     = SRF_stock_int ( SRF.snow_end        )
725    print ( ' 8/8', end='' ) ; sys.stdout.flush ()
726    SRF.mass_lake_end     = ONE_stock_int ( SRF.lakeres_end     )
727
728    print (' -- ') ; sys.stdout.flush ()
729
730    SRF.diff_mass_watveg   = Sprec ( [SRF.mass_watveg_end , -SRF.mass_watveg_beg] )
731    SRF.diff_mass_watsoil  = Sprec ( [SRF.mass_watsoil_end, -SRF.mass_watsoil_beg] )
732    SRF.diff_mass_snow     = Sprec ( [SRF.mass_snow_end   , -SRF.mass_snow_beg] )
733    SRF.diff_mass_lake     = Sprec ( [SRF.mass_lake_end   , -SRF.mass_lake_beg] )
734
735    echo ( '------------------------------------------------------------------------------------' )
736    echo ( f'\nSurface reservoirs  -- {Title} ')
737    echo ( f'SRF.mass_watveg_beg   = {SRF.mass_watveg_beg :12.6e} kg | SRF.mass_watveg_end   = {SRF.mass_watveg_end :12.6e} kg ' )
738    echo ( f'SRF.mass_watsoil_beg  = {SRF.mass_watsoil_beg:12.6e} kg | SRF.mass_watsoil_end  = {SRF.mass_watsoil_end:12.6e} kg ' )
739    echo ( f'SRF.mass_snow_beg     = {SRF.mass_snow_beg   :12.6e} kg | SRF.mass_snow_end     = {SRF.mass_snow_end   :12.6e} kg ' )
740    echo ( f'SRF.mass_lake_beg     = {SRF.mass_lake_beg   :12.6e} kg | SRF.mass_lake_end     = {SRF.mass_lake_end   :12.6e} kg ' )
741
742    prtFlux ( 'dMass (watveg) ', SRF.diff_mass_watveg , 'e' , True )
743    prtFlux ( 'dMass (watsoil)', SRF.diff_mass_watsoil, 'e' , True )
744    prtFlux ( 'dMass (snow)   ', SRF.diff_mass_snow   , 'e' , True )
745    prtFlux ( 'dMass (lake)   ', SRF.diff_mass_lake   , 'e' , True )
746
747    SRF.mass_wat_beg = Sprec ( [SRF.mass_watveg_beg , SRF.mass_watsoil_beg, SRF.mass_snow_beg] )
748    SRF.mass_wat_end = Sprec ( [SRF.mass_watveg_end , SRF.mass_watsoil_end, SRF.mass_snow_end] )
749
750    SRF.diff_mass_wat    = Sprec ( [+SRF.mass_watveg_end , +SRF.mass_watsoil_end, +SRF.mass_snow_end,
751                               -SRF.mass_watveg_beg , -SRF.mass_watsoil_beg, -SRF.mass_snow_beg]  )
752
753    echo ( '------------------------------------------------------------------------------------' )
754    echo ( f'Water content in surface  -- {Title} ' )
755    echo ( f'SRF.mass_wat_beg   = {SRF.mass_wat_beg:12.6e} kg | SRF.mass_wat_end  = {SRF.mass_wat_end:12.6e} kg ')
756    prtFlux ( 'dMass (water srf)', SRF.diff_mass_wat, 'e', True )
757
758    echo ( '------------------------------------------------------------------------------------' )
759    echo ( 'Water content in  ATM + SECHIBA + RUN + LAKE' )
760    echo ( 'mass_wat_beg = {:12.6e} kg | mass_wat_end = {:12.6e} kg '.
761            format (DYN.mass_wat_beg + ATM.mass_sno_beg + SRF.RUN_mass_wat_beg + SRF.mass_wat_beg + SRF.mass_lake_beg ,
762                    DYN.mass_wat_end + ATM.mass_sno_end + SRF.RUN_mass_wat_end + SRF.mass_wat_end + SRF.mass_lake_end ) )
763    prtFlux ( 'dMass (water atm+srf+run+lake)', DYN.diff_mass_wat + ATM.diff_mass_sno + SRF.RUN_diff_mass_wat + SRF.diff_mass_wat + SRF.diff_mass_lake, 'e', True)
764
765echo ( '\n====================================================================================' )
766echo ( f'-- ATM Fluxes  -- {Title} ' )
767
768if ATM_HIS == 'latlon' :
769    echo ( ' latlon case' )
770    ATM.wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1d='cell' )
771    ATM.wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1d='cell' )
772    ATM.wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1d='cell' )
773    ATM.wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1d='cell' )
774    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' )
775    ATM.fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1d='cell' )
776    ATM.fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1d='cell' )
777    ATM.precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1d='cell' )
778    ATM.snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1d='cell' )
779    ATM.evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1d='cell' )
780    ATM.wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1d='cell' )
781    ATM.wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1d='cell' )
782    ATM.wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1d='cell' )
783    ATM.wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1d='cell' )
784    ATM.wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1d='cell' )
785    ATM.wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1d='cell' )
786    ATM.wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1d='cell' )
787    ATM.wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1d='cell' )
788    ATM.wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1d='cell' )
789    ATM.wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1d='cell' )
790    ATM.wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1d='cell' )
791    ATM.wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1d='cell' )
792    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' )
793    echo ( 'End of LATLON case')
794
795if ATM_HIS == 'ico' :
796    echo (' ico case')
797    ATM.wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
798    ATM.wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
799    ATM.wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
800    ATM.wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
801    ATM.runofflic   = rprec (d_ATM_his ['runofflic'])
802    ATM.fqcalving   = rprec (d_ATM_his ['fqcalving'])
803    ATM.fqfonte     = rprec (d_ATM_his ['fqfonte']  )
804    ATM.precip      = rprec (d_ATM_his ['precip']   )
805    ATM.snowf       = rprec (d_ATM_his ['snow']     )
806    ATM.evap        = rprec (d_ATM_his ['evap']     )
807    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
808    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
809    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
810    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
811    ATM.runofflic   = rprec (d_ATM_his ['runofflic'])
812    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
813    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
814    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
815    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
816    ATM.wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
817    ATM.wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
818    ATM.wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
819    ATM.wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
820    ATM.wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
821    ATM.wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
822    ATM.wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
823    ATM.wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
824    echo ( 'End of ico case ')
825
826echo ( 'ATM wprecip_oce' )
827ATM.wprecip_oce = ATM.wrain_oce + ATM.wsnow_oce
828ATM.wprecip_ter = ATM.wrain_ter + ATM.wsnow_ter
829ATM.wprecip_sic = ATM.wrain_sic + ATM.wsnow_sic
830ATM.wprecip_lic = ATM.wrain_lic + ATM.wsnow_lic
831
832ATM.wbilo       = ATM.wbilo_oce   + ATM.wbilo_sic   + ATM.wbilo_ter   + ATM.wbilo_lic
833ATM.wevap       = ATM.wevap_oce   + ATM.wevap_sic   + ATM.wevap_ter   + ATM.wevap_lic
834ATM.wprecip     = ATM.wprecip_oce + ATM.wprecip_sic + ATM.wprecip_ter + ATM.wprecip_lic
835ATM.wsnow       = ATM.wsnow_oce   + ATM.wsnow_sic   + ATM.wsnow_ter   + ATM.wsnow_lic
836ATM.wrain       = ATM.wrain_oce   + ATM.wrain_sic   + ATM.wrain_ter   + ATM.wrain_lic
837ATM.wemp        = ATM.wevap - ATM.wprecip
838ATM.emp         = ATM.evap  - ATM.precip
839
840ATM.wprecip_sea = ATM.wprecip_oce + ATM.wprecip_sic
841ATM.wsnow_sea   = ATM.wsnow_oce   + ATM.wsnow_sic
842ATM.wrain_sea   = ATM.wrain_oce   + ATM.wrain_sic
843ATM.wbilo_sea   = ATM.wbilo_oce   + ATM.wbilo_sic
844ATM.wevap_sea   = ATM.wevap_sic   + ATM.wevap_oce
845
846ATM.wemp_ter    = ATM.wevap_ter - ATM.wprecip_ter
847ATM.wemp_oce    = ATM.wevap_oce - ATM.wprecip_oce
848ATM.wemp_sic    = ATM.wevap_sic - ATM.wprecip_sic
849ATM.wemp_lic    = ATM.wevap_lic - ATM.wprecip_lic
850ATM.wemp_sea    = ATM.wevap_sic - ATM.wprecip_oce
851
852if SECHIBA :
853    if RUN_HIS == 'latlon' :
854        echo ( 'RUN costalflow Grille LATLON' )
855        if TestInterp :
856            echo ( 'RUN runoff TestInterp' )
857            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1d='cell' )
858            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1d='cell' )
859        else :
860            echo ( 'RUN runoff' )
861            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1d='cell' )
862            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1d='cell' )
863
864        SRF.RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1d='cell' )
865        SRF.RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1d='cell' )
866        SRF.RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1d='cell' )
867        SRF.RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1d='cell' )
868        SRF.RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1d='cell' )
869
870    if RUN_HIS == 'ico' :
871        echo ( 'RUN costalflow Grille ICO' )
872        SRF.RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])
873        SRF.RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )
874        SRF.RUN_runoff      =  rprec (d_RUN_his ['runoff']     )
875        SRF.RUN_drainage    =  rprec (d_RUN_his ['drainage']   )
876        SRF.RUN_riversret   =  rprec (d_RUN_his ['riversret']  )
877
878        SRF.RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
879        SRF.RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
880
881    Step = 0
882
883    if SRF_HIS == 'latlon' :
884        if TestInterp :
885            echo ( 'SECHIBA rain TestInterp' )
886            SRF.rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1d='cell')
887            SRF.evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1d='cell')
888            SRF.snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1d='cell')
889            SRF.subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1d='cell')
890            SRF.transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1d='cell' )
891            #SRF.rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs )
892            #SRF.evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs )
893            #SRF.snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs )
894            #SRF.subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs )
895            #SRF.transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs )
896        else :
897            echo ( 'SECHIBA rain' )
898            SRF.rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1d='cell')
899            SRF.evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1d='cell')
900            SRF.snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1d='cell')
901            SRF.subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1d='cell')
902            SRF.transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1d='cell' )
903
904    if SRF_HIS == 'ico' :
905        echo ( 'SECHIBA rain')
906        SRF.rain     = rprec (d_SRF_his ['rain'] )
907        SRF.evap     = rprec (d_SRF_his ['evap'] )
908        SRF.snowf    = rprec (d_SRF_his ['snowf'])
909        SRF.subli    = rprec (d_SRF_his ['subli'])
910        SRF.transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget')
911
912    echo ( 'SECHIBA emp' )
913    SRF.transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
914    SRF.emp      = SRF.evap - SRF.rain - SRF.snowf ; SRF.emp.attrs['units'] = SRF.rain.attrs['units']
915
916    ## Correcting units of SECHIBA variables
917    def mmd2SI ( Var ) :
918        '''Change unit from mm/d or m^3/s to kg/s if needed'''
919        if 'units' in Var.attrs : 
920            if Var.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
921                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg/s'
922            if Var.attrs['units'] == 'mm/d' :
923                Var.values = Var.values  * ATM_RHO * (1e-3/86400.) ;  Var.attrs['units'] = 'kg/s'
924            if Var.attrs['units'] in ['m^3', 'm3', ] :
925                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg'
926        return Var
927   
928    SRF.RUN_coastalflow     = mmd2SI ( SRF.RUN_coastalflow )
929    SRF.RUN_coastalflow_cpl = mmd2SI ( SRF.RUN_coastalflow_cpl )
930    SRF.RUN_drainage        = mmd2SI ( SRF.RUN_drainage )
931    SRF.RUN_riverflow       = mmd2SI ( SRF.RUN_riverflow )
932    SRF.RUN_riverflow_cpl   = mmd2SI ( SRF.RUN_riverflow_cpl )
933    SRF.RUN_riversret       = mmd2SI ( SRF.RUN_riversret )
934    SRF.RUN_runoff          = mmd2SI ( SRF.RUN_runoff )
935           
936    SRF.evap     = mmd2SI ( SRF.evap )
937    SRF.snowf    = mmd2SI ( SRF.snowf )
938    SRF.subli    = mmd2SI ( SRF.subli )
939    SRF.transpir = mmd2SI ( SRF.transpir )
940    SRF.rain     = mmd2SI ( SRF.rain )
941    SRF.emp      = mmd2SI ( SRF.emp )
942
943    echo ( 'RUN input' )
944    SRF.RUN_input  = SRF.RUN_runoff      + SRF.RUN_drainage
945    SRF.RUN_output = SRF.RUN_coastalflow + SRF.RUN_riverflow
946
947echo ( 'ATM flw_wbilo' )
948ATM.flx_wbilo       = ATM_flux_int ( ATM.wbilo      )
949ATM.flx_wevap       = ATM_flux_int ( ATM.wevap      )
950ATM.flx_wprecip     = ATM_flux_int ( ATM.wprecip    )
951ATM.flx_wsnow       = ATM_flux_int ( ATM.wsnow      )
952ATM.flx_wrain       = ATM_flux_int ( ATM.wrain      )
953ATM.flx_wemp        = ATM_flux_int ( ATM.wemp       )
954
955ATM.flx_wbilo_lic   = ATM_flux_int ( ATM.wbilo_lic  )
956ATM.flx_wbilo_oce   = ATM_flux_int ( ATM.wbilo_oce  )
957ATM.flx_wbilo_sea   = ATM_flux_int ( ATM.wbilo_sea  )
958ATM.flx_wbilo_sic   = ATM_flux_int ( ATM.wbilo_sic  )
959ATM.flx_wbilo_ter   = ATM_flux_int ( ATM.wbilo_ter  )
960# Type d'integration a verifier
961ATM.flx_calving     = ATM_flux_int ( ATM.fqcalving  )
962ATM.flx_fqfonte     = ATM_flux_int ( ATM.fqfonte    )
963
964ATM.LIC_flx_calving     = LIC_flux_int ( ATM.fqcalving  )
965ATM.LIC_flx_fqfonte     = LIC_flux_int ( ATM.fqfonte    )
966
967echo ( 'ATM flx precip' )
968ATM.flx_precip      = ATM_flux_int ( ATM.precip     )
969ATM.flx_snowf       = ATM_flux_int ( ATM.snowf      )
970ATM.flx_evap        = ATM_flux_int ( ATM.evap       )
971ATM.flx_runlic      = ATM_flux_int ( ATM.runofflic  )
972
973ATM.LIC_flx_precip      = LIC_flux_int ( ATM.precip     )
974ATM.LIC_flx_snowf       = LIC_flux_int ( ATM.snowf      )
975ATM.LIC_flx_evap        = LIC_flux_int ( ATM.evap       )
976ATM.LIC_flx_runlic      = LIC_flux_int ( ATM.runofflic  )
977
978echo ( 'ATM flx_wrain_ter' )
979ATM.flx_wrain_ter    = ATM_flux_int ( ATM.wrain_ter )
980ATM.flx_wrain_oce    = ATM_flux_int ( ATM.wrain_oce )
981ATM.flx_wrain_lic    = ATM_flux_int ( ATM.wrain_lic )
982ATM.flx_wrain_sic    = ATM_flux_int ( ATM.wrain_sic )
983ATM.flx_wrain_sea    = ATM_flux_int ( ATM.wrain_sea )
984
985ATM.flx_wsnow_ter    = ATM_flux_int ( ATM.wsnow_ter )
986ATM.flx_wsnow_oce    = ATM_flux_int ( ATM.wsnow_oce )
987ATM.flx_wsnow_lic    = ATM_flux_int ( ATM.wsnow_lic )
988ATM.flx_wsnow_sic    = ATM_flux_int ( ATM.wsnow_sic )
989ATM.flx_wsnow_sea    = ATM_flux_int ( ATM.wsnow_sea )
990
991echo ( 'ATM flx_evap_ter' )
992ATM.flx_wevap_ter    = ATM_flux_int ( ATM.wevap_ter )
993ATM.flx_wevap_oce    = ATM_flux_int ( ATM.wevap_oce )
994ATM.flx_wevap_lic    = ATM_flux_int ( ATM.wevap_lic )
995ATM.flx_wevap_sic    = ATM_flux_int ( ATM.wevap_sic )
996ATM.flx_wevap_sea    = ATM_flux_int ( ATM.wevap_sea )
997ATM.flx_wprecip_lic  = ATM_flux_int ( ATM.wprecip_lic )
998ATM.flx_wprecip_oce  = ATM_flux_int ( ATM.wprecip_oce )
999ATM.flx_wprecip_sic  = ATM_flux_int ( ATM.wprecip_sic )
1000ATM.flx_wprecip_ter  = ATM_flux_int ( ATM.wprecip_ter )
1001ATM.flx_wprecip_sea  = ATM_flux_int ( ATM.wprecip_sea )
1002ATM.flx_wemp_lic     = ATM_flux_int ( ATM.wemp_lic )
1003ATM.flx_wemp_oce     = ATM_flux_int ( ATM.wemp_oce )
1004ATM.flx_wemp_sic     = ATM_flux_int ( ATM.wemp_sic )
1005ATM.flx_wemp_ter     = ATM_flux_int ( ATM.wemp_ter )
1006ATM.flx_wemp_sea     = ATM_flux_int ( ATM.wemp_sea )
1007
1008ATM.flx_emp          = ATM_flux_int ( ATM.emp )
1009
1010if SECHIBA :
1011    echo ( 'RUN flx_coastal' )
1012    SRF.RUN_flx_coastal     = ONE_flux_int ( SRF.RUN_coastalflow)
1013    echo ( 'RUN flx_river' )
1014    SRF.RUN_flx_river       = ONE_flux_int ( SRF.RUN_riverflow  )
1015    echo ( 'RUN flx_coastal_cpl' )
1016    SRF.RUN_flx_coastal_cpl = ONE_flux_int ( SRF.RUN_coastalflow_cpl)
1017    echo ( 'RUN flx_river_cpl' )
1018    SRF.RUN_flx_river_cpl   = ONE_flux_int ( SRF.RUN_riverflow_cpl  )
1019    echo ( 'RUN flx_drainage' )
1020    SRF.RUN_flx_drainage    = SRF_flux_int ( SRF.RUN_drainage   )
1021    echo ( 'RUN flx_riversset' )
1022    SRF.RUN_flx_riversret   = SRF_flux_int ( SRF.RUN_riversret  )
1023    echo ( 'RUN flx_runoff' )
1024    SRF.RUN_flx_runoff      = SRF_flux_int ( SRF.RUN_runoff     )
1025    echo ( 'RUN flx_input' )
1026    SRF.RUN_flx_input       = SRF_flux_int ( SRF.RUN_input      )
1027    echo ( 'RUN flx_output' )
1028    SRF.RUN_flx_output      = ONE_flux_int ( SRF.RUN_output     )
1029
1030    echo ( 'RUN flx_bil' ) ; Step += 1
1031    #SRF.RUN_flx_bil    = SRF.RUN_flx_input   - SRF.RUN_flx_output
1032    #SRF.RUN_flx_rivcoa = SRF.RUN_flx_coastal + SRF.RUN_flx_river
1033
1034    SRF.RUN_flx_bil    = ONE_flux_int ( SRF.RUN_input       - SRF.RUN_output)
1035    SRF.RUN_flx_rivcoa = ONE_flux_int ( SRF.RUN_coastalflow + SRF.RUN_riverflow)
1036
1037prtFlux ('wbilo_oce            ', ATM.flx_wbilo_oce     , 'f' )
1038prtFlux ('wbilo_sic            ', ATM.flx_wbilo_sic     , 'f' )
1039prtFlux ('wbilo_sic+oce        ', ATM.flx_wbilo_sea     , 'f' )
1040prtFlux ('wbilo_ter            ', ATM.flx_wbilo_ter     , 'f' )
1041prtFlux ('wbilo_lic            ', ATM.flx_wbilo_lic     , 'f' )
1042prtFlux ('Sum wbilo_*          ', ATM.flx_wbilo         , 'f', True)
1043prtFlux ('E-P                  ', ATM.flx_emp           , 'f', True)
1044prtFlux ('calving              ', ATM.flx_calving       , 'f' )
1045prtFlux ('fqfonte              ', ATM.flx_fqfonte       , 'f' )
1046prtFlux ('precip               ', ATM.flx_precip        , 'f' )
1047prtFlux ('snowf                ', ATM.flx_snowf         , 'f' )
1048prtFlux ('evap                 ', ATM.flx_evap          , 'f' )
1049prtFlux ('runoff lic           ', ATM.flx_runlic        , 'f' )
1050
1051prtFlux ('ATM.flx_wevap*       ', ATM.flx_wevap         , 'f' )
1052prtFlux ('ATM.flx_wrain*       ', ATM.flx_wrain         , 'f' )
1053prtFlux ('ATM.flx_wsnow*       ', ATM.flx_wsnow         , 'f' )
1054prtFlux ('ATM.flx_wprecip*     ', ATM.flx_wprecip       , 'f' )
1055prtFlux ('ATM.flx_wemp*        ', ATM.flx_wemp          , 'f', True )
1056
1057echo ( 'Errors <field> vs. wbil_<field>' )
1058prtFlux ('ERROR evap           ', ATM.flx_wevap   - ATM.flx_evap  , 'e', True )
1059prtFlux ('ERROR precip         ', ATM.flx_wprecip - ATM.flx_precip, 'e', True )
1060prtFlux ('ERROR snow           ', ATM.flx_wsnow   - ATM.flx_snowf , 'e', True )
1061prtFlux ('ERROR emp            ', ATM.flx_wemp    - ATM.flx_emp   , 'e', True )
1062
1063if SECHIBA :
1064    echo ( '\n====================================================================================' )
1065    echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
1066    prtFlux ('coastalflow   ', SRF.RUN_flx_coastal    , 'f' )
1067    prtFlux ('riverflow     ', SRF.RUN_flx_river      , 'f' )
1068    prtFlux ('coastal_cpl   ', SRF.RUN_flx_coastal_cpl, 'f' )
1069    prtFlux ('riverf_cpl    ', SRF.RUN_flx_river_cpl  , 'f' )
1070    prtFlux ('river+coastal ', SRF.RUN_flx_rivcoa     , 'f' )
1071    prtFlux ('drainage      ', SRF.RUN_flx_drainage   , 'f' )
1072    prtFlux ('riversret     ', SRF.RUN_flx_riversret  , 'f' )
1073    prtFlux ('runoff        ', SRF.RUN_flx_runoff     , 'f' )
1074    prtFlux ('river in      ', SRF.RUN_flx_input      , 'f' )
1075    prtFlux ('river out     ', SRF.RUN_flx_output     , 'f' )
1076    prtFlux ('river bil     ', SRF.RUN_flx_bil        , 'f' )
1077
1078ATM.flx_budget = -ATM.flx_wbilo + ATM.flx_calving + ATM.flx_runlic #+# ATM.flx_fqfonte #+ SRF.RUN_flx_river
1079
1080
1081echo ('')
1082#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 ))
1083
1084#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 ))
1085
1086ATM.flx_toSRF = -ATM.flx_wbilo_ter
1087
1088echo (' ')
1089echo ( '\n====================================================================================' )
1090echo ( f'--  Atmosphere  -- {Title} ' )
1091echo ( f'Mass begin = {DYN.mass_wat_beg:12.6e} kg | Mass end = {DYN.mass_wat_end:12.6e} kg' )
1092prtFlux ( 'dmass (atm) ', DYN.diff_mass_wat , 'e', True )
1093prtFlux ( 'Sum wbilo_* ', ATM.flx_wbilo, 'e', True )
1094prtFlux ( 'E-P         ', ATM.flx_emp  , 'e', True )
1095echo ( ' ' )
1096prtFlux ( 'Water loss atm from wbil_*', ATM.flx_wbilo - DYN.diff_mass_wat, 'f', True )
1097echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM.flx_wbilo - DYN.diff_mass_wat)/DYN.diff_mass_wat  ) )
1098
1099echo (' ')
1100prtFlux ( 'Water loss atm from E-P', ATM.flx_emp  - DYN.diff_mass_wat , 'f', True )
1101echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM.flx_emp-DYN.diff_mass_wat)/DYN.diff_mass_wat  ) )
1102echo (' ')
1103
1104ATM.error =  ATM.flx_emp  - DYN.diff_mass_wat
1105
1106
1107echo (' ')
1108echo ( '\n====================================================================================' )
1109
1110ATM.LIC_flx_budget1 = Sprec ( [-ATM.flx_wemp_lic  , -ATM.LIC_flx_calving , -ATM.LIC_flx_fqfonte] )
1111ATM.LIC_flx_budget2 = Sprec ( [-ATM.flx_wbilo_lic , -ATM.LIC_flx_calving , -ATM.LIC_flx_fqfonte] )
1112ATM.LIC_flx_budget3 = Sprec ( [-ATM.flx_wbilo_lic , -ATM.LIC_flx_runlic] )
1113ATM.LIC_flx_budget4 = Sprec ( [-ATM.flx_wemp_lic  , -ATM.LIC_flx_runlic] )
1114
1115echo ( f'--  LIC  -- {Title} ' )
1116echo ( f'Mass total   begin = {ATM.LIC_mass_wat_beg    :12.6e} kg | Mass end = {ATM.LIC_mass_wat_end    :12.6e} kg' )
1117echo ( f'Mass snow    begin = {ATM.LIC_mass_sno_beg    :12.6e} kg | Mass end = {ATM.LIC_mass_sno_end    :12.6e} kg' )
1118echo ( f'Mass qs      begin = {ATM.LIC_mass_qs_beg     :12.6e} kg | Mass end = {ATM.LIC_mass_qs_end     :12.6e} kg' )
1119echo ( f'Mass runlic0 begin = {ATM.LIC_mass_runlic0_beg:12.6e} kg | Mass end = {ATM.LIC_mass_runlic0_end:12.6e} kg' )
1120prtFlux ( 'dmass (LIC sno)       ', ATM.LIC_diff_mass_sno          , 'f', True, width=45 )
1121prtFlux ( 'dmass (LIC qs)        ', ATM.LIC_diff_mass_qs           , 'e', True, width=45 )
1122prtFlux ( 'dmass (LIC wat)       ', ATM.LIC_diff_mass_wat          , 'f', True, width=45 )
1123prtFlux ( 'dmass (LIC runlic0)   ', ATM.LIC_diff_mass_runlic0      , 'e', True, width=45 )
1124prtFlux ( 'dmass (LIC total)     ', ATM.LIC_diff_mass_wat          , 'e', True, width=45 )
1125prtFlux ( 'LIC ATM.flx_wemp_lic  ', ATM.flx_wemp_lic          , 'f', True, width=45 )
1126prtFlux ( 'LIC LIC_flx_fqfonte   ', ATM.LIC_flx_fqfonte       , 'f', True, width=45 )
1127prtFlux ( 'LIC LIC_flx_calving   ', ATM.LIC_flx_calving       , 'f', True, width=45 )
1128prtFlux ( 'LIC LIC_flx_runofflic ', ATM.LIC_flx_runlic        , 'f', True, width=45 )
1129prtFlux ( 'LIC fqfonte + calving ', ATM.LIC_flx_calving+ATM.LIC_flx_fqfonte , 'f', True, width=45 )
1130prtFlux ( 'LIC fluxes 1 ( wemp_lic   - fqcalving - fqfonte)) ', ATM.LIC_flx_budget1              , 'f', True, width=45 )
1131prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)   ', ATM.LIC_flx_budget2              , 'f', True, width=45 )
1132prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)    ', ATM.LIC_flx_budget3              , 'f', True, width=45 )
1133prtFlux ( 'LIC fluxes 3 ( wemp_lic  - runofflic*frac_lic)    ', ATM.LIC_flx_budget4              , 'f', True, width=45 )
1134prtFlux ( 'LIC error 1                                       ', ATM.LIC_flx_budget1-ATM.LIC_diff_mass_wat , 'e', True, width=45 )
1135prtFlux ( 'LIC error 2                                       ', ATM.LIC_flx_budget2-ATM.LIC_diff_mass_wat , 'e', True, width=45 )
1136prtFlux ( 'LIC error 3                                       ', ATM.LIC_flx_budget3-ATM.LIC_diff_mass_wat , 'e', True, width=45 )
1137echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget1-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) )
1138echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget2-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) )
1139echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget3-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) )
1140
1141if SECHIBA :
1142    echo ( '\n====================================================================================' )
1143    echo ( f'-- SECHIBA fluxes  -- {Title} ' )
1144
1145    SRF.flx_rain     = SRF_flux_int ( SRF.rain     )
1146    SRF.flx_evap     = SRF_flux_int ( SRF.evap     )
1147    SRF.flx_snowf    = SRF_flux_int ( SRF.snowf    )
1148    SRF.flx_subli    = SRF_flux_int ( SRF.subli    )
1149    SRF.flx_transpir = SRF_flux_int ( SRF.transpir )
1150    SRF.flx_emp      = SRF_flux_int ( SRF.emp      )
1151
1152    SRF.RUN_flx_torouting   = SRF_flux_int  ( SRF.RUN_runoff      + SRF.RUN_drainage)
1153    SRF.RUN_flx_fromrouting = ONE_flux_int  ( SRF.RUN_coastalflow + SRF.RUN_riverflow )
1154
1155    SRF.flx_all =  SRF_flux_int  ( SRF.rain + SRF.snowf - SRF.evap - SRF.RUN_runoff - SRF.RUN_drainage )
1156
1157    prtFlux ('rain         ', SRF.flx_rain       , 'f' )
1158    prtFlux ('evap         ', SRF.flx_evap       , 'f' )
1159    prtFlux ('snowf        ', SRF.flx_snowf      , 'f' )
1160    prtFlux ('E-P          ', SRF.flx_emp        , 'f' )
1161    prtFlux ('subli        ', SRF.flx_subli      , 'f' )
1162    prtFlux ('transpir     ', SRF.flx_transpir   , 'f' )
1163    prtFlux ('to routing   ', SRF.RUN_flx_torouting  , 'f' )
1164    prtFlux ('budget       ', SRF.flx_all        , 'f', small=True )
1165
1166    echo ( '\n------------------------------------------------------------------------------------' )
1167    echo ( 'Water content in surface ' )
1168    echo ( f'SRF.mass_wat_beg  = {SRF.mass_wat_beg:12.6e} kg | SRF.mass_wat_end  = {SRF.mass_wat_end:12.6e} kg ' )
1169    prtFlux ( 'dMass (water srf)', SRF.diff_mass_wat, 'e', small=True)
1170    prtFlux ( 'Error            ', SRF.flx_all-SRF.diff_mass_wat, 'e', small=True )
1171    echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF.flx_all-SRF.diff_mass_wat)/SRF.diff_mass_wat) )
1172
1173    echo ( '\n====================================================================================' )
1174    echo ( f'-- Check ATM vs. SECHIBA -- {Title} ' )
1175    prtFlux ('E-P ATM       ', ATM.flx_wemp_ter   , 'f' )
1176    prtFlux ('wbilo ter     ', ATM.flx_wbilo_ter  , 'f' )
1177    prtFlux ('E-P SECHIBA   ', SRF.flx_emp        , 'f' )
1178    prtFlux ('SRF/ATM error ', ATM.flx_wbilo_ter - SRF.flx_emp, 'e', True)
1179    echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM.flx_wbilo_ter - SRF.flx_emp)/SRF.flx_emp ) )
1180
1181    echo ('')
1182    echo ( '\n====================================================================================' )
1183    echo ( f'-- RUNOFF fluxes  -- {Title} ' )
1184    SRF.RUN_flx_all = SRF.RUN_flx_torouting - SRF.RUN_flx_river - SRF.RUN_flx_coastal
1185    prtFlux ('runoff    ', SRF.RUN_flx_runoff      , 'f' )
1186    prtFlux ('drainage  ', SRF.RUN_flx_drainage    , 'f' )
1187    prtFlux ('run+drain ', SRF.RUN_flx_torouting   , 'f' )
1188    prtFlux ('river     ', SRF.RUN_flx_river       , 'f' )
1189    prtFlux ('coastal   ', SRF.RUN_flx_coastal     , 'f' )
1190    prtFlux ('riv+coa   ', SRF.RUN_flx_fromrouting , 'f' )
1191    prtFlux ('budget    ', SRF.RUN_flx_all         , 'f' , small=True)
1192
1193    echo ( '\n------------------------------------------------------------------------------------' )
1194    echo ( f'Water content in routing+lake  -- {Title} ' )
1195    echo ( f'SRF.RUN_mass_wat_beg  = {SRF.RUN_mass_wat_beg:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg ' )
1196    prtFlux ( 'dMass (routing)  ', SRF.RUN_diff_mass_wat+SRF.diff_mass_lake, 'f', small=True)
1197    prtFlux ( 'Routing error    ', SRF.RUN_flx_all+SRF.diff_mass_lake-SRF.RUN_diff_mass_wat, 'e', small=True )
1198    echo ( 'Routing error : {:12.3e} (rel)'.format ( (SRF.RUN_flx_all-SRF.diff_mass_lake-SRF.RUN_diff_mass_wat)/(SRF.RUN_diff_mass_wat+SRF.diff_mass_lake) ) )
1199
1200    echo ( '\n------------------------------------------------------------------------------------' )
1201    echo ( f'Water content in routing  -- {Title} ' )
1202    echo ( f'SRF.RUN_mass_wat_beg  = {SRF.RUN_mass_wat_beg:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg ' )
1203    prtFlux ( 'dMass (routing) ', SRF.RUN_diff_mass_wat, 'f', small=True  )
1204    prtFlux ( 'Routing error   ', SRF.RUN_flx_all-SRF.RUN_diff_mass_wat, 'e', small=True)
1205    echo ( 'Routing error : {:12.3e} (rel)'.format ( (SRF.RUN_flx_all-SRF.RUN_diff_mass_wat)/SRF.RUN_diff_mass_wat ) )
1206
1207echo ( ' ' )
1208echo ( f'{Title = }' )
1209
1210echo ( 'SVN Information' )
1211for kk in SVN.keys():
1212    print ( SVN[kk] )
1213
1214## Write the full configuration
1215params_out = open (FullIniFile, 'w', encoding = 'utf-8')
1216params = wu.dict2config ( dpar )
1217params.write ( params_out )
1218params_out.close ()
Note: See TracBrowser for help on using the repository browser.