source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6508

Last change on this file since 6508 was 6508, checked in by omamce, 11 months ago

WaterUtils?.pyVersion qui marche pour :
Water Budget
O.M.

Version qui marche pour :

  • Grille LMDZ et routage SECHIBA
  • Grille ICO avec sorties natives et routage SIMPLE : routage pas très précis.

Ne marche pas pour :

  • Grille LMDZ et routage SIMPLE : pb sur runoff
  • Grille ICO avec sorties interpolées :

# Erreurs relatives

### VALID-CM622-LR.01 - LMDZ : OK ? Sauf LIC

  • LIC : 1.e-2
  • SRF : 4.e-6
  • SRF/ATM : 2e-10
  • RUNOFF : 1.5e-4

### VALID-CM622-SIMPLE-ROUTING - LMDZ : Erreur RUNOFF, LIC

  • LIC : 1.6
  • SRF : 1e-6
  • SRF/ATM : 1.e-10
  • RUNOFF : -7

### TEST-CM72-SIMPLE-ROUTING.13 - ICO : Erreur SRF, RUNOFF, LIC

  • LIC : 150
  • SRF : 0.5
  • SRF/ATM : 4e-2
  • RUNOFF : 700

### CM71v420-LR-pd-02.ini - ICO interpolé : Erreur SRF, RUNOFF, LIC

  • LIC : 15
  • SRF : 0.7
  • SRF/ATM : -5.e-2
  • ROUTING : 3e3

### CM71v420-LR-pd-02.ini - ICO natif : Erreurs faibles RUNOFF, LIC

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