source: TOOLS/WATER_BUDGET/ATM_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: 18.7 KB
Line 
1#!/usr/bin/env python3
2###
3### Script to check water conservation in the IPSL coupled model
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
23if False : 
24    JobName="TEST-CM72-SIMPLE-ROUTING.12" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6" 
25    Freq = 'MO' ; YearBegin = 1970 ; YearEnd = 1979 ; PackFrequency = 10
26    ATM = 'ICO40' ; Routing = 'SIMPLE'
27
28if False :
29    JobName="TEST-CM72-SIMPLE-ROUTING.10" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6" 
30    Freq = 'MO' ; YearBegin = 1860 ; YearEnd = 1869 ; PackFrequency = 10
31    ATM = 'ICO40' ; Routing = 'SIMPLE'
32
33if False : 
34    JobName="VALID-CM622-LR.01" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6" 
35    Freq = 'MO' ; YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10
36    ATM = 'LMD144142' ; Routing = 'ORCHIDEE'
37   
38if True :
39    JobName="CM65v420-LR-SKL-pi-05" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p48ethe" ; Project="gencmip6" 
40    Freq = 'MO' ;  YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10
41    ORCA = 'eORCA1.2'  ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=4.2 ; OCE_relax = False ; OCE_icb = False ; Coupled = True
42   
43Coupled = True
44
45import numpy as np
46##-- Some physical constants
47#-- Earth Radius
48Ra = 40.e6 / (2. * np.pi)
49#-- Gravity
50g = 9.81
51#-- Ice density (kg/m3) in LIM3
52ICE_rho_ice = 917.0
53#-- Snow density (kg/m3) in LIM3
54ICE_rho_sno = 330.0
55#-- Ocean water density (kg/m3) in NEMO
56OCE_rho_liq = 1026.
57
58###
59ICO  = False
60if 'ICO' in ATM : ICO  = True
61LMDZ = False
62if 'LMD' in ATM : LMDZ = True
63
64###
65## Import system modules
66import sys, os, shutil, subprocess
67
68# Where do we run ?
69TGCC = False ; IDRIS = False
70SysName, NodeName, Release, Version, Machine = os.uname()
71if 'irene'   in NodeName : TGCC  = True
72if 'jeanzay' in NodeName : IDRIS = True
73
74## Set site specific libIGCM directories, and other specific stuff
75if TGCC :
76    CPU = subprocess.getoutput ( 'lscpu | grep "Model name"' )
77    if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene'
78    if "AMD"                       in CPU : Machine = 'irene-amd'
79       
80    ARCHIVE     = subprocess.getoutput ( f'ccc_home --cccstore   -d {Project} -u {User}' )
81    STORAGE     = subprocess.getoutput ( f'ccc_home --cccwork    -d {Project} -u {User}' )
82    SCRATCHDIR  = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' )
83    R_IN        = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM')
84    rebuild     = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' )
85
86    ## Specific to run at TGCC.
87    # Needed before importing a NetCDF library (netCDF4, xarray, cmds, etc ...)
88    import mpi4py
89    mpi4py.rc.initialize = False
90       
91    ## Creates output directory
92    TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
93
94if IDRIS :
95    raise Exception("Pour IDRIS : repertoires et chemins a definir") 
96
97## Import specific module
98import nemo, lmdz
99## Now import needed scientific modules
100import xarray as xr
101   
102# Output file
103FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
104f_out = open ( FileOut, mode = 'w' )
105
106# Function to print to stdout *and* output file
107def echo (string) :
108    print ( string  )
109    f_out.write ( string + '\n' )
110    f_out.flush ()
111    return None
112   
113## Set libIGCM directories
114R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT')
115R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT')
116
117L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName)
118R_SAVE      = os.path.join ( R_OUT, L_EXP )
119R_BUFR      = os.path.join ( R_BUF, L_EXP )
120POST_DIR    = os.path.join ( R_BUF, L_EXP, 'Out' )
121REBUILD_DIR = os.path.join ( R_BUF, L_EXP, 'REBUILD' )
122R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
123R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
124
125#if os.path.isdir (TmpDir) : shutil.rmtree ( TmpDir )
126if not os.path.isdir (TmpDir) : os.mkdir (TmpDir)
127TmpDirOCE = os.path.join (TmpDir, 'OCE')
128TmpDirICE = os.path.join (TmpDir, 'ICE')
129if not os.path.exists (TmpDirOCE) : os.mkdir (TmpDirOCE )
130if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE )
131
132echo ( f'Working in TMPDIR : {TmpDir}' )
133
134echo ( f'\nDealing with {L_EXP}' )
135
136#-- Model output directories
137if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' )
138if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' )
139dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
140dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
141
142echo ( f'The analysis relies on files from the following model output directories : ' )
143echo ( f'{dir_ATM_his}' )
144echo ( f'{dir_SRF_his}' )
145
146#-- Files Names
147if Freq == 'MO' : File = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M'
148if Freq == 'SE' : File = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M'
149
150echo ('\nOpen history files' )
151file_ATM_his = os.path.join (  dir_ATM_his, f'{File}_histmth.nc' )
152file_SRF_his = os.path.join (  dir_SRF_his, f'{File}_sechiba_history.nc' )
153
154d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
155d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
156
157echo ( file_ATM_his )
158echo ( file_SRF_his )
159
160## Compute run length
161dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
162echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
163dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
164
165## Compute length of each period
166dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
167echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
168dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
169dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
170
171#-- Open restart files
172YearRes = YearBegin - 1        # Year of the restart of beginning of simulation
173YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
174
175echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
176
177file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' )
178file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' )
179
180echo ( f'{file_restart_beg}' )
181echo ( f'{file_restart_end}' )
182
183file_ATM_beg = f'ATM_{JobName}_{YearRes}1231_restartphy.nc'
184file_ATM_end = f'ATM_{JobName}_{YearEnd}1231_restartphy.nc'
185if LMDZ :
186    file_DYN_beg = f'ATM_{JobName}_{YearRes}1231_restart.nc'
187    file_DYN_end = f'ATM_{JobName}_{YearEnd}1231_restart.nc'
188if ICO :
189    file_DYN_beg = f'ICO_{JobName}_{YearRes}1231_restart.nc'
190    file_DYN_end = f'ICO_{JobName}_{YearEnd}1231_restart.nc'
191file_SRF_beg = f'SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'
192file_SRF_end = f'SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc'
193liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg, ]
194liste_end = [file_ATM_end, file_DYN_end, file_SRF_end, ]
195echo ( f'{file_ATM_beg}')
196echo ( f'{file_ATM_end}')
197echo ( f'{file_SRF_beg}')
198echo ( f'{file_SRF_end}')
199
200if Routing == 'SIMPLE' : 
201    file_RUN_beg = f'SRF_{JobName}_{YearRes}1231_routing_restart.nc'
202    file_RUN_end = f'SRF_{JobName}_{YearEnd}1231_routing_restart.nc'
203    liste_beg.append ( file_RUN_beg )
204    liste_end.append ( file_RUN_end )
205    echo ( f'{file_RUN_beg}')
206    echo ( f'{file_RUN_end}')
207
208echo ('\nExtract restart files from tar : ATM, ICO and SRF')
209for resFile in liste_beg :
210    if not os.path.exists ( os.path.join (TmpDir, resFile) ) :
211        command =  f'cd {TmpDir} ; tar xf {file_restart_beg} {resFile}'
212        echo ( command )
213        os.system ( command )
214       
215for resFile in liste_end :
216    if not os.path.exists ( os.path.join (TmpDir, resFile) ) :
217        command =  f'cd {TmpDir} ; tar xf {file_restart_end} {resFile}'
218        echo ( command )
219        os.system ( command )
220
221echo ('\nOpening ATM SRF and ICO restart files')
222d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze()
223d_ATM_end = xr.open_dataset ( os.path.join (TmpDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze()
224d_SRF_beg = xr.open_dataset ( os.path.join (TmpDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze()
225d_SRF_end = xr.open_dataset ( os.path.join (TmpDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze()
226d_DYN_beg = xr.open_dataset ( os.path.join (TmpDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze()
227d_DYN_end = xr.open_dataset ( os.path.join (TmpDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze()
228
229for var in d_SRF_beg.variables :
230    #d_SRF_beg[var].attrs['_FillValue'] = 1.e20
231    #d_SRF_end[var].attrs['_FillValue'] = 1.e20
232    d_SRF_beg[var] = d_SRF_beg[var].where (  d_SRF_beg[var]<1.e20, 0.)
233    d_SRF_end[var] = d_SRF_end[var].where (  d_SRF_end[var]<1.e20, 0.)
234
235if ICO :
236    d_RUN_beg = xr.open_dataset ( os.path.join (TmpDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze()
237    d_RUN_end = xr.open_dataset ( os.path.join (TmpDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze()
238   
239echo ( file_ATM_beg )
240echo ( file_ATM_end )
241echo ( file_DYN_beg )
242echo ( file_DYN_end )
243echo ( file_SRF_beg )
244echo ( file_SRF_end )
245if Routing == 'SIMPLE' :
246    echo ( file_RUN_beg )
247    echo ( file_RUN_end )
248
249# ATM grid with cell surfaces
250if ICO : 
251    ATM_aire = d_ATM_his ['aire'][0]
252    ATM_fsea = d_ATM_his ['fract_ter'][Ø] + d_ATM_his ['fract_sic'][Ø]
253if LMDZ :
254    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0] )
255    ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_sic'][0] )
256    ATM_aire[0]  = np.sum ( d_ATM_his ['aire'][0, 0] )
257    ATM_aire[-1] = np.sum ( d_ATM_his ['aire'][0,-1] )
258
259if ICO :
260    file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
261    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze()
262    d_DYN_aire = d_DYN_aire.rename( {'cell':'cell_mesh'})
263    DYN_aire   = d_DYN_aire['aire']
264if LMDZ :
265    DYN_aire = ATM_aire
266
267#if LMDZ :
268#    d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} )
269
270ATM_aire_tot = np.sum (ATM_aire).values.item()
271ATM_aire_sea_tot = np.sum (ATM_aire*ATM_fsea).values.item()
272
273echo ( '\n------------------------------------------------------------------------------------' )
274echo ( '-- LMDZ changes in stores (for the records)' )
275
276#-- Change in precipitable water from the atmosphere daily and monthly files
277#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
278
279# ATM vertical grid
280ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
281ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
282klevp1 = ATM_Ahyb.shape[0]
283
284# Surface pressure
285if ICO : 
286    ATM_ps_beg = d_DYN_beg['ps']
287    ATM_ps_end = d_DYN_end['ps']
288   
289if LMDZ : 
290    ATM_ps_beg =  lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) )
291    ATM_ps_end =  lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) )
292   
293# 3D Pressure
294ATM_p_beg = ATM_Ahyb + ATM_Bhyb * ATM_ps_beg
295ATM_p_end = ATM_Ahyb + ATM_Bhyb * ATM_ps_end
296
297# Layer thickness
298ATM_sigma_beg = ATM_p_beg[0:-1]*0.
299ATM_sigma_end = ATM_p_end[0:-1]*0. 
300
301for k in np.arange (klevp1-1) :
302    ATM_sigma_beg[k,:] = (ATM_p_beg[k,:] - ATM_p_beg[k+1,:]) / g
303    ATM_sigma_end[k,:] = (ATM_p_end[k,:] - ATM_p_end[k+1,:]) / g
304
305ATM_sigma_beg = ATM_sigma_beg.rename ( {'klevp1':'sigs'} )
306ATM_sigma_end = ATM_sigma_end.rename ( {'klevp1':'sigs'} )
307
308##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases
309if LMDZ :
310    try : 
311        ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ) )
312        ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ) )
313    except :
314        ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) )
315        ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) )
316if ICO :
317    ATM_wat_beg = ( d_DYN_beg['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} )
318    ATM_wat_end = ( d_DYN_end['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} )
319   
320ATM_mas_wat_beg = np.sum (ATM_sigma_beg * ATM_wat_beg * ATM_aire).values.item()
321ATM_mas_wat_end = np.sum (ATM_sigma_end * ATM_wat_end * ATM_aire).values.item()
322
323dATM_mas_wat = ATM_mas_wat_end - ATM_mas_wat_beg
324
325echo ( 'Variation du contenu en eau atmosphere ' )
326echo ( 'ATM_mass_beg = {:12.6e} kg - ATM_mass_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) )
327echo ( 'dMass(atm)   = {:12.3e} kg '.format (dATM_mas_wat) )
328echo ( 'dMass(atm)   = {:12.3e} Sv '.format (dATM_mas_wat/dtime_sec*1.E-9) )
329echo ( 'dMass(atm)   = {:12.3e}m   '.format (dATM_mas_wat/ATM_aire_sea_tot*1E-3) )
330
331ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER']+d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']+d_ATM_beg['SNOW03']*d_ATM_beg['FOCE']+d_ATM_beg['SNOW04']*d_ATM_beg['FSIC']
332ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER']+d_ATM_end['SNOW02']*d_ATM_end['FLIC']+d_ATM_end['SNOW03']*d_ATM_end['FOCE']+d_ATM_end['SNOW04']*d_ATM_end['FSIC']
333if ICO :
334   ATM_sno_beg = ATM_sno_beg.rename ( {'points_physiques':'cell_mesh'} )
335   ATM_sno_end = ATM_sno_end.rename ( {'points_physiques':'cell_mesh'} )
336
337ATM_mas_sno_beg = np.sum ( ATM_sno_beg * DYN_aire ).values.item()
338ATM_mas_sno_end = np.sum ( ATM_sno_end * DYN_aire ).values.item()
339
340dATM_mas_sno = ATM_mas_sno_end - ATM_mas_sno_beg
341
342echo ( 'Variation du contenu en neige atmosphere ' )
343echo ( 'ATM_mas_sno_beg  = {:12.6e} kg  - ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
344echo ( 'dMass(neige atm) = {:12.3e} kg '.format (dATM_mas_sno) )
345echo ( 'dMass(neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1E-6/ICE_rho_ice) )
346echo ( 'dMass(neige atm) = {:12.3e} m  '.format (dATM_mas_sno/ATM_aire_sea_tot*1E-3) )
347
348echo ( '\nVariation du contenu en eau+neige atmosphere ' )
349echo ( 'dMass(eau + neige atm) = {:12.3e} kg '.format (  dATM_mas_wat + dATM_mas_sno) )
350echo ( 'dMass(eau + neige atm) = {:12.3e} Sv '.format ( (dATM_mas_wat + dATM_mas_sno)/dtime_sec*1E-9) )
351echo ( 'dMass(eau + neige atm) = {:12.3e} m  '.format ( (dATM_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot*1E-3) )
352
353echo ( '\n------------------------------------------------------------------------------------' )
354echo ( '-- SRF changes in routing reservoirs' )
355
356if Routing == 'SIMPLE' :
357    RUN_mas_wat_beg = np.sum ( d_RUN_beg ['fast_reservoir'] +  d_RUN_beg ['slow_reservoir'] +  d_RUN_beg ['stream_reservoir']).values.item()
358    RUN_mas_wat_end = np.sum ( d_RUN_end ['fast_reservoir'] +  d_RUN_end ['slow_reservoir'] +  d_RUN_end ['stream_reservoir']).values.item()
359
360if Routing == 'ORCHIDEE' : 
361    RUN_mas_wat_beg = np.sum ( d_SRF_beg['fastres']  + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \
362                             + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] ) .values.item()
363    RUN_mas_wat_end = np.sum ( d_SRF_end['fastres']  + d_SRF_end['slowres'] + d_SRF_end['streamres'] \
364                             + d_SRF_end['floodres'] + d_SRF_end['lakeres'] + d_SRF_end['pondres'] ) .values.item()
365
366dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg
367   
368echo ( '\nWater content in routing ' )
369echo ( 'RUN_mas_wat_beg = {:12.6e} kg '.format (RUN_mas_wat_beg) )
370echo ( 'RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end) )
371echo ( 'dMass(routing)  = {:12.3e} kg '.format(dRUN_mas_wat) )
372echo ( 'dMass(routing)  = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) )
373echo ( 'dMass(routing)  = {:12.3e} m  '.format(dRUN_mas_wat/ATM_aire_sea_tot*1E-3) )
374   
375SRF_mas_wat_beg = d_SRF_beg['tot_watveg_beg']+d_SRF_beg['tot_watsoil_beg'] + d_SRF_beg['snow_beg']
376SRF_mas_wat_end = d_SRF_end['tot_watveg_beg']+d_SRF_end['tot_watsoil_beg'] + d_SRF_end['snow_beg']
377
378if LMDZ : 
379    SRF_mas_wat_beg = SRF_mas_wat_beg.rename ( {'y':'points_phyiques'} )
380    SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'points_phyiques'} )
381if ICO : 
382    SRF_mas_wat_beg = SRF_mas_wat_beg.rename ( {'y':'cell_mesh'} )
383    SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'cell_mesh'} ) 
384
385SRF_mas_wat_beg = np.sum ( SRF_mas_wat_beg * SRF_aire ).values.item()
386SRF_mas_wat_end = np.sum ( SRF_mas_wat_end * SRF_aire ).values.item()
387
388dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg
389
390echo ( '\nWater content in and surface ' )
391echo ( 'SRF_mas_wat_beg  = {:12.6e} kg - SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) )
392echo ( 'dMass(water srf) = {:12.3e} kg '.format (dSRF_mas_wat) )
393echo ( 'dMass(water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-9) )
394echo ( 'dMass(water srf) = {:12.3e} m  '.format (dSRF_mas_wat/ATM_aire_sea_tot*1E-3) )
395
396echo ( '\n------------------------------------------------------------------------------------' )
397echo ( '\nWater content in  ATM + SRF + RUN ' )
398echo ( 'mas_wat_beg              = {:12.6e} kg '.format (ATM_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg) )
399echo ( 'mas_wat_end              = {:12.6e} kg '.format (ATM_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) )
400echo ( 'dMass(water atm+srf+run) = {:12.6e} kg '.format ( dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat) )
401echo ( 'dMass(water atm+srf+run) = {:12.3e} Sv '.format ((dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat)/dtime_sec*1E-9) )
402echo ( 'dMass(water atm+srf+run) = {:12.3e} m  '.format ((dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat)/ATM_aire_sea_tot*1E-3) )
403
404
405echo ( "\n-----------------" )
406echo ( " Atm -> Oce fluxes" )
407
408ATM_wbilo_sea = np.sum ( lmdz.geo2point (d_ATM_his['wbilo_oce'] + d_ATM_his['wbilo_sic'])*dtime_per_sec*ATM_aire ).values.item()
409ATM_runoff    = np.sum ( lmdz.geo2point (d_ATM_his['runofflic'])*ATM_aire*dtime_per_sec ).values.item()
410ATM_calving   = np.sum ( lmdz.geo2point (d_ATM_his['fqcalving'])*ATM_aire*dtime_per_sec ).values.item()
411
412echo (' wbilo oce+sic         = {:12.5e} (kg) '.format ( ATM_wbilo_sea ) )
413echo (' runoff liq            = {:12.5e} (kg) '.format ( ATM_runoff    ) )
414echo (' calving               = {:12.5e} (kg) '.format ( ATM_calving   ) )
415echo (' total                 = {:12.5e} (kg) '.format ( ATM_wbilo_sea - ATM_runoff - ATM_calving ) )
Note: See TracBrowser for help on using the repository browser.