source: TOOLS/WATER_BUDGET/CPL_waterbudget.py @ 6803

Last change on this file since 6803 was 6688, checked in by omamce, 5 months ago

O.M. :

WATER_BUDGET : put function common to ATM and OCE in WAterUtils library

use of dictionary for all input parameters
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 27.9 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
15SVN = {
16    'Author'  : "$Author$",
17    'Date'    : "$Date$",
18    'Revision': "$Revision$",
19    'Id'      : "$Id$",
20    'HeadURL' : "$HeadUrl:  $"
21    }
22###
23###
24## Import system modules
25import sys
26import os
27import configparser
28
29## Import needed scientific modules
30import numpy as np
31import xarray as xr
32
33## Import local modules
34import WaterUtils as wu
35import nemo, lmdz
36
37## Read command line arguments
38## ---------------------------
39print ( "Name of Python script:", sys.argv[0] )
40IniFile = sys.argv[1]
41
42# Test existence of IniFile
43if not os.path.exists (IniFile ) :
44    raise FileExistsError ( f"File not found : {IniFile = }" )
45
46if 'full' in IniFile or 'ATM' in IniFile :
47    FullIniFile = IniFile
48else :
49    FullIniFile = 'ATM_' + IniFile
50
51print ("Output file : ", FullIniFile )
52
53## Experiment parameters
54## --------------------
55dpar = wu.ReadConfig ( IniFile )
56
57## Configure all needed parameter from existant parameters
58## -------------------------------------------------------
59dpar = wu.SetDatesAndFiles ( dpar )
60
61## Output file with water budget diagnostics
62## -----------------------------------------
63f_out = dpar['Files']['f_out']
64
65## Put dpar values in local namespace
66## ----------------------------------
67for Section in dpar.keys () : 
68    print ( f'\nReading [{Section}]' )
69    for VarName in dpar[Section].keys() :
70        locals()[VarName] = dpar[Section][VarName]
71        print ( f'    {VarName:21} set to : {locals()[VarName]}' )
72       
73## Debuging and timer
74Timer = wu.functools.partial (wu.Timer, debug=Debug, timer=Timing)
75   
76## Useful functions
77## ----------------
78if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float :
79    def rprec (ptab) :
80        '''This version does nothing
81
82        rprec may be use to reduce floating precision when reading history files
83        '''
84        return ptab
85else                 :
86    def rprec (ptab) :
87        '''Returns float with a different precision'''
88        return ptab.astype(readPrec).astype(float)
89   
90def kg2Sv    (val, rho=ATM_RHO) :
91    '''From kg to Sverdrup'''
92    return val/dtime_sec*1.0e-6/rho
93
94def kg2myear (val, rho=ATM_RHO) :
95    '''From kg to m/year'''
96    return val/ATM.aire_sea_tot/rho/NbYear
97
98def var2prt (var, small=False, rho=ATM_RHO) :
99    '''Formats value for printing'''
100    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000
101    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
102
103def prtFlux (Desc, var, Form='F', small=False, rho=ATM_RHO, width=15) :
104    '''Pretty print of formattd value'''
105    if small :
106        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
107        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
108    else :
109        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
110        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
111    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) )
112    return None
113
114def echo (string, end='\n') :
115    '''Function to print to stdout *and* output file'''
116    print ( str(string), end=end  )
117    sys.stdout.flush ()
118    f_out.write ( str(string) + end )
119    f_out.flush ()
120    return None
121       
122d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
123if SECHIBA : 
124    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
125    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
126    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
127
128d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
129d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
130#d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
131d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
132if NEMO == '3.6' : d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
133
134echo ( f'{file_OCE_his = }' )
135echo ( f'{file_ICE_his = }' )
136echo ( f'{file_OCE_sca = }' )
137echo ( f'{file_OCE_srf = }' )
138
139## Compute run length
140## ------------------
141dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
142echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
143dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
144
145## Compute length of each period
146dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
147echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
148dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
149dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
150dtime_per_sec.attrs['unit'] = 's'
151
152# Number of years
153NbYear = dtime_sec / YEAR_LENGTH
154
155## Write the full configuration
156## ----------------------------
157params_out = open (FullIniFile, 'w', encoding = 'utf-8')
158params = wu.dict2config ( dpar )
159params.write ( params_out )
160params_out.close ()
161
162## Compute aire, fractions, etc ...
163## --------------------------------
164if ICO :
165    if not file_DYN_aire :
166        file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ResolAtm+'_grid.nc' )
167    dpar['Files']['file_DYN_aire'] = file_DYN_aire
168    echo ( f'{file_DYN_aire = }' )
169    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
170else :
171    d_DYN_aire = None
172       
173dpar, ATM = wu.ComputeGridATM ( dpar, d_ATM_his) 
174dpar, SRF = wu.ComputeGridSRF ( dpar, d_SRF_his)
175dpar, OCE = wu.ComputeGridOCE ( dpar, d_OCE_his, nperio=nperio )
176
177def ATM_flux_int (flux) :
178    '''Integrate (* time * surface) flux on atmosphere grid'''
179    zATM_flux_int  = wu.P1sum ( flux * dtime_per_sec * ATM.aire )
180    return zATM_flux_int
181
182def SRF_flux_int (flux) :
183    '''Integrate (* time * surface) flux on land grid'''
184    zSRF_flux_int  = wu.P1sum ( flux * dtime_per_sec * SRF.aire )
185    return zSRF_flux_int
186
187def ONE_stock_int (stock) :
188    '''Sum stock'''
189    zONE_stock_int =  wu.P1sum ( stock )
190    return zONE_stock_int
191
192def OCE_flux_int (flux) :
193    '''Integrate flux on oce grid'''
194    zOCE_flux_int = wu.P1sum ( flux * OCE.OCE_aire * dtime_per_sec )
195    return zOCE_flux_int
196
197def ONE_flux_int (flux) :
198    '''Integrate flux on oce grid'''
199    zOCE_flux_int = wu.P1sum ( flux * dtime_per_sec )
200    return zOCE_flux_int
201
202echo ( '\n====================================================================================' )
203echo ( f'-- ATM Fluxes  -- {Title} ' )
204
205if ATM_HIS == 'latlon' :
206    echo ( ' latlon case' )
207    ATM.wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' )
208    ATM.wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' )
209    ATM.wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' )
210    ATM.wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' )
211    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
212    ATM.fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' )
213    ATM.fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' )
214    ATM.precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' )
215    ATM.snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' )
216    ATM.evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' )
217    ATM.wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' )
218    ATM.wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' )
219    ATM.wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' )
220    ATM.wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' )
221    ATM.wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' )
222    ATM.wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' )
223    ATM.wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' )
224    ATM.wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' )
225    ATM.wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' )
226    ATM.wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' )
227    ATM.wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' )
228    ATM.wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' )
229    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
230    echo ( f'End of LATLON case')
231   
232if ATM_HIS == 'ico' :
233    echo (' ico case')
234    ATM.wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
235    ATM.wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
236    ATM.wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
237    ATM.wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
238    ATM.runofflic   = rprec (d_ATM_his ['runofflic'])
239    ATM.fqcalving   = rprec (d_ATM_his ['fqcalving'])
240    ATM.fqfonte     = rprec (d_ATM_his ['fqfonte']  )
241    ATM.precip      = rprec (d_ATM_his ['precip']   )
242    ATM.snowf       = rprec (d_ATM_his ['snow']     )
243    ATM.evap        = rprec (d_ATM_his ['evap']     )
244    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
245    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
246    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
247    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
248    ATM.runofflic   = rprec (d_ATM_his ['runofflic'])
249    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
250    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
251    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
252    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
253    ATM.wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
254    ATM.wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
255    ATM.wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
256    ATM.wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
257    ATM.wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
258    ATM.wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
259    ATM.wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
260    ATM.wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
261    echo ( f'End of ico case ')
262   
263echo ( 'ATM wprecip_oce' )
264ATM.wprecip_oce = ATM.wrain_oce + ATM.wsnow_oce
265ATM.wprecip_ter = ATM.wrain_ter + ATM.wsnow_ter
266ATM.wprecip_sic = ATM.wrain_sic + ATM.wsnow_sic
267ATM.wprecip_lic = ATM.wrain_lic + ATM.wsnow_lic
268
269ATM.wbilo       = ATM.wbilo_oce   + ATM.wbilo_sic   + ATM.wbilo_ter   + ATM.wbilo_lic
270ATM.wevap       = ATM.wevap_oce   + ATM.wevap_sic   + ATM.wevap_ter   + ATM.wevap_lic
271ATM.wprecip     = ATM.wprecip_oce + ATM.wprecip_sic + ATM.wprecip_ter + ATM.wprecip_lic
272ATM.wsnow       = ATM.wsnow_oce   + ATM.wsnow_sic   + ATM.wsnow_ter   + ATM.wsnow_lic
273ATM.wrain       = ATM.wrain_oce   + ATM.wrain_sic   + ATM.wrain_ter   + ATM.wrain_lic
274ATM.wemp        = ATM.wevap - ATM.wprecip
275ATM.emp         = ATM.evap  - ATM.precip
276
277ATM.wprecip_sea = ATM.wprecip_oce + ATM.wprecip_sic
278ATM.wsnow_sea   = ATM.wsnow_oce   + ATM.wsnow_sic
279ATM.wrain_sea   = ATM.wrain_oce   + ATM.wrain_sic
280ATM.wbilo_sea   = ATM.wbilo_oce   + ATM.wbilo_sic
281ATM.wevap_sea   = ATM.wevap_sic   + ATM.wevap_oce
282
283ATM.wemp_ter    = ATM.wevap_ter - ATM.wprecip_ter
284ATM.wemp_oce    = ATM.wevap_oce - ATM.wprecip_oce
285ATM.wemp_sic    = ATM.wevap_sic - ATM.wprecip_sic
286ATM.wemp_lic    = ATM.wevap_lic - ATM.wprecip_lic
287ATM.wemp_sea    = ATM.wevap_sic - ATM.wprecip_oce
288
289if SECHIBA : 
290    if RUN_HIS == 'latlon' :
291        echo ( f'RUN costalflow Grille LATLON' )
292        if TestInterp :
293            echo ( f'RUN runoff TestInterp' )
294            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' )
295            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' )
296        else : 
297            echo ( f'RUN runoff' )
298            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' )
299            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' )
300           
301        SRF.RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' ) 
302        SRF.RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' )
303        SRF.RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' )
304        SRF.RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' ) 
305        SRF.RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' )
306       
307    if RUN_HIS == 'ico' :
308        echo ( f'RUN costalflow Grille ICO' )
309       
310        SRF.RUN_coastalflow     =  rprec (d_RUN_his ['coastalflow'])
311        SRF.RUN_riverflow       =  rprec (d_RUN_his ['riverflow']  )
312        SRF.RUN_runoff          =  rprec (d_RUN_his ['runoff']     )
313        SRF.RUN_drainage        =  rprec (d_RUN_his ['drainage']   )
314        SRF.RUN_riversret       =  rprec (d_RUN_his ['riversret']  )
315        SRF.RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
316        SRF.RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
317       
318    ## Correcting units of SECHIBA variables
319    def mmd2SI ( Var ) :
320        '''Change unit from mm/d or m^3/s to kg/s if needed'''
321        if 'units' in Var.attrs : 
322            if Var.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
323                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg/s'
324            if Var.attrs['units'] == 'mm/d' :
325                Var.values = Var.values  * ATM_RHO * (1e-3/86400.) ;  Var.attrs['units'] = 'kg/s'
326            if Var.attrs['units'] in ['m^3', 'm3', ] :
327                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg'
328        return Var
329           
330    SRF.RUN_coastalflow     = mmd2SI ( SRF.RUN_coastalflow )
331    SRF.RUN_coastalflow_cpl = mmd2SI ( SRF.RUN_coastalflow_cpl )
332    SRF.RUN_drainage        = mmd2SI ( SRF.RUN_drainage )
333    SRF.RUN_riverflow       = mmd2SI ( SRF.RUN_riverflow )
334    SRF.RUN_riverflow_cpl   = mmd2SI ( SRF.RUN_riverflow_cpl )
335    SRF.RUN_riversret       = mmd2SI ( SRF.RUN_riversret )
336    SRF.RUN_runoff          = mmd2SI ( SRF.RUN_runoff )
337
338    SRF.RUN_input  = SRF.RUN_runoff      + SRF.RUN_drainage
339    SRF.RUN_output = SRF.RUN_coastalflow + SRF.RUN_riverflow
340   
341echo ( f'ATM flw_wbilo' )
342ATM.flx_wbilo       = ATM_flux_int ( ATM.wbilo      )
343ATM.flx_wevap       = ATM_flux_int ( ATM.wevap      )
344ATM.flx_wprecip     = ATM_flux_int ( ATM.wprecip    )
345ATM.flx_wsnow       = ATM_flux_int ( ATM.wsnow      )
346ATM.flx_wrain       = ATM_flux_int ( ATM.wrain      )
347ATM.flx_wemp        = ATM_flux_int ( ATM.wemp       )
348
349ATM.flx_wbilo_lic   = ATM_flux_int ( ATM.wbilo_lic  )
350ATM.flx_wbilo_oce   = ATM_flux_int ( ATM.wbilo_oce  )
351ATM.flx_wbilo_sea   = ATM_flux_int ( ATM.wbilo_sea  )
352ATM.flx_wbilo_sic   = ATM_flux_int ( ATM.wbilo_sic  )
353ATM.flx_wbilo_ter   = ATM_flux_int ( ATM.wbilo_ter  )
354ATM.flx_calving     = ATM_flux_int ( ATM.fqcalving  )
355ATM.flx_fqfonte     = ATM_flux_int ( ATM.fqfonte    )
356
357ATM.LIC_flx_calving     = ATM_flux_int ( ATM.fqcalving  )
358ATM.LIC_flx_fqfonte     = ATM_flux_int ( ATM.fqfonte    )
359
360ATM.flx_precip      = ATM_flux_int ( ATM.precip     )
361ATM.flx_snowf       = ATM_flux_int ( ATM.snowf      )
362ATM.flx_evap        = ATM_flux_int ( ATM.evap       )
363ATM.flx_runlic      = ATM_flux_int ( ATM.runofflic  )
364
365ATM.flx_wrain_ter   = ATM_flux_int ( ATM.wrain_ter )
366ATM.flx_wrain_oce   = ATM_flux_int ( ATM.wrain_oce )
367ATM.flx_wrain_lic   = ATM_flux_int ( ATM.wrain_lic )
368ATM.flx_wrain_sic   = ATM_flux_int ( ATM.wrain_sic )
369ATM.flx_wrain_sea   = ATM_flux_int ( ATM.wrain_sea )
370
371ATM.flx_wsnow_ter   = ATM_flux_int ( ATM.wsnow_ter )
372ATM.flx_wsnow_oce   = ATM_flux_int ( ATM.wsnow_oce )
373ATM.flx_wsnow_lic   = ATM_flux_int ( ATM.wsnow_lic )
374ATM.flx_wsnow_sic   = ATM_flux_int ( ATM.wsnow_sic )
375ATM.flx_wsnow_sea   = ATM_flux_int ( ATM.wsnow_sea )
376
377ATM.flx_wevap_ter   = ATM_flux_int ( ATM.wevap_ter )
378ATM.flx_wevap_oce   = ATM_flux_int ( ATM.wevap_oce )
379ATM.flx_wevap_lic   = ATM_flux_int ( ATM.wevap_lic )
380ATM.flx_wevap_sic   = ATM_flux_int ( ATM.wevap_sic )
381ATM.flx_wevap_sea   = ATM_flux_int ( ATM.wevap_sea )
382ATM.flx_wprecip_lic = ATM_flux_int ( ATM.wprecip_lic )
383ATM.flx_wprecip_oce = ATM_flux_int ( ATM.wprecip_oce )
384ATM.flx_wprecip_sic = ATM_flux_int ( ATM.wprecip_sic )
385ATM.flx_wprecip_ter = ATM_flux_int ( ATM.wprecip_ter )
386ATM.flx_wprecip_sea = ATM_flux_int ( ATM.wprecip_sea )
387ATM.flx_wemp_lic    = ATM_flux_int ( ATM.wemp_lic )
388ATM.flx_wemp_oce    = ATM_flux_int ( ATM.wemp_oce )
389ATM.flx_wemp_sic    = ATM_flux_int ( ATM.wemp_sic )
390ATM.flx_wemp_ter    = ATM_flux_int ( ATM.wemp_ter )
391ATM.flx_wemp_sea    = ATM_flux_int ( ATM.wemp_sea )
392
393ATM.flx_emp          = ATM_flux_int ( ATM.emp )
394
395if SECHIBA : 
396    SRF.RUN_flx_coastal     = ONE_flux_int ( SRF.RUN_coastalflow)
397    SRF.RUN_flx_river       = ONE_flux_int ( SRF.RUN_riverflow  )
398    SRF.RUN_flx_coastal_cpl = ONE_flux_int ( SRF.RUN_coastalflow_cpl)
399    SRF.RUN_flx_river_cpl   = ONE_flux_int ( SRF.RUN_riverflow_cpl  )
400    SRF.RUN_flx_drainage    = SRF_flux_int ( SRF.RUN_drainage   )
401    SRF.RUN_flx_riversret   = SRF_flux_int ( SRF.RUN_riversret  )
402    SRF.RUN_flx_runoff      = SRF_flux_int ( SRF.RUN_runoff     )
403    SRF.RUN_flx_input       = SRF_flux_int ( SRF.RUN_input      )
404    SRF.RUN_flx_output      = ONE_flux_int ( SRF.RUN_output     )
405   
406    SRF.RUN_flx_bil    = ONE_flux_int ( SRF.RUN_input       - SRF.RUN_output)
407    SRF.RUN_flx_rivcoa = ONE_flux_int ( SRF.RUN_coastalflow + SRF.RUN_riverflow)
408   
409prtFlux ('wbilo_oce            ', ATM.flx_wbilo_oce     , 'f' )         
410prtFlux ('wbilo_sic            ', ATM.flx_wbilo_sic     , 'f' )         
411prtFlux ('wbilo_sic+oce        ', ATM.flx_wbilo_sea     , 'f' )         
412prtFlux ('wbilo_ter            ', ATM.flx_wbilo_ter     , 'f' )         
413prtFlux ('wbilo_lic            ', ATM.flx_wbilo_lic     , 'f' )         
414prtFlux ('Sum wbilo_*          ', ATM.flx_wbilo         , 'f', True) 
415prtFlux ('E-P                  ', ATM.flx_emp           , 'f', True) 
416prtFlux ('calving              ', ATM.flx_calving       , 'f' ) 
417prtFlux ('fqfonte              ', ATM.flx_fqfonte       , 'f' )       
418prtFlux ('precip               ', ATM.flx_precip        , 'f' )       
419prtFlux ('snowf                ', ATM.flx_snowf         , 'f' )       
420prtFlux ('evap                 ', ATM.flx_evap          , 'f' )
421prtFlux ('runoff lic           ', ATM.flx_runlic        , 'f' )
422
423prtFlux ('ATM.flx_wevap*       ', ATM.flx_wevap         , 'f' )
424prtFlux ('ATM.flx_wrain*       ', ATM.flx_wrain         , 'f' )
425prtFlux ('ATM.flx_wsnow*       ', ATM.flx_wsnow         , 'f' )
426prtFlux ('ATM.flx_wprecip*     ', ATM.flx_wprecip       , 'f' )
427prtFlux ('ATM.flx_wemp*        ', ATM.flx_wemp          , 'f', True )
428
429prtFlux ('ERROR evap           ', ATM.flx_wevap   - ATM.flx_evap  , 'e', True )
430prtFlux ('ERROR precip         ', ATM.flx_wprecip - ATM.flx_precip, 'e', True )
431prtFlux ('ERROR snow           ', ATM.flx_wsnow   - ATM.flx_snowf , 'e', True )
432prtFlux ('ERROR emp            ', ATM.flx_wemp    - ATM.flx_emp   , 'e', True )
433
434if SECHIBA : 
435    echo ( '\n====================================================================================' )
436    echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
437    prtFlux ('coastalflow   ', SRF.RUN_flx_coastal    , 'f' ) 
438    prtFlux ('riverflow     ', SRF.RUN_flx_river      , 'f' )       
439    prtFlux ('coastal_cpl   ', SRF.RUN_flx_coastal_cpl, 'f' ) 
440    prtFlux ('riverf_cpl    ', SRF.RUN_flx_river_cpl  , 'f' )   
441    prtFlux ('river+coastal ', SRF.RUN_flx_rivcoa     , 'f' )   
442    prtFlux ('drainage      ', SRF.RUN_flx_drainage   , 'f' )   
443    prtFlux ('riversret     ', SRF.RUN_flx_riversret  , 'f' )   
444    prtFlux ('runoff        ', SRF.RUN_flx_runoff     , 'f' )   
445    prtFlux ('river in      ', SRF.RUN_flx_input      , 'f' )   
446    prtFlux ('river out     ', SRF.RUN_flx_output     , 'f' )   
447    prtFlux ('river bil     ', SRF.RUN_flx_bil        , 'f' ) 
448   
449echo ( '\n====================================================================================' )
450echo ( f'-- OCE Fluxes  -- {Title} ' )
451
452# Read variable and computes integral over space and time
453OCE.OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE.OCE_mass_empmr    = OCE_flux_int ( OCE.OCE_empmr    )
454OCE.OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE.OCE_mass_wfob     = OCE_flux_int ( OCE.OCE_wfob     )
455OCE.OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE.OCE_mass_emp_oce  = OCE_flux_int ( OCE.OCE_emp_oce  )
456OCE.OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE.OCE_mass_emp_ice  = OCE_flux_int ( OCE.OCE_emp_ice  )
457OCE.OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE.OCE_mass_iceshelf = OCE_flux_int ( OCE.OCE_iceshelf )
458OCE.OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE.OCE_mass_calving  = OCE_flux_int ( OCE.OCE_calving  )
459OCE.OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE.OCE_mass_iceberg  = OCE_flux_int ( OCE.OCE_iceberg  )
460OCE.OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE.OCE_mass_friver   = OCE_flux_int ( OCE.OCE_friver   )
461OCE.OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE.OCE_mass_runoffs  = OCE_flux_int ( OCE.OCE_runoffs  )
462if NEMO == 4.0 or NEMO == 4.2 :
463    OCE.OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE.OCE_mass_wfxice   = OCE_flux_int ( OCE.OCE_wfxice )
464    OCE.OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE.OCE_mass_wfxsnw   = OCE_flux_int ( OCE.OCE_wfxsnw )
465    OCE.OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE.OCE_mass_wfxsub   = OCE_flux_int ( OCE.OCE_wfxsub )
466if NEMO == 3.6 :
467    OCE.OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*OCE.ICE_RHO_ICE ; OCE.OCE_mass_wfxice   = OCE_flux_int ( OCE.OCE_wfxice )
468    OCE.OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*OCE.ICE_RHO_SNO ; OCE.OCE_mass_wfxsnw   = OCE_flux_int ( OCE.OCE_wfxsnw )
469    OCE.OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_RHO_SNO     ; OCE.OCE_mass_wfxsub   = OCE_flux_int ( OCE.OCE_wfxsub )
470# Additional checks
471OCE.OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE.OCE_mass_evap_oce   = OCE_flux_int ( OCE.OCE_evap_oce )
472OCE.ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; OCE.ICE_mass_evap_ice   = OCE_flux_int ( OCE.ICE_evap_ice )
473OCE.OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE.OCE_mass_snow_oce   = OCE_flux_int ( OCE.OCE_snow_oce )
474OCE.OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE.OCE_mass_snow_ice   = OCE_flux_int ( OCE.OCE_snow_ice )
475OCE.OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE.OCE_mass_rain       = OCE_flux_int ( OCE.OCE_rain     )
476OCE.ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; OCE.ICE_mass_wfxsub_err = OCE_flux_int ( OCE.ICE_wfxsub_err )
477if NEMO == 4.0 or NEMO == 4.2 :
478    OCE.ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; OCE.ICE_mass_wfxpnd     = OCE_flux_int ( OCE.ICE_wfxpnd     )
479    OCE.ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; OCE.ICE_mass_wfxsnw_sub = OCE_flux_int ( OCE.ICE_wfxsnw_sub )
480    OCE.ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; OCE.ICE_mass_wfxsnw_pre = OCE_flux_int ( OCE.ICE_wfxsnw_pre )
481if NEMO == 3.6 :
482    OCE.ICE_wfxpnd = 0.0 ; OCE.ICE_mass_wfxpnd = 0.0
483    OCE.ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_RHO_SNO  ; OCE.ICE_mass_wfxsnw_sub = OCE_flux_int ( OCE.ICE_wfxsnw_sub )
484    OCE.ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_RHO_SNO  ; OCE.ICE_mass_wfxsnw_pre = OCE_flux_int ( OCE.ICE_wfxsnw_pre     )
485
486OCE.OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE.OCE_mass_wfcorr  = OCE_flux_int ( OCE.OCE_wfcorr )
487if OCE_relax  :
488    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
489    OCE.OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE.OCE_mass_vflx_fwb   = OCE_flux_int ( OCE.OCE_vflx_fwb )
490    OCE.OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE.OCE_mass_vflx_ssr   = OCE_flux_int ( OCE.OCE_vflx_ssr )
491else : 
492    OCE.OCE_fwb = 0.0 ; OCE.OCE_mass_fwb = 0.0
493    OCE.OCE_ssr = 0.0 ; OCE.OCE_mass_ssr = 0.0
494if OCE_icb :
495    OCE.OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE.OCE_mass_berg_icb = OCE_flux_int ( OCE.OCE_berg_icb    )
496    OCE.OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE.OCE_mass_calv_icb = OCE_flux_int ( OCE.OCE_calving_icb )
497else :
498    OCE.OCE_berg_icb = 0. ; OCE.OCE_mass_berg_icb = 0.
499    OCE.OCE_calv_icb = 0. ; OCE.OCE_mass_calv_icb = 0.
500
501OCE.OCE_mass_emp = OCE.OCE_mass_emp_oce - OCE.OCE_mass_wfxice  - OCE.OCE_mass_wfxsnw  - OCE.ICE_mass_wfxpnd - OCE.ICE_mass_wfxsub_err
502OCE.OCE_mass_all = OCE.OCE_mass_emp_oce + OCE.OCE_mass_emp_ice - OCE.OCE_mass_runoffs - OCE.OCE_mass_iceshelf
503
504prtFlux ('OCE+ICE budget ', OCE.OCE_mass_all       , 'e', True)
505prtFlux ('  EMPMR        ', OCE.OCE_mass_empmr     , 'e', True)
506prtFlux ('  WFOB         ', OCE.OCE_mass_wfob      , 'e', True)
507prtFlux ('  EMP_OCE      ', OCE.OCE_mass_emp_oce   , 'e', True)
508prtFlux ('  EMP_ICE      ', OCE.OCE_mass_emp_ice   , 'e', True)
509prtFlux ('  EMP          ', OCE.OCE_mass_emp       , 'e', True)
510prtFlux ('  ICEBERG      ', OCE.OCE_mass_iceberg   , 'e',     )
511prtFlux ('  ICESHELF     ', OCE.OCE_mass_iceshelf  , 'e', True)
512prtFlux ('  CALVING      ', OCE.OCE_mass_calving   , 'e', True)
513prtFlux ('  FRIVER       ', OCE.OCE_mass_friver    , 'e',     ) 
514prtFlux ('  RUNOFFS      ', OCE.OCE_mass_runoffs   , 'e', True)
515prtFlux ('  WFXICE       ', OCE.OCE_mass_wfxice    , 'e', True)
516prtFlux ('  WFXSNW       ', OCE.OCE_mass_wfxsnw    , 'e', True)
517prtFlux ('  WFXSUB       ', OCE.OCE_mass_wfxsub    , 'e', True)
518prtFlux ('  WFXPND       ', OCE.ICE_mass_wfxpnd    , 'e', True)
519prtFlux ('  WFXSNW_SUB   ', OCE.ICE_mass_wfxsnw_sub, 'e', True)
520prtFlux ('  WFXSNW_PRE   ', OCE.ICE_mass_wfxsnw_pre, 'e', True)
521prtFlux ('  WFXSUB_ERR   ', OCE.ICE_mass_wfxsub_err, 'e', True)
522prtFlux ('  EVAP_OCE     ', OCE.OCE_mass_evap_oce  , 'e'      )
523prtFlux ('  EVAP_ICE     ', OCE.ICE_mass_evap_ice  , 'e', True)
524prtFlux ('  SNOW_OCE     ', OCE.OCE_mass_snow_oce  , 'e', True)
525prtFlux ('  SNOW_ICE     ', OCE.OCE_mass_snow_ice  , 'e', True)
526prtFlux ('  RAIN         ', OCE.OCE_mass_rain      , 'e'      )
527prtFlux ('  FWB          ', OCE.OCE_mass_fwb       , 'e', True)
528prtFlux ('  SSR          ', OCE.OCE_mass_ssr       , 'e', True)
529prtFlux ('  WFCORR       ', OCE.OCE_mass_wfcorr    , 'e', True)
530prtFlux ('  BERG_ICB     ', OCE.OCE_mass_berg_icb  , 'e', True)
531prtFlux ('  CALV_ICB     ', OCE.OCE_mass_calv_icb  , 'e', True)
532
533
534echo (' ')
535
536prtFlux ( 'wbilo sea  ', ATM_flux_int (ATM.wbilo_sea), 'e',   )
537if SECHIBA : 
538    prtFlux ( 'costalflow ', ONE_flux_int (SRF.RUN_coastalflow), 'e',   )
539    prtFlux ( 'riverflow  ', SRF.RUN_flx_river  , 'e',   )
540    prtFlux ( 'costalflow ', SRF.RUN_flx_coastal, 'e',   )
541    prtFlux ( 'runoff     ', SRF.RUN_flx_river+SRF.RUN_flx_coastal, 'e',   )
542
543#ATM.to_OCE   =  ATM_flux_int (ATM.wbilo_sea) - SRF.RUN_flx_river - SRF.RUN_flx_coastal - ATM.flx_calving
544ATM.to_OCE   =  ATM_flux_int (ATM.wbilo_sea) - SRF.RUN_flx_river - SRF.RUN_flx_coastal - ATM.flx_calving
545#OCE.OCE_from_ATM = -OCE.OCE_mass_emp_oce - OCE.OCE_mass_emp_ice + OCE.OCE_mass_runoffs + OCE.OCE_mass_iceberg + OCE.OCE_mass_calving + OCE.OCE_mass_iceshelf
546OCE.OCE_from_ATM = OCE.OCE_mass_all
547
548prtFlux ( 'ATM.to_OCE      ', ATM.to_OCE      , 'e', True )
549prtFlux ( 'OCE.OCE_from_ATM', OCE.OCE_from_ATM, 'e', True )
550
551echo ( ' ' )
552echo ( f'{Title = }' )
553
554echo ( 'SVN Information' )
555for kk in SVN.keys():
556    print ( SVN[kk] )
557
558## Write the full configuration
559## ----------------------------
560params_out = open (FullIniFile, 'w', encoding = 'utf-8')
561params = wu.dict2config ( dpar )
562params.write ( params_out )
563params_out.close ()
Note: See TracBrowser for help on using the repository browser.