source: TOOLS/WATER_BUDGET/OCE_waterbudget.py @ 6264

Last change on this file since 6264 was 6264, checked in by omamce, 21 months ago

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