source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6652

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

O.M. : TOOLS/MOSAIX - Cosmetic changes
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 45.5 KB
Line 
1#!/usr/bin/env python3
2###
3### Script to check water conservation in NEMO
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##
14SVN = {
15    'Author'  : "$Author$",
16    'Date'    :   "$Date$",
17    'Revision': "$Revision$",
18    'Id'      : "$Id: ATM_waterbudget.py 6508 2023-06-13 10:58:38Z omamce $",
19    'HeadURL' : "$HeadUrl: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/WATER_BUDGET/ATM_waterbudget.py $"
20    }
21
22###
23## Import system modules
24import sys, os, shutil#, subprocess, platform
25import configparser, re
26from pathlib import Path
27
28## Import needed scientific modules
29import numpy as np, xarray as xr
30
31# Check python version
32if sys.version_info < (3, 8, 0) :
33    print ( f'Python version : {platform.python_version()}' ) 
34    raise Exception ( "Minimum Python version is 3.8" )
35
36## Import local module
37import WaterUtils as wu
38import libIGCM_sys
39import nemo, lmdz
40
41from WaterUtils import VarInt, Rho, Ra, Grav, ICE_rho_ice, ICE_rho_sno, OCE_rho_liq, ATM_rho, SRF_rho, RUN_rho, ICE_rho_pnd, YearLength
42
43## Creates parser for reading .ini input file
44## -------------------------------------------
45config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
46config.optionxform = str # To keep capitals
47
48## Experiment parameters
49## ---------------------
50ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
51OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None
52TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
53YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None
54
55##
56ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
57FileDir=None ; FileOut=None
58dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
59FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
60file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
61tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
62file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
63file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
64file_ICE_beg=None ; file_OCE_beg=None
65file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
66tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
67tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
68tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
69tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
70ContinueOnError=False ; ErrorCount=0
71
72##
73## Precision of history file reading
74## ---------------------------------
75# Default is float (full precision). Degrade the precision by using np.float32
76# Restart file are always read at the full precision
77readPrec=float
78
79## Read command line arguments
80## ---------------------------
81print ( "Name of Python script:", sys.argv[0] )
82IniFile = sys.argv[1]
83
84# Test existence of IniFile
85if not os.path.exists (IniFile ) :
86    raise FileExistsError ( f"File not found : {IniFile = }" )
87
88if 'full' in IniFile : FullIniFile = IniFile
89else                 : FullIniFile = 'full_' + IniFile
90   
91print ("Input file : ", IniFile )
92config.read (IniFile)
93FullIniFile = 'full_' + IniFile
94
95# Reading config.card if possible
96## -------------------------------
97ConfigCard = None
98
99if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
100    if 'ConfigCard' in config['Experiment'].keys () :
101        ConfigCard = config['Experiment']['ConfigCard']
102        print ( f'{ConfigCard=}' )
103
104if ConfigCard : ## Read config card if it exists
105    # Text existence of ConfigCard
106    if os.path.exists ( ConfigCard ) :
107        print ( f'Reading Config Card : {ConfigCard}' )
108        ## Creates parser for reading .ini input file
109        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
110        MyReader.optionxform = str # To keep capitals
111       
112        MyReader.read (ConfigCard)
113
114        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
115            if VarName in MyReader['UserChoices'].keys() :
116                locals()[VarName] = MyReader['UserChoices'][VarName]
117                exec  ( f'{VarName} = wu.setBool ({VarName})' )
118                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
119                exec  ( f'{VarName} = wu.setNone ({VarName})' )
120                exec  ( f'wu.{VarName} = {VarName}' )
121                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
122               
123        for VarName in ['PackFrequency'] :
124            if VarName in MyReader['Post'].keys() :
125                locals()[VarName] = MyReader['Post'][VarName]
126                exec  ( f'{VarName} = wu.setBool ({VarName})' )
127                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
128                exec  ( f'{VarName} = wu.setNone ({VarName})' )
129                exec  ( f'wu.{VarName} = {VarName}' )
130                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
131    else :
132        raise FileExistsError ( f"File not found : {ConfigCard = }" )
133
134## Reading config file
135## -------------------
136for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
137   if Section in config.keys () : 
138        print ( f'\nReading [{Section}]' )
139        for VarName in config[Section].keys() :
140            locals()[VarName] = config[Section][VarName]
141            exec ( f'{VarName} = wu.setBool ({VarName})' )
142            exec ( f'{VarName} = wu.setNum  ({VarName})' )
143            exec ( f'{VarName} = wu.setNone ({VarName})' )
144            exec ( f'wu.{VarName} = {VarName}' )
145            print ( f'    {VarName:21} set to : {locals()[VarName]}' )
146
147print ( f'\nConfig file readed : {IniFile} ' )
148
149##
150## Reading prec
151if wu.unDefined ( 'readPrec' )  :
152    readPrec = np.float64
153else :
154    if readPrec in ["float", "float64", "r8", "double", "<class 'float'>"         ] : readPrec = float
155    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32
156    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16
157             
158## Some physical constants
159## =======================
160if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m)
161if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s
162if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3
163if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3
164if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO
165if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3)
166if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3)
167if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3)
168if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3)
169if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s)
170
171## Set libIGCM and machine dependant values
172## ----------------------------------------
173if not 'Files' in config.keys () : config['Files'] = {}
174
175config['Physics'] = { 'Ra':str(Ra), 'Grav':str(Grav), 'ICE_rho_ice':str(ICE_rho_ice), 'ICE_rho_sno':str(ICE_rho_sno),
176                      'OCE_rho_liq':str(OCE_rho_liq), 'ATM_rho':str(ATM_rho), 'SRF_rho':str(SRF_rho), 'RUN_rho':str(RUN_rho)}
177
178config['Config'] = { 'ContinueOnError':ContinueOnError, 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) }
179
180## --------------------------
181ICO  = ( 'ICO' in wu.ATM )
182LMDZ = ( 'LMD' in wu.ATM )
183
184mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
185             ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None, 
186             R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None )
187globals().update(mm)
188
189config['Files']['TmpDir'] = TmpDir
190config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild }
191
192## Defines begining and end of experiment
193## --------------------------------------
194if wu.unDefined ( 'DateBegin' ) :
195    DateBegin = f'{YearBegin}0101'
196    config['Experiment']['DateBegin'] = DateBegin
197else :
198    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin )
199    DateBegin = wu.FormatToGregorian ( DateBegin)
200    config['Experiment']['YearBegin'] = str(YearBegin)
201
202if wu.unDefined ( 'DateEnd' )   :
203    DateEnd   = f'{YearEnd}1231'
204    config['Experiment']['DateEnd']   = str(DateEnd)
205else :
206    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd )
207    DateEnd   = wu.FormatToGregorian ( DateEnd)
208    config['Experiment']['DateEnd']   = str(DateEnd)
209
210if wu.unDefined ( 'PackFrequency' ) :
211    PackFrequency = YearEnd - YearBegin +1
212    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
213
214if type ( PackFrequency ) == str : 
215    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
216    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
217    PackFrequency = int ( PackFrequency )
218
219print ( f'{YearBegin=} {DateBegin=}' )
220print ( f'{YearEnd  =} {DateEnd  =}' )
221print ( f'{PackFrequency=}' )
222   
223## Output file with water budget diagnostics
224## -----------------------------------------
225if wu.unDefined ( 'FileOut' ) :
226    FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}'
227    if readPrec == np.float32 : FileOut = f'{FileOut}_float32'
228    FileOut = f'{FileOut}.out'
229
230    config['Files']['FileOut'] = FileOut
231
232f_out = open ( FileOut, mode = 'w' )
233
234## Useful functions
235## ----------------
236
237# Degrades precision
238if readPrec == float :
239    def rprec (tab) : return tab
240else : 
241    def rprec (tab) : return tab.astype(readPrec).astype(float)
242       
243def kg2Sv  (val, rho=ATM_rho) :
244    '''From kg to Sverdrup'''
245    return val/dtime_sec*1.0e-6/rho
246
247def kg2myear (val, rho=ATM_rho) :
248    '''From kg to m/year'''
249    return val/OCE_aire_tot/rho/NbYear
250
251def var2prt (var, small=False, rho=ATM_rho) :
252    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
253    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
254
255def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
256    if small :
257        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
258        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
259    else :
260        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
261        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
262    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) )
263    return None
264
265def echo (string, end='\n') :
266    '''Function to print to stdout *and* output file'''
267    print ( str(string), end=end  )
268    sys.stdout.flush ()
269    f_out.write ( str(string) + end )
270    f_out.flush ()
271    return None
272   
273   
274## Set libIGCM directories
275## -----------------------
276if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
277if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
278if wu.unDefined ('L_EXP'      ) : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName)
279if wu.unDefined ('R_SAVE'     ) : R_SAVE      = os.path.join ( R_OUT, L_EXP )
280if wu.unDefined ('R_BUFR'     ) : R_BUFR      = os.path.join ( R_BUF, L_EXP )
281if wu.unDefined ('POST_DIR'   ) : POST_DIR    = os.path.join ( R_BUFR, 'Out' )
282if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' )
283if wu.unDefined ('R_BUF_KSH'  ) : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
284if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
285
286config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
287                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
288
289## Set directory to extract files
290## ------------------------------
291if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
292config['Files']['FileDir'] = FileDir
293
294if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
295
296##- Set directories to rebuild ocean and ice restart files
297if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
298if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
299if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
300if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
301
302echo (' ')
303echo ( f'JobName    : {JobName}'   )
304echo ( f'Comment    : {Comment}'   )
305echo ( f'TmpDir     : {TmpDir}'    )
306echo ( f'FileDir    : {FileDir}'    )
307echo ( f'FileDirOCE : {FileDirOCE}' )
308echo ( f'FileDirICE : {FileDirICE}' )
309
310echo ( f'\nDealing with {L_EXP}'  )
311
312## Creates model output directory names
313## ------------------------------------
314if Freq == "MO" : FreqDir = os.path.join ( 'Output' , 'MO' )
315if Freq == "SE" : FreqDir = os.path.join ( 'Analyse', 'SE' )
316if wu.unDefined ( 'dir_OCE_his' ) :
317    dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir )
318    config['Files']['dir_OCE_his'] = dir_OCE_his
319if wu.unDefined ( 'dir_ICE_his' ) :
320    dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir )
321    config['Files']['dir_OCE_his'] = dir_OCE_his
322 
323echo ( f'The analysis relies on files from the following model output directories : ' )
324echo ( f'{dir_OCE_his}' )
325echo ( f'{dir_ICE_his}' )
326
327#-- Creates files names
328if wu.unDefined ( 'Period ' ) :
329    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
330    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
331    config['Files']['Period'] = Period
332
333config['Files']['DateBegin'] = DateBegin
334config['Files']['DateBegin'] = DateEnd
335       
336echo ( f'Period   : {Period}' )
337
338if wu.unDefined ( 'FileCommon' ) :
339    FileCommon = f'{JobName}_{Period}'
340    config['Files']['FileCommon'] = FileCommon
341
342if wu.unDefined ( 'Title' ) :
343    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
344    config['Files']['Title'] = Title
345       
346echo ('\nOpen history files' )
347if wu.unDefined ('file_OCE_his' ) :
348    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' )
349    file_OCE_his = file_OCE_his
350if wu.unDefined ('file_OCE_sca' ) :   
351    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' )
352    config['Files']['file_OCE_sca'] = file_OCE_sca
353if wu.unDefined ( 'file_ICE_hi' ) : 
354    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' )
355    config['Files']['file_ICE_his'] = file_ICE_his
356
357d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
358d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
359d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
360if NEMO == 3.6 :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
361
362echo ( f'{file_OCE_his = }' )
363echo ( f'{file_OCE_sca = }' )
364echo ( f'{file_ICE_his = }' )
365
366## Compute run length
367## ------------------
368dtime = ( d_OCE_his.time_counter_bounds.max() - d_OCE_his.time_counter_bounds.min() )
369echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' )
370dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
371
372## Compute length of each period
373## -----------------------------
374dtime_per = ( d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] )
375echo ( f'\nPeriods lengths (days) : ')
376echo ( f' {(dtime_per/np.timedelta64(1, "D")).values}' )
377dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
378dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_OCE_his.time_counter,] )
379dtime_per_sec.attrs['unit'] = 's'
380
381## Number of years
382NbYear = dtime_sec / YearLength
383
384##-- Extract restart files from tar
385
386if wu.unDefined ('TarRestartDate_beg' ) : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin )
387if wu.unDefined ('TarRestartDate_end' ) : TarRestartDate_end = wu.FormatToGregorian ( DateEnd )
388
389if wu.unDefined ( 'TarRestartPeriod_beg' ) :
390
391    TarRestartPeriod_beg_DateEnd = TarRestartDate_beg
392    TarRestartPeriod_beg_DateBeg = wu.DateAddYear ( TarRestartPeriod_beg_DateEnd, -PackFrequency )
393    TarRestartPeriod_beg_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_beg_DateBeg )
394   
395    TarRestartPeriod_beg = f'{TarRestartPeriod_beg_DateBeg}_{TarRestartPeriod_beg_DateEnd}'
396    echo (f'Tar period for initial restart : {TarRestartPeriod_beg}')
397    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
398
399if wu.unDefined ( 'TarRestartPeriod_end' ) :
400
401    TarRestartPeriod_end_DateEnd = TarRestartDate_end
402    TarRestartPeriod_end_DateBeg = wu.DateAddYear ( TarRestartPeriod_end_DateEnd, -PackFrequency )
403    TarRestartPeriod_end_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_end_DateBeg )
404     
405    TarRestartPeriod_end = f'{TarRestartPeriod_end_DateBeg}_{TarRestartPeriod_end_DateEnd}'
406    echo (f'Tar period for final restart : {TarRestartPeriod_end}')
407    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
408
409echo (f'Restart dates - Start : {TarRestartPeriod_beg}  /  End : {TarRestartPeriod_end}')
410
411if wu.unDefined ( 'tar_restart_beg' ) :
412    tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' )
413    config['Files']['tar_restart_beg'] = tar_restart_beg
414if wu.unDefined ( 'tar_restart_end' ) :
415    tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' )
416    config['Files']['tar_restart_end'] = tar_restart_end
417
418echo ( f'{tar_restart_beg = }' )
419echo ( f'{tar_restart_end = }' )
420   
421#
422if wu.unDefined ('tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg
423if wu.unDefined ('tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg
424if wu.unDefined ('tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg
425if wu.unDefined ('tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg
426if wu.unDefined ('tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg
427if wu.unDefined ('tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg
428
429if wu.unDefined ('tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end
430if wu.unDefined ('tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end
431if wu.unDefined ('tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end
432if wu.unDefined ('tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end
433if wu.unDefined ('tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end
434if wu.unDefined ('tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end
435
436if wu.unDefined ( 'file_OCE_beg' ) : 
437    file_OCE_beg = f'{FileDir}/OCE_{JobName}_{TarRestartDate_beg}_restart.nc'
438    config['Files']['file_OCE_beg'] = file_OCE_beg
439if wu.unDefined ( 'file_OCE_end' ) :
440    file_OCE_end = f'{FileDir}/OCE_{JobName}_{TarRestartDate_end}_restart.nc'
441    config['Files']['file_OCE_end'] = file_OCE_end
442if wu.unDefined ( 'file_ICE_beg' ) :
443    file_ICE_beg = f'{FileDir}/ICE_{JobName}_{TarRestartDate_beg}_restart_icemod.nc'
444    config['Files']['file_ICE_beg'] = file_ICE_beg
445if wu.unDefined ( 'file_ICE_end' ) : 
446    file_ICE_end = f'{FileDir}/ICE_{JobName}_{TarRestartDate_end}_restart_icemod.nc'
447    config['Files']['file_ICE_end'] = file_ICE_end
448
449echo ( f'{file_OCE_beg}' )
450echo ( f'{file_OCE_end}' )
451echo ( f'{file_ICE_beg}' )
452echo ( f'{file_ICE_end}' )
453
454liste_beg = [file_OCE_beg, file_ICE_beg ]
455liste_end = [file_ATM_end, file_ICE_end ]
456
457echo ('\nExtract and rebuild OCE and ICE restarts')
458def get_ndomain (zfile) :
459    #d_zfile = xr.open_dataset (zfile, decode_times=False).squeeze()
460    #ndomain_opa = d_zfile.attrs['DOMAIN_number_total']
461    #d_zfile.close ()
462    ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) #.format()
463    return int (ndomain_opa)
464
465def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, FileDirComp=FileDirOCE, ErrorCount=ErrorCount ) :
466    '''Extract restart file from tar. Rebuild ocean files if needed'''
467    echo ( f'----------')
468    if os.path.exists ( file_name ) :
469        echo ( f'-- File ready : {file_name = }' )
470    else : 
471        echo ( f'-- Extracting {file_name = }' )
472        base_resFile = Path (file_name).stem # basename, and remove suffix
473        # Try to extract the rebuilded file
474        if os.path.exists ( tar_restart ) :
475            command =  f'cd {FileDirComp} ; tar xf {tar_restart} {base_resFile}.nc'
476            echo ( command )
477            try : 
478                os.system ( command )
479            except :
480                if not os.path.exists ( os.path.join (FileDir, f'{base_resFile}_0000.nc') ):
481                    command =  f'cd {FileDirComp} ; tar xf {tar_restart_end} {base_file}_*.nc'
482                    echo ( command )
483                    ierr = os.system ( command )
484                    if ierr == 0 : echo ( f'tar done : {base_resFile}.nc')
485                    else         : raise Exception ( f'command failed : {command}' )
486                    echo ( f'extract ndomain' )
487                ndomain_opa = get_ndomain ( os.path.join (FileDir, f'{base_file}') )
488                command = f'cd {FileDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {FileDir}'
489                echo ( command )
490                ierr = os.system ( command )
491                if ierr == 0 : echo ( f'Rebuild done : {base_resFile}.nc')
492                else         : raise Exception ( f'command failed : {command}' )
493            else : 
494                echo ( f'tar done : {base_resFile}')
495                command = f'cd {FileDirComp} ; mv {base_resFile}.nc {FileDir}'
496                ierr = os.system ( command )
497                if ierr == 0 : echo ( f'command done : {command}' )
498                else         : raise Exception ( f'command failed : {command = }' )                   
499        else :
500            echo ( f'****** Tar restart file {tar_restart = } not found ' )
501            if ContinueOnError :
502                ErrorCount += 1
503            else : 
504                raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' )
505    return ErrorCount
506
507           
508ErrorCount += extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, FileDirComp=FileDirOCE )
509ErrorCount += extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, FileDirComp=FileDirOCE )
510ErrorCount += extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, FileDirComp=FileDirICE )
511ErrorCount += extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, FileDirComp=FileDirICE )
512
513echo ('Opening OCE and ICE restart files')
514if NEMO == 3.6 : 
515    d_OCE_beg = xr.open_dataset ( os.path.join (FileDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
516    d_OCE_end = xr.open_dataset ( os.path.join (FileDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()
517    d_ICE_beg = xr.open_dataset ( os.path.join (FileDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()
518    d_ICE_end = xr.open_dataset ( os.path.join (FileDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()
519if NEMO == 4.0 or NEMO == 4.2 : 
520    d_OCE_beg = xr.open_dataset ( os.path.join (FileDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
521    d_OCE_end = xr.open_dataset ( os.path.join (FileDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
522    d_ICE_beg = xr.open_dataset ( os.path.join (FileDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
523    d_ICE_end = xr.open_dataset ( os.path.join (FileDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
524
525## Write the full configuration
526config_out = open (FullIniFile, 'w')
527config.write (config_out )
528config_out.close ()
529   
530# Get mask and surfaces
531sos = d_OCE_his ['sos'][0].squeeze()
532OCE_msk = nemo.lbc_mask ( xr.where ( sos>0., 1., 0. ), cd_type = 'T', sval = 0. )
533
534so = sos = d_OCE_his ['sos'][0].squeeze()
535OCE_msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', sval = 0. )
536
537# lbc_mask removes the duplicate points (periodicity and north fold)
538OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE_msk, cd_type = 'T', sval = 0.0 )
539ICE_aire = OCE_aire
540
541def OCE_stock_int (stock) :
542    '''Integrate stock on ocean grid'''
543    OCE_stock_int = wu.Psum (  (stock * OCE_aire ).to_masked_array().ravel() )
544    return OCE_stock_int
545
546def ONE_stock_int (stock) :
547    '''Sum stock'''
548    ONE_stock_int =  wu.Psum ( stock.to_masked_array().ravel() )
549    return ONE_stock_int
550
551def OCE_flux_int (flux) :
552    '''Integrate flux on oce grid'''
553    OCE_flux_int =  wu.Psum  ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel() )
554    return OCE_flux_int
555
556def ONE_flux_int (flux) :
557    '''Integrate flux on oce grid'''
558    OCE_flux_int =  wu.Psum ( (flux * dtime_per_sec).to_masked_array().ravel() )
559    return OCE_flux_int
560
561OCE_aire_tot = ONE_stock_int ( OCE_aire )
562ICE_aire_tot = ONE_stock_int ( ICE_aire )
563   
564echo ( '\n------------------------------------------------------------------------------------' )
565echo ( '-- NEMO change in stores (for the records)' )
566#-- Note that the total number of days of the run should be diagnosed rather than assumed
567#-- Here the result is in Sv
568#
569#-- Change in ocean volume in freshwater equivalent
570
571OCE_ssh_beg = d_OCE_beg['sshn']
572OCE_ssh_end = d_OCE_end['sshn']
573
574OCE_sum_ssh_beg = OCE_stock_int ( OCE_ssh_beg )
575OCE_sum_ssh_end = OCE_stock_int ( OCE_ssh_end )
576
577if NEMO == 3.6 :
578    OCE_e3tn_beg = d_OCE_beg['fse3t_n']
579    OCE_e3tn_end = d_OCE_end['fse3t_n']
580    OCE_sum_e3tn_beg = OCE_stock_int ( OCE_e3tn_beg * OCE_msk3)
581    OCE_sum_e3tn_end = OCE_stock_int ( OCE_e3tn_end * OCE_msk3)
582
583echo ( f'OCE_sum_ssh_beg = {OCE_sum_ssh_beg:12.6e} m^3  - OCE_sum_ssh_end = {OCE_sum_ssh_end:12.6e} m^3' )
584dOCE_ssh_vol = ( OCE_sum_ssh_end - OCE_sum_ssh_beg )
585dOCE_ssh_mas = dOCE_ssh_vol * OCE_rho_liq
586
587if NEMO == 3.6 :
588    echo ( f'OCE_sum_e3tn_beg = {OCE_sum_e3tn_beg:12.6e} m^3  - OCE_sum_e3tn_end = {OCE_sum_e3tn_end:12.6e} m^3' )
589    dOCE_e3tn_vol = ( OCE_sum_e3tn_end - OCE_sum_e3tn_beg )
590    dOCE_e3tn_mas = dOCE_e3tn_vol * OCE_rho_liq
591
592dOCE_vol_wat = dOCE_ssh_vol ; dOCE_mas_wat = dOCE_ssh_mas
593
594echo ( f'dOCE vol    = {dOCE_vol_wat             :12.3e} m^3' )
595echo ( f'dOCE ssh    = {dOCE_vol_wat/OCE_aire_tot:12.3e} m  ' )
596prtFlux ( 'dOCE mass ', dOCE_mas_wat, 'e' )
597
598if NEMO == 3.6 :
599    echo ( f'dOCE e3tn   vol    = {dOCE_e3tn_vol:12.3e} m^3' )
600    prtFlux ( 'dOCE e3tn   mass', dOCE_e3tn_mas, 'e' )
601
602## Glace et neige
603if NEMO == 3.6 :       
604    ICE_ice_beg = d_ICE_beg['v_i_htc1']+d_ICE_beg['v_i_htc2']+d_ICE_beg['v_i_htc3']+d_ICE_beg['v_i_htc4']+d_ICE_beg['v_i_htc5']
605    ICE_ice_end = d_ICE_end['v_i_htc1']+d_ICE_end['v_i_htc2']+d_ICE_end['v_i_htc3']+d_ICE_end['v_i_htc4']+d_ICE_end['v_i_htc5']
606
607    ICE_sno_beg = d_ICE_beg['v_s_htc1']+d_ICE_beg['v_s_htc2']+d_ICE_beg['v_s_htc3']+d_ICE_beg['v_s_htc4']+d_ICE_beg['v_s_htc5']
608    ICE_sno_end = d_ICE_end['v_s_htc1']+d_ICE_end['v_s_htc2']+d_ICE_end['v_s_htc3']+d_ICE_end['v_s_htc4']+d_ICE_end['v_s_htc5']
609   
610    ICE_pnd_beg = 0.0 ; ICE_pnd_end = 0.0 ; ICE_fzl_beg = 0.0 ; ICE_fzl_end = 0.0
611
612    ICE_mas_wat_beg = ICE_ice_beg*ICE_rho_ice + ICE_sno_beg*ICE_rho_sno
613    ICE_mas_wat_end = ICE_ice_end*ICE_rho_ice + ICE_sno_end*ICE_rho_sno
614   
615if NEMO == 4.0 or NEMO == 4.2 :
616    ICE_ice_beg = d_ICE_beg['v_i']  ; ICE_ice_end = d_ICE_end['v_i']
617    ICE_sno_beg = d_ICE_beg['v_s']  ; ICE_sno_end = d_ICE_end['v_s'] 
618    ICE_pnd_beg = d_ICE_beg['v_ip'] ; ICE_pnd_end = d_ICE_end['v_ip'] # Ice ponds
619    ICE_fzl_beg = d_ICE_beg['v_il'] ; ICE_fzl_end = d_ICE_end['v_il'] # Frozen Ice Ponds
620   
621    ICE_mas_wat_beg =  OCE_stock_int ( d_ICE_beg['snwice_mass'] )
622    ICE_mas_wat_end =  OCE_stock_int ( d_ICE_end['snwice_mass'] )
623   
624ICE_vol_ice_beg = OCE_stock_int ( ICE_ice_beg )
625ICE_vol_ice_end = OCE_stock_int ( ICE_ice_end )
626
627ICE_vol_sno_beg = OCE_stock_int ( ICE_sno_beg )
628ICE_vol_sno_end = OCE_stock_int ( ICE_sno_end )
629
630ICE_vol_pnd_beg = OCE_stock_int ( ICE_pnd_beg )
631ICE_vol_pnd_end = OCE_stock_int ( ICE_pnd_end )
632
633ICE_vol_fzl_beg = OCE_stock_int ( ICE_fzl_beg )
634ICE_vol_fzl_end = OCE_stock_int ( ICE_fzl_end )
635
636#-- Converting to freswater volume
637dICE_vol_ice = ICE_vol_ice_end - ICE_vol_ice_beg
638dICE_mas_ice = dICE_vol_ice * ICE_rho_ice
639
640dICE_vol_sno = ICE_vol_sno_end - ICE_vol_sno_beg
641dICE_mas_sno = dICE_vol_sno * ICE_rho_sno
642
643dICE_vol_pnd = ICE_vol_pnd_end - ICE_vol_pnd_beg
644dICE_mas_pnd = dICE_vol_pnd * ICE_rho_pnd
645
646dICE_vol_fzl= ICE_vol_fzl_end - ICE_vol_fzl_beg
647dICE_mas_fzl = dICE_vol_fzl * ICE_rho_pnd
648
649if NEMO == 3.6 :
650    dICE_mas_wat = dICE_mas_ice + dICE_mas_sno 
651    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_ice + dICE_mas_sno
652
653if NEMO == 4.0 or NEMO == 4.2 :
654    dICE_mas_wat = ICE_mas_wat_end - ICE_mas_wat_beg
655    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat
656
657echo ( f'ICE_vol_ice_beg = {ICE_vol_ice_beg:12.6e} m^3 | ICE_vol_ice_end = {ICE_vol_ice_end:12.6e} m^3' )
658echo ( f'ICE_vol_sno_beg = {ICE_vol_sno_beg:12.6e} m^3 | ICE_vol_sno_end = {ICE_vol_sno_end:12.6e} m^3' )
659echo ( f'ICE_vol_pnd_beg = {ICE_vol_pnd_beg:12.6e} m^3 | ICE_vol_pnd_end = {ICE_vol_pnd_end:12.6e} m^3' )
660echo ( f'ICE_vol_fzl_beg = {ICE_vol_fzl_beg:12.6e} m^3 | ICE_vol_fzl_end = {ICE_vol_fzl_end:12.6e} m^3' )
661
662echo ( f'dICE_vol_ice   = {dICE_vol_ice:12.3e} m^3' )
663echo ( f'dICE_vol_sno   = {dICE_vol_sno:12.3e} m^3' )
664echo ( f'dICE_vol_pnd   = {dICE_vol_pnd:12.3e} m^3' )
665echo ( f'dICE_mas_ice   = {dICE_mas_ice:12.3e} m^3' )
666echo ( f'dICE_mas_sno   = {dICE_mas_sno:12.3e} m^3' ) 
667echo ( f'dICE_mas_pnd   = {dICE_mas_pnd:12.3e} m^3' ) 
668echo ( f'dICE_mas_fzl   = {dICE_mas_fzl:12.3e} m^3' ) 
669
670echo ( '\n------------------------------------------------------------')
671echo ( 'Change in water content (ocean + ice) ' )
672prtFlux ( 'dMass (ocean)', dSEA_mas_wat, 'e', True )
673
674
675###
676### And now, working on fluxes !!
677
678# if coupled:
679# emp_oce = evap - precip (including blowing snow) - calving
680# emp_ice = evap - precip (excluding blowing snow)
681# emp_ice = wfx_spr(<0) + wfx_sub(>0) + wfx_snw_sub(v4.0.6) - wfx_err_sub(<0)
682# runoffs = rivers + icebergs
683# empmr = emp_oce - wfx_ice - wfx_snw - wfx_pnd(v4.0.6) - wfx_err_sub - runoffs
684# doce+ice = - evap + precip + calving (emp_oce)
685#            + rivers + icebergs       (friver+iceberg ou runoffs)
686#            + iceshelf                (iceshelf)
687#            - emp_ice                 (emp_ice)
688# dice = - emp_ice - wfx_snw - wfx_ice - wfx_pnd(v4.0.6) - wfx_err_sub
689
690#OCE_emp = evap - precip (including blowing snow) - calving
691#ICE_emp = wfx_spr(<0) + wfx_sub(>0) + wfx_snw_sub(v4.0.6) - wfx_err_sub(<0)
692
693echo ( '\n------------------------------------------------------------------------------------' )
694echo ( '-- Checks in NEMO - from budget_modipsl.sh (Clément Rousset)' )
695
696# Read variable and computes integral over space and time
697OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
698OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
699OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
700OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
701OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
702OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
703OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
704OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
705OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
706if NEMO == 4.0 or NEMO == 4.2 :
707    OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
708    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
709    OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
710if NEMO == 3.6 :
711    OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
712    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
713    OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
714# Additional checks
715OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
716ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
717OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
718OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
719OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
720ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
721if NEMO == 4.0 or NEMO == 4.2 :
722    ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
723    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
724    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
725if NEMO == 3.6 :
726    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
727    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
728    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
729
730OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
731if OCE_relax  :
732    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
733    OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
734    OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
735else : 
736    OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0
737    OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0
738if OCE_icb :
739    OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
740    OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
741else :
742    OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0.
743    OCE_calv_icb = 0. ; OCE_mas_calv_icb = 0.
744
745OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err
746OCE_mas_all = OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
747
748echo ('\n-- Fields:' ) 
749prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True)
750prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True)
751prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True)
752prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True)
753prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True)
754prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True)
755prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e', True)
756prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True)
757prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True)
758prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e', True) 
759prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True)
760prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True)
761prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True)
762prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True)
763prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True)
764prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True)
765prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True)
766prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True)
767prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e'      )
768prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True)
769prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True)
770prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True)
771prtFlux ('  RAIN         ', OCE_mas_rain      , 'e'      )
772prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True)
773prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True)
774prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True)
775prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True)
776prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
777
778echo (' ')
779if Coupled :
780    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceshelf, 'e', True )
781else :
782    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf, 'e', True )
783
784
785echo     ('\n===>  Comparing ===>' ) 
786echo     ('  WFX OCE                     = -empmr + iceshelf                                            = {:12.5e} (kg) '.format ( -OCE_mas_empmr + OCE_mas_iceshelf) )
787echo     ('     versus d OCE                                                                            = {:12.5e} (kg) '.format ( dOCE_mas_wat) )
788echo     ('  WFX ICE+SNW+PND 1           = emp_ice - wfxice - wfxsnw - wfxpnd - wfxsub_err              = {:12.5e} (kg) '.format ( -OCE_mas_emp_ice - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err) )
789echo     ('  WFX ICE+SNW+PND 2           = -emp_ice + empmr - emp_oce + runoffs                         = {:12.5e} (kg) '.format ( -OCE_mas_emp_ice + OCE_mas_empmr - OCE_mas_emp_oce + OCE_mas_runoffs )) 
790echo     ('     versus d ICE+SNW+PND                                                                    = {:12.5e} (kg) '.format ( dICE_mas_wat))  # Manque PND ?
791if Coupled : 
792    echo ('  WFX OCE+ICE+SNW+PND         = -emp_oce - emp_ice + runoffs + iceshelf                      = {:12.5e} (kg) '.format ( -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceshelf))
793else :
794    echo ('  WFX OCE+ICE+SNW+PND         = -emp_oce - emp_ice + runoffs + iceberg + calving + iceshelf  = {:12.5e} (kg) '.format ( -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf))
795echo     ('     versus d OCE+ICE+SNW+PND                                                                = {:12.5e} (kg) '.format ( dSEA_mas_wat  )) # Manque PND
796
797echo ('\n===> Leaks ===>')
798echo ('   Leak OCE             = dOCE_mas_wat + empmr - iceshelf                                                 = {:12.3e} (kg) '.format ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) ) 
799echo ('                        = (dOCE_mas_wat + empmr - iceshelf)/abs(dOCE_mas_wat)                             = {:12.1e} (-)  '.format ( (dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf )  / abs (dOCE_mas_wat)  ) )
800echo ('   Leak ICE+SNW+PND 1   = dICE_mas_wat + emp_ice + wfxice + wfxsnw + wfxpnd + wfxsub_err                  = {:12.3e} (kg) '.format ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) )
801echo ('                        = (dICE_mas_wat + emp_ice + wfxice + wfxsnw + wfxpnd + wfxsub_err )/dICE_mas_wat  = {:12.1e} (-)  '.format ( (dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err)/abs(dICE_mas_wat)) )
802echo ('   Leak ICE+SNW+PND 2   = dICE_mas_wat + emp_ice - empmr + emp_oce - runoffs                              = {:12.3e} (kg) '.format ( dICE_mas_wat + OCE_mas_emp_ice - OCE_mas_empmr + OCE_mas_emp_oce - OCE_mas_runoffs ))
803echo ('                        = (dICE_mas_wat - empmr  + emp_oce - runoffs)/abs(dICE_mas_wat)                   = {:12.1e} (-)  '.format ( (dICE_mas_wat - OCE_mas_empmr  + OCE_mas_emp_oce - OCE_mas_runoffs)/abs(dICE_mas_wat)) )
804echo ('   Leak OCE+ICE+SNW+PND =  d(ICE+OCE)_mas_wat + emp_oce + emp_ice - runoffs - iceshelf                    = {:12.3e} (kg) '.format ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ))
805echo ('                        = (dSEA_mas_wat + emp_oce + emp_ice - runoffs - iceshelf)*/abs(dSEA_mas_wat)      = {:12.1e} (-)  '.format ( (dSEA_mas_wat + OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf)/abs(dSEA_mas_wat) ) )
806
807prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e', True )
808prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e', True )
809prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e', True )
810
811
812
813
814# check if emp     = emp_ice + emp_oce - calving
815#          emp_ice = wfxsub + wfxsnw_sub + wfxsnw_pre - wfxsub_err
816
817echo     ('\n===> Checks ===>' )
818echo     ('   Check EMPMR           = empmr - emp_oce + runoffs + wfxice + wfxsnw + wfxpnd + wfxsub_err = {:12.5e} (kg) '.format ( OCE_mas_empmr - OCE_mas_emp_oce + OCE_mas_runoffs + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ))
819if Coupled : 
820    echo ('   Check EMP_OCE         = emp_oce + fwb + ssr - evap_oce + rain + snow_oce + calving        = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_fwb + OCE_mas_ssr - OCE_mas_evap_oce + OCE_mas_rain + OCE_mas_snow_oce + OCE_mas_calving ))
821else :
822    echo ('   Check EMP_OCE         = emp_oce + ssr + fwb - evap_oce + rain + snow_oce                  = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_ssr + OCE_mas_fwb - OCE_mas_evap_oce + OCE_mas_rain + OCE_mas_snow_oce ))
823echo     ('   Check EMP_ICE         = emp_ice - evap_ice + snow_ice                                     = {:12.5e} (kg) '.format ( OCE_mas_emp_ice - ICE_mas_evap_ice + OCE_mas_snow_ice ))
824echo     ('   Check EMP_ICE vs WFX  = emp_ice - wfxsub - wfxsnw_sub - wfxsnw_pre + wfxsub_err           = {:12.5e} (kg) '.format ( OCE_mas_emp_ice - OCE_mas_wfxsub - ICE_mas_wfxsnw_sub - ICE_mas_wfxsnw_pre + ICE_mas_wfxsub_err ))
825
826echo ( '\n------------------------------------------------------------------------------------' )
827echo ( '-- Computations in the PDF note of Clément Rousset')
828echo ( 'Freshwater budget of the ice-ocean system          = emp_oce + emp_ice - runoffs - iceberg - iceshelf                = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceberg - OCE_mas_iceshelf ))
829echo ( 'Freshwater budget of the ice-ocean system          = emp_oce + emp_ice - friver  - iceberg - iceshelf                = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_friver  - OCE_mas_iceberg - OCE_mas_iceshelf ))
830echo ( 'Freshwater budget in the ocean = ocean mass change = emp_oce - vfxice - ( vfxsnw + vfxsub_err ) - runoffs - iceshelf = {:12.5e} (kg) '.format ( OCE_mas_emp_oce - OCE_mas_wfxice - ( OCE_mas_wfxsnw + ICE_mas_wfxsub_err ) - OCE_mas_runoffs - OCE_mas_iceshelf ))
831echo ( 'Freshwater budget in the ocean = ocean mass change = emp_oce - vfxice - ( vfxsnw + vfxsub_err ) - friver  - iceshelf = {:12.5e} (kg) '.format ( OCE_mas_emp_oce - OCE_mas_wfxice - ( OCE_mas_wfxsnw + ICE_mas_wfxsub_err ) - OCE_mas_friver  - OCE_mas_iceshelf ))
832echo ( 'Freshwater budget in the ice = ice mass change     = vfxice + ( vfxsnw + vfxsub + vfxspr )                           = {:12.5e} (kg) '.format ( OCE_mas_wfxice  + OCE_mas_wfxsnw + OCE_mas_wfxsub + ICE_mas_wfxsnw_pre ))
833echo ( 'Freshwater flux at the interface [ice/snow]-ocean  = vfxice + ( vfxsnw + vfxsub_err )                                = {:12.5e} (kg) '.format ( OCE_mas_wfxsub  + ICE_mas_wfxsnw_pre ))
834echo ( 'Freshwater flux at the interface [ice/snow]-atm    = ( vfxsub + vfxspr )                                             = {:12.5e} (kg) '.format ( OCE_mas_emp_ice + ICE_mas_wfxsub_err ))
835echo ( 'Freshwater flux at the interface [ice/snow]-atm    = emp_ice + vfxsub_err                                            = {:12.5e} (kg) '.format ( OCE_mas_emp_ice + ICE_mas_wfxsub_err ))
836echo ( 'Freshwater flux at the interface ocean-atm         = emp_oce + calving - vfxsub_err                                  = {:12.5e} (kg) '.format ( OCE_mas_emp_oce + OCE_mas_calving - ICE_mas_wfxsub_err ))
837
838echo ( "scsshtot   : global_average_sea_level_change                            = {:12.3e} (m) ".format ( np.sum (d_OCE_sca['scsshtot']  ).values  ) )
839echo ( "scsshtot   : global_average_sea_level_change                            = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['scsshtot']  ).values * OCE_aire_tot*OCE_rho_liq  ) )
840echo ( "bgvolssh   : drift in global mean ssh volume wrt timestep 1             = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['bgvolssh']  ).values * 1e9 * OCE_rho_liq  ) )
841echo ( "bgvole3t   : drift in global mean volume variation (e3t) wrt timestep 1 = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['bgvole3t']  ).values * 1e9 * OCE_rho_liq  ) )
842echo ( "bgfrcvol   : global mean volume from forcing                            = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['bgfrcvol']  ).values * 1e9 * OCE_rho_liq  ) )
843echo ( "ibgvol_tot : global mean ice volume                                     = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['ibgvol_tot']).values * 1e9 * OCE_rho_liq  ) )
844echo ( "sbgvol_tot : global mean snow volume                                    = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['sbgvol_tot']).values * 1e9 * OCE_rho_liq  ) )
845echo ( "ibgvolume  : drift in ice/snow volume (equivalent ocean volume)         = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['ibgvolume'] ).values * 1e9 * OCE_rho_liq  ) )
846
Note: See TracBrowser for help on using the repository browser.