source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6669

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

O.M. : WATER BUDGET

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