source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6271

Last change on this file since 6271 was 6271, checked in by omamce, 19 months ago

O.M. : Water Bugget

Various improvments
OCE budget is OK
Use .ini file as input

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