source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6446

Last change on this file since 6446 was 6277, checked in by omamce, 18 months ago

O.M. : WATER_BUDGET

Switch to .ini file to read configuratio
Correct SECHIBA budget
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 38.8 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##
14## SVN information
15#  $Author$
16#  $Date$
17#  $Revision$
18#  $Id$
19#  $HeadURL$
20
21###
22## Import system modules
23import sys, os, shutil, subprocess, platform
24import numpy as np
25import configparser, re
26
27## Creates parser
28config = configparser.ConfigParser()
29config.optionxform = str # To keep capitals
30
31## Read experiment parameters
32ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ;  TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
33
34##
35ARCHIVE=None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
36TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None
37file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None 
38file_restart_beg=None ; file_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
39file_RUN_beg=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
40
41# Arguments passed
42print ( "Name of Python script:", sys.argv[0] )
43IniFile = sys.argv[1]
44print ("Input file : ", IniFile )
45config.read (IniFile)
46if 'full' in IniFile : FullIniFile = IniFile
47else                 : FullIniFile = 'full_' + IniFile
48
49def setBool (chars) :
50    '''Convert specific char string in boolean if possible'''
51    setBool = chars
52    for key in configparser.ConfigParser.BOOLEAN_STATES.keys () :
53        if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key]
54    return setBool
55
56def setNum (chars) :
57    '''Convert specific char string in integer or real if possible'''
58    if type (chars) == str :
59        realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$")
60        isReal = realnum.match(chars.strip()) != None
61        isInt  = chars.strip().isdigit()
62        if isReal :
63            if isInt : setNum = int   (chars)
64            else     : setNum = float (chars)
65        else : setNum = chars
66    else : setNum = chars
67    return setNum
68
69def setNone (chars) :
70    '''Convert specific char string to None if possible'''
71    if type (chars) == str :
72        if chars in ['None', 'NONE', 'none'] : setNone = None
73        else : setNone = chars
74    else : setNone = chars
75    return setNone
76
77## Reading config
78for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] :
79    if Section in config.keys() : 
80        print ( f'[{Section}]' )
81        for VarName in config[Section].keys() :
82            locals()[VarName] = config[Section][VarName]
83            locals()[VarName] = setBool (locals()[VarName])
84            locals()[VarName] = setNum  (locals()[VarName])
85            print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) )
86
87if not 'Files' in config.keys() : config['Files'] = {}
88
89def unDefined (char) :
90    if char in globals () :
91        if char == None : return True
92        else : return False
93    else : return True
94
95##-- Some physical constants
96#-- Earth Radius
97if not 'Ra'            in locals () : Ra = 6366197.7236758135
98#-- Gravity
99if not 'Grav'          in locals () : Grav = 9.81
100#-- Ice volumic mass (kg/m3) in LIM3
101if not 'ICE_rho_ice'   in locals () : ICE_rho_ice = 917.0
102#-- Snow volumic mass (kg/m3) in LIM3
103if not 'ICE_rho_sno'   in locals () : ICE_rho_sno = 330.0
104    #-- Water density in ice pounds in SI3
105if not 'ICE_rho_pnd'   in locals () : ICE_rho_pnd = 1000.
106#-- Ocean water volumic mass (kg/m3) in NEMO
107if not 'OCE_rho_liq'   in locals () : OCE_rho_liq = 1026.
108#-- Water volumic mass in atmosphere
109if not 'ATM_rho'       in locals () : ATM_rho = 1000.
110#-- Water volumic mass in surface reservoirs
111if not 'SRF_rho'       in locals () : SRF_rho = 1000.
112#-- Water volumic mass of rivers
113if not 'RUN_rho'       in locals () : RUN_rho = 1000.
114#-- Year length
115if not 'YearLength'    in locals () : YearLength = 365.25 * 24. * 60. * 60.
116
117config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ICE_rho_pnd':ICE_rho_pnd,
118                          'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
119       
120
121
122# Where do we run ?
123SysName, NodeName, Release, Version, Machine = os.uname()
124TGCC  = ( 'irene'   in NodeName )
125IDRIS = ( 'jeanzay' in NodeName )
126
127## Set site specific libIGCM directories, and other specific stuff
128if TGCC :
129    CPU = subprocess.getoutput ( 'lscpu | grep "Model name"' )
130    if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene'
131    if "AMD"                       in CPU : Machine = 'irene-amd'
132       
133    ARCHIVE     = subprocess.getoutput ( f'ccc_home --cccstore   -d {Project} -u {User}' )
134    STORAGE     = subprocess.getoutput ( f'ccc_home --cccwork    -d {Project} -u {User}' )
135    SCRATCHDIR  = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' )
136    R_IN        = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg'), 'IGCM')
137    rebuild     = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg'), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' )
138
139    ## Specific to run at TGCC.
140    # Needed before importing a NetCDF library (netCDF4, xarray, cmds, etc ...)
141    import mpi4py
142    mpi4py.rc.initialize = False
143       
144    ## Creates output directory
145    TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
146
147if IDRIS :
148    raise Exception ("Pour IDRIS : repertoires et chemins a definir") 
149
150config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU,
151                    'Program'  : "Generated by : " + str(sys.argv), 
152                    'HOSTNAME' : platform.node (), 'LOGNAME'  : os.getlogin (),
153                    'Python'   : f'{platform.python_version ()}',
154                    'OS'       : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}',
155                    'SVN_Author'   : "$Author$", 
156                    'SVN_Date'     : "$Date$",
157                    'SVN_Revision' : "$Revision$",
158                    'SVN_Id'       : "$Id$",
159                    'SVN_HeadURL'  : "$HeadURL$"}
160
161if libIGCM :
162    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
163
164## Import specific module
165import nemo
166## Now import needed scientific modules
167import xarray as xr
168   
169# Output file
170if FileOut == None : 
171    FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
172    config['Files']['FileOut'] = FileOut
173
174f_out = open ( FileOut, mode = 'w' )
175   
176# Function to print to stdout *and* output file
177def echo (string, end='\n') :
178    print ( str(string), end=end  )
179    sys.stdout.flush ()
180    f_out.write ( str(string) + end )
181    f_out.flush ()
182    return None
183   
184## Set libIGCM directories
185R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT')
186R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT')
187
188L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName)
189R_SAVE      = os.path.join ( R_OUT, L_EXP )
190R_BUFR      = os.path.join ( R_BUF, L_EXP )
191POST_DIR    = os.path.join ( R_BUF, L_EXP, 'Out' )
192REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' )
193R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
194R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
195   
196#if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir )
197if not os.path.isdir (TmpDir) : os.mkdir (TmpDir)
198TmpDirOCE = os.path.join (TmpDir, 'OCE')
199TmpDirICE = os.path.join (TmpDir, 'ICE')
200if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE )
201if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE )
202
203echo (' ')
204echo ( f'JobName : {JobName}' )
205echo (Comment)
206echo ( f'Working in TMPDIR : {TmpDir}' )
207
208echo ( f'\nDealing with {L_EXP}' )
209
210#-- Model output directories
211if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' )
212if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' )
213if dir_OCE_his == None :
214    dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir )
215    config['Files']['dir_OCE_his'] = dir_OCE_his
216if dir_ICE_his == None :
217    dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir )
218    config['Files']['dir_OCE_his'] = dir_OCE_his
219 
220echo ( f'The analysis relies on files from the following model output directories : ' )
221echo ( f'{dir_OCE_his}' )
222echo ( f'{dir_ICE_his}' )
223
224#-- Files Names
225if Period == None :
226    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
227    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
228    config['Files']['Period'] = Period
229if FileCommon == None :
230    FileCommon = f'{JobName}_{Period}'
231    config['Files']['FileCommon'] = FileCommon
232
233if Title == None :
234    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
235    config['Files']['Title'] = Title
236       
237       
238echo ('\nOpen history files' )
239if file_OCE_his == None :
240    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' )
241    file_OCE_his = file_OCE_his
242if file_OCE_sca == None :   
243    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' )
244    config['Files']['file_OCE_sca'] = file_OCE_sca
245if file_ICE_his == None : 
246    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' )
247    config['Files']['file_ICE_his'] = file_ICE_his
248
249d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
250d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
251d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
252if NEMO == 3.6 :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
253
254echo ( file_OCE_his )
255echo ( file_OCE_sca )
256echo ( file_ICE_his )
257
258## Compute run length
259dtime = ( d_OCE_his.time_counter_bounds.max() - d_OCE_his.time_counter_bounds.min() )
260echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
261dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
262
263## Compute length of each period
264dtime_per = ( d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] )
265echo ('\nPeriods lengths (days) : ')
266echo (' {:}'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
267dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
268dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_OCE_his.time_counter,] )
269dtime_per_sec.attrs['unit'] = 's'
270
271## Number of years
272NbYear = dtime_sec / YearLength
273#-- Open restart files
274YearRes = YearBegin - 1              # Year of the restart of beginning of simulation
275YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
276
277echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
278
279if TarRestartPeriod_beg == None : 
280    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
281    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231'
282    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
283
284if TarRestartPeriod_end == None : 
285    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
286    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
287    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231'
288    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
289
290if file_restart_beg == None :
291    file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' )
292    config['Files']['file_restart_beg'] = file_restart_beg
293if file_restart_end == None :
294    file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' )
295    config['Files']['file_restart_end'] = file_restart_end
296
297echo ( f'{file_restart_beg}' )
298echo ( f'{file_restart_end}' )
299
300if file_OCE_beg == None : 
301    file_OCE_beg = f'{TmpDir}/OCE_{JobName}_{YearRes}1231_restart.nc'
302    config['Files']['file_OCE_beg'] = file_OCE_beg
303if file_OCE_end == None :
304    file_OCE_end = f'{TmpDir}/OCE_{JobName}_{YearEnd}1231_restart.nc'
305    config['Files']['file_OCE_end'] = file_OCE_end
306if file_ICE_beg == None :
307    file_ICE_beg = f'{TmpDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc'
308    config['Files']['file_ICE_beg'] = file_ICE_beg
309if file_ICE_end == None : 
310    file_ICE_end = f'{TmpDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc'
311    config['Files']['file_ICE_end'] = file_ICE_end
312
313echo ( f'{file_OCE_beg}' )
314echo ( f'{file_OCE_end}' )
315echo ( f'{file_ICE_beg}' )
316echo ( f'{file_ICE_end}' )
317
318echo ('\nExtract and rebuild OCE and ICE restarts')
319def get_ndomain (zfile) :
320    #d_zfile = xr.open_dataset (zfile, decode_times=False).squeeze()
321    #ndomain_opa = d_zfile.attrs['DOMAIN_number_total']
322    #d_zfile.close ()
323    ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ).format()
324    return int (ndomain_opa)
325
326if not os.path.exists ( file_OCE_beg) :
327    echo ( f'Extracting {file_OCE_beg}' )
328    base_file_OCE_beg = os.path.basename (file_OCE_beg)
329    if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart_0000.nc') ) :
330        command =  f'cd {TmpDirOCE} ; tar xf {file_restart_beg}  OCE_{JobName}_{YearRes}1231_restart_*.nc'
331        echo ( command )
332        os.system ( command )
333    echo ('extract ndomain' )
334    ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') )
335    command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {base_file_OCE_beg} {TmpDir}'
336    echo ( command )
337    os.system ( command )
338    echo ( f'Rebuild done : {file_OCE_beg}')
339   
340if not os.path.exists ( file_OCE_end) :
341    echo ( f'Extracting {file_OCE_end}' )
342    base_file_OCE_end = os.path.basename (file_OCE_end)
343    if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart_0000.nc') ):
344        command =  f'cd {TmpDirOCE} ; tar xf {file_restart_end}  OCE_{JobName}_{YearEnd}1231_restart_*.nc'
345        echo ( command )
346        os.system ( command )
347    echo ('extract ndomain' )
348    ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart') )
349    command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {base_file_OCE_end} {TmpDir}'
350    echo ( command )
351    os.system ( command )
352    echo ( f'Rebuild done : {file_OCE_end}')
353
354if not os.path.exists ( file_ICE_beg) :
355    echo ( f'Extracting {file_ICE_beg}' )
356    base_file_ICE_beg = os.path.basename (file_ICE_beg)
357    if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod_0000.nc') ):
358        command =  f'cd {TmpDirICE} ; tar xf {file_restart_beg}  ICE_{JobName}_{YearRes}1231_restart_icemod_*.nc'
359        echo ( command )
360        os.system ( command )
361    echo ('extract ndomain' )
362    ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') )
363    command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_beg} {TmpDir} '
364    echo ( command )
365    os.system ( command )
366    echo ( f'Rebuild done : {file_OCE_beg}')
367   
368if not os.path.exists ( file_ICE_end ) :
369    echo ( f'Extracting {file_ICE_end}' )
370    base_file_ICE_end = os.path.basename (file_ICE_end)
371    if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearEnd}1231_restart_icemod_0000.nc') ):
372        command =  f'cd {TmpDirICE} ; tar xf {file_restart_end} ICE_{JobName}_{YearEnd}1231_restart_icemod_*.nc'
373        echo ( command )
374        os.system ( command )
375    echo ('extract ndomain' )
376    ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') )
377    command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_end} {TmpDir}'
378    echo ( command )
379    os.system ( command )
380    echo ( f'Rebuild done : {file_ICE_end}')
381
382echo ('Opening OCE and ICE restart files')
383if NEMO == 3.6 : 
384    d_OCE_beg = xr.open_dataset ( os.path.join (TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
385    d_OCE_end = xr.open_dataset ( os.path.join (TmpDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()
386    d_ICE_beg = xr.open_dataset ( os.path.join (TmpDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()
387    d_ICE_end = xr.open_dataset ( os.path.join (TmpDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()
388if NEMO == 4.0 or NEMO == 4.2 : 
389    d_OCE_beg = xr.open_dataset ( os.path.join (TmpDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
390    d_OCE_end = xr.open_dataset ( os.path.join (TmpDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
391    d_ICE_beg = xr.open_dataset ( os.path.join (TmpDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
392    d_ICE_end = xr.open_dataset ( os.path.join (TmpDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
393
394##
395config_out = open (FullIniFile, 'w')
396config.write (config_out )
397config_out.close ()
398
399   
400def kg2Sv  (val, rho=OCE_rho_liq) :
401    '''From kg to Sverdrup'''
402    return val/dtime_sec*1.0e-6/rho
403
404def kg2myear (val, rho=OCE_rho_liq) :
405    '''From kg to m/year'''
406    return val/OCE_aire_tot/rho/NbYear
407
408def var2prt (var, small=False) :
409    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000.
410    else     :  return var , kg2Sv(var)      , kg2myear(var)
411
412def prtFlux (Desc, var, Form='F', small=False) :
413    if small :
414        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
415        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
416    else :
417        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
418        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
419    echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) )
420    return None
421
422# Get mask and surfaces
423sos = d_OCE_his ['sos'][0].squeeze()
424OCE_msk = nemo.lbc_mask ( xr.where ( sos>0., 1., 0. ), cd_type = 'T', sval = 0. )
425
426so = sos = d_OCE_his ['sos'][0].squeeze()
427OCE_msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', sval = 0. )
428
429# lbc_mask removes the duplicate points (periodicity and north fold)
430OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE_msk, cd_type = 'T', sval = 0.0 )
431ICE_aire = OCE_aire
432
433def OCE_stock_int (stock) :
434    '''Integrate stock on ocean grid'''
435    OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) )
436    return OCE_stock_int
437
438def ONE_stock_int (stock) :
439    '''Sum stock'''
440    ONE_stock_int = np.sum (  np.sort ( (stock ).to_masked_array().ravel()) )
441    return ONE_stock_int
442
443def OCE_flux_int (flux) :
444    '''Integrate flux on oce grid'''
445    OCE_flux_int = np.sum (  np.sort ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel()) )
446    return OCE_flux_int
447
448def ONE_flux_int (flux) :
449    '''Integrate flux on oce grid'''
450    OCE_flux_int = np.sum (  np.sort ( (flux * dtime_per_sec).to_masked_array().ravel()) )
451    return OCE_flux_int
452
453OCE_aire_tot = ONE_stock_int ( OCE_aire )
454ICE_aire_tot = ONE_stock_int ( ICE_aire )
455   
456echo ( '\n------------------------------------------------------------------------------------' )
457echo ( '-- NEMO change in stores (for the records)' )
458#-- Note that the total number of days of the run should be diagnosed rather than assumed
459#-- Here the result is in Sv
460#
461#-- Change in ocean volume in freshwater equivalent
462
463OCE_ssh_beg = d_OCE_beg['sshn']
464OCE_ssh_end = d_OCE_end['sshn']
465
466OCE_sum_ssh_beg = OCE_stock_int ( OCE_ssh_beg )
467OCE_sum_ssh_end = OCE_stock_int ( OCE_ssh_end )
468
469if NEMO == 3.6 :
470    OCE_e3tn_beg = d_OCE_beg['fse3t_n']
471    OCE_e3tn_end = d_OCE_end['fse3t_n']
472    OCE_sum_e3tn_beg = OCE_stock_int ( OCE_e3tn_beg * OCE_msk3)
473    OCE_sum_e3tn_end = OCE_stock_int ( OCE_e3tn_end * OCE_msk3)
474
475echo ( 'OCE_sum_ssh_beg = {:12.6e} m^3  - OCE_sum_ssh_end = {:12.6e} m^3'.format (OCE_sum_ssh_beg, OCE_sum_ssh_end) )
476dOCE_ssh_vol = ( OCE_sum_ssh_end - OCE_sum_ssh_beg )
477dOCE_ssh_mas = dOCE_ssh_vol * OCE_rho_liq
478
479if NEMO == 3.6 :
480    echo ( 'OCE_sum_e3tn_beg = {:12.6e} m^3  - OCE_sum_e3tn_end = {:12.6e} m^3'.format (OCE_sum_e3tn_beg, OCE_sum_e3tn_end) )
481    dOCE_e3tn_vol = ( OCE_sum_e3tn_end - OCE_sum_e3tn_beg )
482    dOCE_e3tn_mas = dOCE_e3tn_vol * OCE_rho_liq
483
484dOCE_vol_wat = dOCE_ssh_vol ; dOCE_mas_wat = dOCE_ssh_mas
485
486echo ( 'dOCE vol    = {:12.3e} m^3'.format ( dOCE_vol_wat) )
487echo ( 'dOCE ssh    = {:12.3e} m  '.format ( dOCE_vol_wat/OCE_aire_tot) )
488prtFlux ( 'dOCE mass ', dOCE_mas_wat, 'e' )
489
490if NEMO == 3.6 :
491    echo ( 'dOCE e3tn   vol    = {:12.3e} m^3'.format ( dOCE_e3tn_vol) )
492    prtFlux ( 'dOCE e3tn   mass', dOCE_e3tn_mas, 'e' )
493
494## Glace et neige
495if NEMO == 3.6 :       
496    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']
497    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']
498
499    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']
500    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']
501   
502    ICE_pnd_beg = 0.0 ; ICE_pnd_end = 0.0 ; ICE_fzl_beg = 0.0 ; ICE_fzl_end = 0.0
503
504    ICE_mas_wat_beg = ICE_ice_beg*ICE_rho_ice + ICE_sno_beg*ICE_rho_sno
505    ICE_mas_wat_end = ICE_ice_end*ICE_rho_ice + ICE_sno_end*ICE_rho_sno
506   
507if NEMO == 4.0 or NEMO == 4.2 :
508    ICE_ice_beg = d_ICE_beg['v_i']  ; ICE_ice_end = d_ICE_end['v_i']
509    ICE_sno_beg = d_ICE_beg['v_s']  ; ICE_sno_end = d_ICE_end['v_s'] 
510    ICE_pnd_beg = d_ICE_beg['v_ip'] ; ICE_pnd_end = d_ICE_end['v_ip'] # Ice ponds
511    ICE_fzl_beg = d_ICE_beg['v_il'] ; ICE_fzl_end = d_ICE_end['v_il'] # Frozen Ice Ponds
512   
513    ICE_mas_wat_beg =  OCE_stock_int ( d_ICE_beg['snwice_mass'] )
514    ICE_mas_wat_end =  OCE_stock_int ( d_ICE_end['snwice_mass'] )
515   
516ICE_vol_ice_beg = OCE_stock_int ( ICE_ice_beg )
517ICE_vol_ice_end = OCE_stock_int ( ICE_ice_end )
518
519ICE_vol_sno_beg = OCE_stock_int ( ICE_sno_beg )
520ICE_vol_sno_end = OCE_stock_int ( ICE_sno_end )
521
522ICE_vol_pnd_beg = OCE_stock_int ( ICE_pnd_beg )
523ICE_vol_pnd_end = OCE_stock_int ( ICE_pnd_end )
524
525ICE_vol_fzl_beg = OCE_stock_int ( ICE_fzl_beg )
526ICE_vol_fzl_end = OCE_stock_int ( ICE_fzl_end )
527
528#-- Converting to freswater volume
529dICE_vol_ice = ICE_vol_ice_end - ICE_vol_ice_beg
530dICE_mas_ice = dICE_vol_ice * ICE_rho_ice
531
532dICE_vol_sno = ICE_vol_sno_end - ICE_vol_sno_beg
533dICE_mas_sno = dICE_vol_sno * ICE_rho_sno
534
535dICE_vol_pnd = ICE_vol_pnd_end - ICE_vol_pnd_beg
536dICE_mas_pnd = dICE_vol_pnd * ICE_rho_pnd
537
538dICE_vol_fzl= ICE_vol_fzl_end - ICE_vol_fzl_beg
539dICE_mas_fzl = dICE_vol_fzl * ICE_rho_pnd
540
541if NEMO == 3.6 :
542    dICE_mas_wat = dICE_mas_ice + dICE_mas_sno 
543    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_ice + dICE_mas_sno
544
545if NEMO == 4.0 or NEMO == 4.2 :
546    dICE_mas_wat = ICE_mas_wat_end - ICE_mas_wat_beg
547    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat
548
549echo ( 'ICE_vol_ice_beg = {:12.6e} m^3 | ICE_vol_ice_end = {:12.6e} m^3'.format (ICE_vol_ice_beg, ICE_vol_ice_end) )
550echo ( 'ICE_vol_sno_beg = {:12.6e} m^3 | ICE_vol_sno_end = {:12.6e} m^3'.format (ICE_vol_sno_beg, ICE_vol_sno_end) )
551echo ( 'ICE_vol_pnd_beg = {:12.6e} m^3 | ICE_vol_pnd_end = {:12.6e} m^3'.format (ICE_vol_pnd_beg, ICE_vol_pnd_end) )
552echo ( 'ICE_vol_fzl_beg = {:12.6e} m^3 | ICE_vol_fzl_end = {:12.6e} m^3'.format (ICE_vol_fzl_beg, ICE_vol_fzl_end) )
553
554echo ( 'dICE_vol_ice   = {:12.3e} m^3'.format (dICE_vol_ice) )
555echo ( 'dICE_vol_sno   = {:12.3e} m^3'.format (dICE_vol_sno) )
556echo ( 'dICE_vol_pnd   = {:12.3e} m^3'.format (dICE_vol_pnd) )
557echo ( 'dICE_mas_ice   = {:12.3e} m^3'.format (dICE_mas_ice) )
558echo ( 'dICE_mas_sno   = {:12.3e} m^3'.format (dICE_mas_sno) ) 
559echo ( 'dICE_mas_pnd   = {:12.3e} m^3'.format (dICE_mas_pnd) ) 
560echo ( 'dICE_mas_fzl   = {:12.3e} m^3'.format (dICE_mas_fzl) ) 
561
562echo ( '\n------------------------------------------------------------')
563echo ( 'Variation du contenu en eau ocean + glace ' )
564prtFlux ( 'dMass (ocean)', dSEA_mas_wat, 'e', True )
565
566
567###
568### And now, working on fluxes !!
569
570# if coupled:
571# emp_oce = evap - precip (including blowing snow) - calving
572# emp_ice = evap - precip (excluding blowing snow)
573# emp_ice = wfx_spr(<0) + wfx_sub(>0) + wfx_snw_sub(v4.0.6) - wfx_err_sub(<0)
574# runoffs = rivers + icebergs
575# empmr = emp_oce - wfx_ice - wfx_snw - wfx_pnd(v4.0.6) - wfx_err_sub - runoffs
576# doce+ice = - evap + precip + calving (emp_oce)
577#            + rivers + icebergs       (friver+iceberg ou runoffs)
578#            + iceshelf                (iceshelf)
579#            - emp_ice                 (emp_ice)
580# dice = - emp_ice - wfx_snw - wfx_ice - wfx_pnd(v4.0.6) - wfx_err_sub
581
582#OCE_emp = evap - precip (including blowing snow) - calving
583#ICE_emp = wfx_spr(<0) + wfx_sub(>0) + wfx_snw_sub(v4.0.6) - wfx_err_sub(<0)
584
585echo ( '\n------------------------------------------------------------------------------------' )
586echo ( '-- Checks in NEMO - from budget_modipsl.sh (Clement Rousset)' )
587
588# Read variable and computes integral over space and time
589OCE_empmr      = d_OCE_his['wfo']         ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
590OCE_wfob       = d_OCE_his['wfob']        ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
591OCE_emp_oce    = d_OCE_his['emp_oce']     ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
592OCE_emp_ice    = d_OCE_his['emp_ice']     ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
593OCE_iceshelf   = d_OCE_his['iceshelf']    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
594OCE_calving    = d_OCE_his['calving']     ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
595OCE_iceberg    = d_OCE_his['iceberg']     ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
596OCE_friver     = d_OCE_his['friver']      ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
597OCE_runoffs    = d_OCE_his['runoffs']     ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
598if NEMO == 4.0 or NEMO == 4.2 :
599    OCE_wfxice     = d_OCE_his['vfxice'] ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
600    OCE_wfxsnw     = d_OCE_his['vfxsnw'] ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
601    OCE_wfxsub     = d_OCE_his['vfxsub'] ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
602if NEMO == 3.6 :
603    OCE_wfxice     = d_OCE_his['vfxice']/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
604    OCE_wfxsnw     = d_OCE_his['vfxsnw']/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
605    OCE_wfxsub     = d_OCE_his['vfxsub']/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
606# Additional checks
607OCE_evap_oce   = d_OCE_his['evap_ao_cea'] ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
608ICE_evap_ice   = d_OCE_his['subl_ai_cea'] ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
609OCE_snow_oce   = d_OCE_his['snow_ao_cea'] ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
610OCE_snow_ice   = d_OCE_his['snow_ai_cea'] ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
611OCE_rain       = d_OCE_his['rain']        ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
612ICE_wfxsub_err = d_ICE_his['vfxsub_err']  ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
613if NEMO == 4.0 or NEMO == 4.2 :
614    ICE_wfxpnd     = d_ICE_his['vfxpnd']     ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
615    ICE_wfxsnw_sub = d_ICE_his['vfxsnw_sub'] ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
616    ICE_wfxsnw_pre = d_ICE_his['vfxsnw_pre'] ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
617if NEMO == 3.6 :
618    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
619    ICE_wfxsnw_sub = d_ICE_his['vfxsub']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
620    ICE_wfxsnw_pre = d_ICE_his['vfxspr']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
621
622OCE_wfcorr    = d_OCE_his['wfcorr'] ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
623if OCE_relax  :
624    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
625    OCE_vflx_fwb = d_OCE_his['vflx_fwb'] ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
626    OCE_vflx_ssr = d_OCE_his['vflx_ssr'] ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
627else : 
628    OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0
629    OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0
630if OCE_icb :
631    OCE_berg_icb    = d_OCE_his['berg_floating_melt'] ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
632    OCE_calving_icb = d_OCE_his['calving_icb']        ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
633else :
634    OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0.
635    OCE_calv_icb = 0. ; OCE_mas_calv_icb = 0.
636
637OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err
638
639OCE_mas_all = OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
640
641echo ('\n-- Fields:' ) 
642prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True)
643prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True)
644prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True)
645prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True)
646prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True)
647prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True)
648prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e', True)
649prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True)
650prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True)
651prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e', True) 
652prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True)
653prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True)
654prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True)
655prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True)
656prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True)
657prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True)
658prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True)
659prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True)
660prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e' )
661prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True)
662prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True)
663prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True)
664prtFlux ('  RAIN         ', OCE_mas_rain      , 'e' )
665prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True)
666prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True)
667prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True)
668prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True)
669prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
670
671echo     ('\n===>  Comparing ===>' ) 
672echo     ('  WFX OCE                     = -empmr + iceshelf                                            = {:12.5e} (kg) '.format ( -OCE_mas_empmr + OCE_mas_iceshelf) )
673echo     ('     versus d OCE                                                                            = {:12.5e} (kg) '.format ( dOCE_mas_wat) )
674echo     ('  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) )
675echo     ('  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 )) 
676echo     ('     versus d ICE+SNW+PND                                                                    = {:12.5e} (kg) '.format ( dICE_mas_wat))  # Manque PND ?
677if Coupled : 
678    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))
679else :
680    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))
681echo     ('     versus d OCE+ICE+SNW+PND                                                                = {:12.5e} (kg) '.format ( dSEA_mas_wat  )) # Manque PND
682
683echo ('\n===> Leaks ===>')
684echo ('   Leak OCE             = dOCE_mas_wat + empmr - iceshelf                                                 = {:12.3e} (kg) '.format ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) ) 
685echo ('                        = (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)  ) )
686echo ('   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 ) )
687echo ('                        = (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)) )
688echo ('   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 ))
689echo ('                        = (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)) )
690echo ('   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 ))
691echo ('                        = (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) ) )
692
693
694# check if emp     = emp_ice + emp_oce - calving
695#          emp_ice = wfxsub + wfxsnw_sub + wfxsnw_pre - wfxsub_err
696
697echo     ('\n===> Checks ===>' )
698echo     ('   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 ))
699if Coupled : 
700    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 ))
701else :
702    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 ))
703echo     ('   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 ))
704echo     ('   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 ))
705
706echo ( '\n------------------------------------------------------------------------------------' )
707echo ( '-- Calculs dans la note PDF de Clement')
708echo ( '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 ))
709echo ( '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 ))
710echo ( '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 ))
711echo ( '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 ))
712echo ( '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 ))
713echo ( 'Freshwater flux at the interface [ice/snow]-ocean  = vfxice + ( vfxsnw + vfxsub_err )                                = {:12.5e} (kg) '.format ( OCE_mas_wfxsub  + ICE_mas_wfxsnw_pre ))
714echo ( 'Freshwater flux at the interface [ice/snow]-atm    = ( vfxsub + vfxspr )                                             = {:12.5e} (kg) '.format ( OCE_mas_emp_ice + ICE_mas_wfxsub_err ))
715echo ( 'Freshwater flux at the interface [ice/snow]-atm    = emp_ice + vfxsub_err                                            = {:12.5e} (kg) '.format ( OCE_mas_emp_ice + ICE_mas_wfxsub_err ))
716echo ( '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 ))
717
718echo ( "scsshtot   : global_average_sea_level_change                            = {:12.3e} (m) ".format ( np.sum (d_OCE_sca['scsshtot']  ).values  ) )
719echo ( "scsshtot   : global_average_sea_level_change                            = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['scsshtot']  ).values * OCE_aire_tot*OCE_rho_liq  ) )
720echo ( "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  ) )
721echo ( "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  ) )
722echo ( "bgfrcvol   : global mean volume from forcing                            = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['bgfrcvol']  ).values * 1e9 * OCE_rho_liq  ) )
723echo ( "ibgvol_tot : global mean ice volume                                     = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['ibgvol_tot']).values * 1e9 * OCE_rho_liq  ) )
724echo ( "sbgvol_tot : global mean snow volume                                    = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['sbgvol_tot']).values * 1e9 * OCE_rho_liq  ) )
725echo ( "ibgvolume  : drift in ice/snow volume (equivalent ocean volume)         = {:12.3e} (kg)".format ( np.sum (d_OCE_sca['ibgvolume'] ).values * 1e9 * OCE_rho_liq  ) )
726
Note: See TracBrowser for help on using the repository browser.