Index: TOOLS/WATER_BUDGET/OCE_waterbudget.py
===================================================================
--- TOOLS/WATER_BUDGET/OCE_waterbudget.py	(revision 6620)
+++ TOOLS/WATER_BUDGET/OCE_waterbudget.py	(revision 6647)
@@ -39,4 +39,6 @@
 import nemo, lmdz
 
+from WaterUtils import VarInt, Rho, Ra, Grav, ICE_rho_ice, ICE_rho_sno, OCE_rho_liq, ATM_rho, SRF_rho, RUN_rho, ICE_rho_pnd, YearLength
+
 ## Creates parser for reading .ini input file
 ## -------------------------------------------
@@ -47,6 +49,7 @@
 ## ---------------------
 ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
-OCE_icb=False ; Coupled=False ; Routing=None
+OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None
 TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
+YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None
 
 ##
@@ -66,4 +69,11 @@
 tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
 ContinueOnError=False ; ErrorCount=0 ; SortIco = False
+
+##
+## Precision of history file reading
+## ---------------------------------
+# Default is float (full precision). Degrade the precision by using np.float32
+# Restart file are always read at the full precision
+readPrec=float
 
 ## Read command line arguments
@@ -137,5 +147,14 @@
 
 print ( f'\nConfig file readed : {IniFile} ' )
-            
+
+##
+## Reading prec
+if wu.unDefined ( 'readPrec' )  :
+    readPrec = np.float64
+else :
+    if readPrec in ["float", "float64", "r8", "double", "<class 'float'>"         ] : readPrec = float
+    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32
+    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16
+              
 ## Some physical constants
 ## =======================
@@ -155,8 +174,8 @@
 if not 'Files' in config.keys () : config['Files'] = {}
 
-config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno,
-                      'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
-
-config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco }
+config['Physics'] = { 'Ra':str(Ra), 'Grav':str(Grav), 'ICE_rho_ice':str(ICE_rho_ice), 'ICE_rho_sno':str(ICE_rho_sno),
+                      'OCE_rho_liq':str(OCE_rho_liq), 'ATM_rho':str(ATM_rho), 'SRF_rho':str(SRF_rho), 'RUN_rho':str(RUN_rho)}
+
+config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco, 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) }
 
 ## --------------------------
@@ -178,14 +197,15 @@
     config['Experiment']['DateBegin'] = DateBegin
 else :
-    YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin )
+    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin )
     DateBegin = wu.FormatToGregorian ( DateBegin)
-    config['Experiment']['YearBegin'] = YearBegin
+    config['Experiment']['YearBegin'] = str(YearBegin)
 
 if wu.unDefined ( 'DateEnd' )   :
     DateEnd   = f'{YearEnd}1231'
-    config['Experiment']['DateEnd']   = DateEnd
+    config['Experiment']['DateEnd']   = str(DateEnd)
 else :
-    YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd )
+    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd )
     DateEnd   = wu.FormatToGregorian ( DateEnd)
+    config['Experiment']['DateEnd']   = str(DateEnd)
 
 if wu.unDefined ( 'PackFrequency' ) :
@@ -204,6 +224,10 @@
 ## Output file with water budget diagnostics
 ## -----------------------------------------
-if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out'
-config['Files']['FileOut'] = FileOut
+if wu.unDefined ( 'FileOut' ) :
+    FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}'
+    if readPrec == np.float32 : FileOut = f'{FileOut}_float32'
+    FileOut = f'{FileOut}.out'
+
+    config['Files']['FileOut'] = FileOut
 
 f_out = open ( FileOut, mode = 'w' )
@@ -211,17 +235,24 @@
 ## Useful functions
 ## ----------------
-def kg2Sv  (val, rho=OCE_rho_liq) :
+
+# Degrades precision 
+if readPrec == float :
+    def rprec (tab) : return tab
+else : 
+    def rprec (tab) : return tab.astype(readPrec).astype(float)
+        
+def kg2Sv  (val, rho=ATM_rho) :
     '''From kg to Sverdrup'''
     return val/dtime_sec*1.0e-6/rho
 
-def kg2myear (val, rho=OCE_rho_liq) :
+def kg2myear (val, rho=ATM_rho) :
     '''From kg to m/year'''
     return val/OCE_aire_tot/rho/NbYear
 
-def var2prt (var, small=False, rho=OCE_rho_liq) :
+def var2prt (var, small=False, rho=ATM_rho) :
     if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
     else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
 
-def prtFlux (Desc, var, Form='F', small=False, rho=OCE_rho_liq, width=15) :
+def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
     if small :
         if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
@@ -330,7 +361,7 @@
 if NEMO == 3.6 :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
 
-echo ( file_OCE_his )
-echo ( file_OCE_sca )
-echo ( file_ICE_his )
+echo ( f'{file_OCE_his = }' )
+echo ( f'{file_OCE_sca = }' )
+echo ( f'{file_ICE_his = }' )
 
 ## Compute run length
@@ -666,48 +697,48 @@
 
 # Read variable and computes integral over space and time
-OCE_empmr      = d_OCE_his['wfo']         ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
-OCE_wfob       = d_OCE_his['wfob']        ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
-OCE_emp_oce    = d_OCE_his['emp_oce']     ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
-OCE_emp_ice    = d_OCE_his['emp_ice']     ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
-OCE_iceshelf   = d_OCE_his['iceshelf']    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
-OCE_calving    = d_OCE_his['calving']     ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
-OCE_iceberg    = d_OCE_his['iceberg']     ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
-OCE_friver     = d_OCE_his['friver']      ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
-OCE_runoffs    = d_OCE_his['runoffs']     ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
+OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
+OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
+OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
+OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
+OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
+OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
+OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
+OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
+OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
 if NEMO == 4.0 or NEMO == 4.2 :
-    OCE_wfxice     = d_OCE_his['vfxice'] ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
-    OCE_wfxsnw     = d_OCE_his['vfxsnw'] ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
-    OCE_wfxsub     = d_OCE_his['vfxsub'] ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
+    OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
+    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
+    OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
 if NEMO == 3.6 :
-    OCE_wfxice     = d_OCE_his['vfxice']/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
-    OCE_wfxsnw     = d_OCE_his['vfxsnw']/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
-    OCE_wfxsub     = d_OCE_his['vfxsub']/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
+    OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
+    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
+    OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
 # Additional checks
-OCE_evap_oce   = d_OCE_his['evap_ao_cea'] ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
-ICE_evap_ice   = d_OCE_his['subl_ai_cea'] ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
-OCE_snow_oce   = d_OCE_his['snow_ao_cea'] ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
-OCE_snow_ice   = d_OCE_his['snow_ai_cea'] ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
-OCE_rain       = d_OCE_his['rain']        ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
-ICE_wfxsub_err = d_ICE_his['vfxsub_err']  ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
+OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
+ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
+OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
+OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
+OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
+ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
 if NEMO == 4.0 or NEMO == 4.2 :
-    ICE_wfxpnd     = d_ICE_his['vfxpnd']     ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
-    ICE_wfxsnw_sub = d_ICE_his['vfxsnw_sub'] ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
-    ICE_wfxsnw_pre = d_ICE_his['vfxsnw_pre'] ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
+    ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
+    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
+    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
 if NEMO == 3.6 :
     ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
-    ICE_wfxsnw_sub = d_ICE_his['vfxsub']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
-    ICE_wfxsnw_pre = d_ICE_his['vfxspr']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
-
-OCE_wfcorr    = d_OCE_his['wfcorr'] ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
+    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
+    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
+
+OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
 if OCE_relax  :
     # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
-    OCE_vflx_fwb = d_OCE_his['vflx_fwb'] ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
-    OCE_vflx_ssr = d_OCE_his['vflx_ssr'] ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
+    OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
+    OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
 else : 
     OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0
     OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0
 if OCE_icb :
-    OCE_berg_icb    = d_OCE_his['berg_floating_melt'] ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
-    OCE_calving_icb = d_OCE_his['calving_icb']        ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
+    OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
+    OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
 else :
     OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0.
@@ -715,6 +746,5 @@
 
 OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err
-
-OCE_mas_all = OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
+OCE_mas_all = OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
 
 echo ('\n-- Fields:' ) 
@@ -748,4 +778,11 @@
 prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
 
+echo (' ')
+if Coupled :
+    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceshelf, 'e', True )
+else :
+    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf, 'e', True )
+
+
 echo     ('\n===>  Comparing ===>' ) 
 echo     ('  WFX OCE                     = -empmr + iceshelf                                            = {:12.5e} (kg) '.format ( -OCE_mas_empmr + OCE_mas_iceshelf) )
@@ -770,7 +807,7 @@
 echo ('                        = (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) ) )
 
-prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e' )
-prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e' )
-prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e' )
+prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e', True )
+prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e', True )
+prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e', True )
 
 
