source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6265

Last change on this file since 6265 was 6265, checked in by omamce, 20 months ago

O.M. : WATER_BUDGET :

Add nemo and lmdz libraries.
Small changes in routines : ATM water budget still not correct
OCE water budget matches fluxes for NEMO 4.2

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