Changeset 6647 for TOOLS/WATER_BUDGET/CPL_waterbudget.py

Timestamp:

10/10/23 12:58:04 (7 months ago)

Author:

omamce

Message:

O.M. :

New version of WATER_BUDGET

• Conservation in NEMO are OK
• Conservation in Sechiba are OK, for both ICO and latlon grids
• Problems in atmosphere, LIC surface and ocean/atmosphere coherence

File:

• TOOLS/WATER_BUDGET/CPL_waterbudget.py (modified) (19 diffs)

TOOLS/WATER_BUDGET/CPL_waterbudget.py

@@ -12 +12 @@
 ##
 ##
-## SVN information
+# SVN information
 #  $Author$
 #  $Date$
@@ -19 +19 @@
 #  $HeadURL$
 
-
+# SVN Information
+SVN = {
+    'Author'  : "$Author$",
+    'Date'    : "$Date$",
+    'Revision': "$Revision$",
+    'Id'      : "$Id$",
+    'HeadURL' : "$HeadUrl:  $"
+    }
 ###
 ## Import system modules
@@ -27 +34 @@
 from pathlib import Path
 
+###
+## Import system modules
+import sys, os, shutil#, subprocess, platform
+import configparser, re
+
+## Import needed scientific modules
+import numpy as np, xarray as xr
+
 # Check python version
 if sys.version_info < (3, 8, 0) :
@@ -32 +47 @@
     raise Exception ( "Minimum Python version is 3.8" )
 
-## Import local module
+## Import local modules
 import WaterUtils as wu
+import libIGCM_sys
+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
-config = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
+## -------------------------------------------
+config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
 config.optionxform = str # To keep capitals
 
 ## Experiment parameters
-ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ;  TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
+## ---------------------
+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 ; TestInterp=None
+TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
+YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None
 
 ##
-ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
-TmpDir=None ; RunDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None
-file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
-tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
-file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_beg=None ; file_OCE_end=None ; file_OCE_srf=None ; file_OCE_sca=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
-
-# Read command line arguments
+ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
+FileDir=None ; FileOut=None
+dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
+FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
+file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None ; file_OCE_srf=None
+tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
+file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
+file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
+file_ICE_beg=None ; file_OCE_beg=None
+file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
+tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
+tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
+tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
+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
+## ---------------------------
 print ( "Name of Python script:", sys.argv[0] )
 IniFile = sys.argv[1]
-# Text existence of IniFile
+
+# Test existence of IniFile
 if not os.path.exists (IniFile ) :
     raise FileExistsError ( f"File not found : {IniFile = }" )
@@ -63 +106 @@
 FullIniFile = 'full_' + IniFile
 
+## Reading config.card if possible
+## -------------------------------
+ConfigCard = None
+
+if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
+    if 'ConfigCard' in config['Experiment'].keys () :
+        ConfigCard = config['Experiment']['ConfigCard']
+        print ( f'{ConfigCard=}' )
+
+if ConfigCard : ## Read config card if it exists
+    # Text existence of ConfigCard
+    if os.path.exists ( ConfigCard ) :
+        print ( f'Reading Config Card : {ConfigCard}' )
+        ## Creates parser for reading .ini input file
+        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
+        MyReader.optionxform = str # To keep capitals
+        
+        MyReader.read (ConfigCard)
+
+        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
+            if VarName in MyReader['UserChoices'].keys() :
+                locals()[VarName] = MyReader['UserChoices'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+                
+        for VarName in ['PackFrequency'] :
+            if VarName in MyReader['Post'].keys() :
+                locals()[VarName] = MyReader['Post'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+    else :
+        raise FileExistsError ( f"File not found : {ConfigCard = }" )
+
+    
 ## Reading config file
-for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] :
+## -------------------
+for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
    if Section in config.keys () : 
-        print ( f'Reading [{Section}]' )
+        print ( f'\nReading [{Section}]' )
         for VarName in config[Section].keys() :
             locals()[VarName] = config[Section][VarName]
-            exec ( f'{VarName} = wu.setBool ({VarName})' )
-            exec ( f'{VarName} = wu.setNum  ({VarName})' )
-            exec ( f'{VarName} = wu.setNone ({VarName})' )
-            exec (  f'wu.{VarName} = {VarName}' )
-            print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) )
+            exec  ( f'{VarName} = wu.setBool ({VarName})' )
+            exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+            exec  ( f'{VarName} = wu.setNone ({VarName})' )
+            exec  ( f'wu.{VarName} = {VarName}' )
+            print ( f'    {VarName:21} set to : {locals()[VarName]}' )
             #exec ( f'del {VarName}' )
 
-#-- Some physical constants
-if wu.unDefined ( 'Ra' )           : Ra          = 6366197.7236758135 #-- Earth Radius (m)
-if wu.unDefined ( 'Grav' )         : Grav        = 9.81               #-- Gravity (m^2/s
-if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = 917.0              #-- Ice volumic mass (kg/m3) in LIM3
-if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3
-if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = 1026.0             #-- Ocean water volumic mass (kg/m3) in NEMO
-if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = 1000.0             #-- Water volumic mass in atmosphere (kg/m^3)
-if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = 1000.0             #-- Water volumic mass in surface reservoir (kg/m^3)
-if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = 1000.0             #-- Water volumic mass of rivers (kg/m^3)
-if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = 1000.              #-- Water volumic mass in ice ponds (kg/m^3)
-if wu.unDefined ( 'YearLength' )   : YearLength  = 365.25 * 24. * 60. * 60. #-- Year length (s) - Use only to compute drif in approximate unit 
-
-# Set libIGCM and machine dependant values
-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}
+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
+## =======================
+if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m)
+if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s
+if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO
+if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3)
+if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3)
+if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3)
+if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3)
+if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s)
+
+## Set libIGCM and machine dependant values
+## ----------------------------------------
+if not 'Files' in config.keys () : config['Files'] = {}
+
+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':str(ContinueOnError), 'SortIco':str(SortIco), 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) }
 
 ## --------------------------
@@ -97 +197 @@
 LMDZ = ( 'LMD' in wu.ATM )
 
-with open ('SetLibIGCM.py') as f: exec ( f.read() )
+mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
+             ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None, 
+             R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None )
+globals().update(mm)
+
 config['Files']['TmpDir'] = TmpDir
-
-if libIGCM :
-    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
-
-# Import specific module
-import nemo, lmdz
-## Now import needed scientific modules
-import xarray as xr
-    
-# Output file with water budget diagnostics
-if wu.unDefined ( 'FileOut' ) : FileOut = f'CPL_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
+config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild }
+    
+## Defines begining and end of experiment
+## --------------------------------------
+if wu.unDefined ( 'DateBegin' ) :
+    DateBegin = f'{YearBegin}0101'
+    config['Experiment']['DateBegin'] = str(DateBegin)
+else :
+    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin )
+    DateBegin = wu.FormatToGregorian (DateBegin)
+    config['Experiment']['YearBegin'] = str(YearBegin)
+
+if wu.unDefined ( 'DateEnd' )   :
+    DateEnd   = f'{YearEnd}1231'
+    config['Experiment']['DateEnd'] = str(DateEnd)
+else :
+    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd )
+    DateEnd   = wu.FormatToGregorian (DateEnd)
+    config['Experiment']['DateEnd'] = str(DateEnd)
+
+if wu.unDefined ( 'PackFrequency' ) :
+    PackFrequency = YearEnd - YearBegin + 1
+    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
+
+if type ( PackFrequency ) == str : 
+    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
+    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
+    PackFrequency = int ( PackFrequency )
+    
+## Output file with water budget diagnostics
+## -----------------------------------------
+if wu.unDefined ( 'FileOut' ) :
+    FileOut = f'CPL_waterbudget_{JobName}_{YearBegin}_{YearEnd}'
+    if ICO : 
+        if ATM_HIS == 'latlon' : FileOut = f'{FileOut}_LATLON'
+        if ATM_HIS == 'ico'    : FileOut = f'{FileOut}_ICO'
+    if readPrec == np.float32  : FileOut = f'{FileOut}_float32'
+    FileOut = f'{FileOut}.out'
+
 config['Files']['FileOut'] = FileOut
 
 f_out = open ( FileOut, mode = 'w' )
     
-# Useful functions
+## Useful functions
+## ----------------
+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'''
@@ -124 +262 @@
 
 def var2prt (var, small=False, rho=ATM_rho) :
-    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
+    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=ATM_rho, width=15) :
     if small :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
     else :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
-    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) )
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
+    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) )
     return None
 
@@ -144 +282 @@
     f_out.flush ()
     return None
-    
+
+echo ( f'{ContinueOnError = }' )
+echo ( f'{SortIco         = }' )
+echo ( f'{readPrec        = }' )
+
+echo ( f'{JobName         = }' )
+echo ( f'{ConfigCard      = }' )
+echo ( f'{libIGCM         = }' )     
+echo ( f'{User            = }' )       
+echo ( f'{Group           = }' )       
+echo ( f'{Freq            = }' )       
+echo ( f'{YearBegin       = }' )     
+echo ( f'{YearEnd         = }' )     
+echo ( f'{DateBegin       = }' )
+echo ( f'{DateEnd         = }' )
+echo ( f'{PackFrequency   = }' ) 
+echo ( f'{ATM             = }' )       
+echo ( f'{Routing         = }' )       
+echo ( f'{ORCA            = }' )      
+echo ( f'{NEMO            = }' )      
+echo ( f'{Coupled         = }' )      
+echo ( f'{ATM_HIS         = }' )      
+echo ( f'{SRF_HIS         = }' )      
+echo ( f'{RUN_HIS         = }' )
+
 ## Set libIGCM directories
+## -----------------------
 if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
 if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
@@ -156 +319 @@
 if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
 
-config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR }
-
-# Set directory to extract filesa
-if RunDir == None : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
-config['Files']['RunDir'] = RunDir
-
-if not os.path.isdir ( RunDir ) : os.makedirs ( RunDir )
-
-# Set directories to rebuild ocean and ice restart files
-RunDirOCE = os.path.join ( RunDir, 'OCE' )
-RunDirICE = os.path.join ( RunDir, 'ICE' )
-#RunDirATM = os.path.join ( RunDir, 'ATM' )
-#RunDirSRF = os.path.join ( RunDir, 'SRF' )
-#RunDirRUN = os.path.join ( RunDir, 'RUN' )
-#RunDirDYN = os.path.join ( RunDir, 'DYN' )
-if not os.path.exists ( RunDirOCE ) : os.mkdir ( RunDirOCE )
-if not os.path.exists ( RunDirICE ) : os.mkdir ( RunDirICE )
+config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
+                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
+
+## Set directory to extract files
+## ------------------------------
+if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
+config['Files']['FileDir'] = FileDir
+
+if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
+
+##- Set directories to rebuild ocean and ice restart files
+if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
+if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
+if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
+if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
+
+echo (' ')
+echo ( f'JobName     : {JobName}'    )
+echo ( f'Comment     : {Comment}'    )
+echo ( f'TmpDir      : {TmpDir}'     )
+echo ( f'FileDir     : {FileDir}'    )
+echo ( f'FileDirOCE  : {FileDirOCE}' )
+echo ( f'FileDirICE  : {FileDirICE}' )
+
+echo ( f'\nDealing with {L_EXP}'  )
 
 echo (' ')
@@ -178 +349 @@
 echo ( f'Comment   : {Comment}'   )
 echo ( f'TmpDir    : {TmpDir}'    )
-echo ( f'RunDirOCE : {RunDirOCE}' )
-echo ( f'RunDirICE : {RunDirICE}' )
 
 echo ( f'\nDealing with {L_EXP}'  )
 
-#-- Model output directories
+## Creates model output directory names
+## ------------------------------------
 if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' )
 if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' )
@@ -205 +375 @@
 echo ( f'{dir_SRF_his}' )
 
-#-- Creates files names
-if Period == None :
-    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
-    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
+##-- Creates files names
+if wu.unDefined ( 'Period' ) :
+    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
+    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
     config['Files']['Period'] = Period
-if FileCommon == None :
+
+config['Files']['DateBegin'] = DateBegin
+config['Files']['DateBegin'] = DateEnd
+
+echo ( f'Period   : {Period}' )
+
+if wu.unDefined ( 'FileCommon' ) :
     FileCommon = f'{JobName}_{Period}'
     config['Files']['FileCommon'] = FileCommon
 
-if Title == None :
-    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
+if wu.unDefined ( 'Title' ) :
+    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
     config['Files']['Title'] = Title
-
 echo ('\nOpen history files' )
-if file_ATM_his == None : 
-    file_ATM_his = os.path.join (  dir_ATM_his, f'{FileCommon}_histmth.nc' )
+if wu.unDefined ( 'file_ATM_his' ) :
+    if ATM_HIS == 'latlon' :
+        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' )
+    if ATM_HIS == 'ico' :
+        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' )
     config['Files']['file_ATM_his'] = file_ATM_his
-if file_SRF_his == None :
-    file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+if wu.unDefined ( 'file_SRF_his' ) :
+    if ATM_HIS == 'latlon' :
+        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+    if ATM_HIS == 'ico' :
+        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
     config['Files']['file_SRF_his'] = file_SRF_his
-#if Routing == 'SECHIBA'  :
-#     file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+
 if Routing == 'SIMPLE' :
-    if file_RUN_his == None : 
-        file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+    if file_RUN_his == None :
+        if ATM_HIS == 'latlon' :
+            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+        if ATM_HIS == 'ico' :
+            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
         config['Files']['file_RUN_his'] = file_RUN_his
-if file_OCE_his == None :
-    file_OCE_his = os.path.join (  dir_OCE_his, f'{FileCommon}_grid_T.nc' )
-    config['Files']['file_OCE_his'] = file_OCE_his
-if file_OCE_sca == None :
-    file_OCE_sca = os.path.join (  dir_OCE_his, f'{FileCommon}_scalar.nc' )
-    config['Files']['file_OCE_sca'] = file_OCE_sca
-if file_ICE_his == None :
-    file_ICE_his = os.path.join (  dir_ICE_his, f'{FileCommon}_icemod.nc' )
-    config['Files']['file_ICE_his'] = file_ICE_his
-if file_OCE_srf == None :
-    file_OCE_srf = os.path.join (  dir_OCE_his, f'{FileCommon}_grid_T.nc' )
-    config['Files']['file_OCE_srf'] = file_OCE_srf
-
-d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-if NEMO == '3.6' :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
-d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
-if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
-
-    
+
 echo ( f'{file_ATM_his = }' )
 echo ( f'{file_SRF_his = }' )
 if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
+        
+d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
+if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+    
+if wu.unDefined ('file_OCE_his' ) :
+    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' )
+    file_OCE_his = file_OCE_his
+if wu.unDefined ('file_OCE_sca' ) :    
+    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' )
+    config['Files']['file_OCE_sca'] = file_OCE_sca
+if wu.unDefined ('file_OCE_srf' ) :    
+    file_OCE_srf = os.path.join ( dir_OCE_his, f'{FileCommon}_sbc.nc' )
+    config['Files']['file_OCE_srf'] = file_OCE_srf
+if wu.unDefined ( 'file_ICE_hi' ) : 
+    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' )
+    config['Files']['file_ICE_his'] = file_ICE_his
+
+d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+#d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
+if NEMO == '3.6' :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
+
 echo ( f'{file_OCE_his = }' )
 echo ( f'{file_ICE_his = }' )
@@ -264 +449 @@
 
 ## Compute run length
+## ------------------
 dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
 echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
@@ -278 +464 @@
 NbYear = dtime_sec / YearLength
 
-#-- Open restart files
-YearRes = YearBegin - 1        # Year of the restart of beginning of simulation
-YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-
-config['Files']['YearPre'] = f'{YearBegin}'
-
-echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-
-if TarRestartPeriod_beg == None : 
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231'
-    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
-
-if TarRestartPeriod_end == None : 
-    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231'
-    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
-
-if tar_restart_beg == None :
-    tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' )
-    config['Files']['tar_restart_beg'] = tar_restart_beg
-if tar_restart_end == None :
-    tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' )
-    config['Files']['tar_restart_end'] = tar_restart_end
-
-echo ( f'{tar_restart_beg}' )
-echo ( f'{tar_restart_end}' )
-
-if file_ATM_beg == None : 
-    file_ATM_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc'
-    config['Files']['file_ATM_beg'] = file_ATM_beg
-if file_ATM_end == None :
-    file_ATM_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc'
-    config['Files']['file_ATM_end'] = file_ATM_end
-
-if file_DYN_beg == None :
-    if LMDZ : file_DYN_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restart.nc'
-    if ICO  : file_DYN_beg = f'{RunDir}/ICO_{JobName}_{YearRes}1231_restart.nc'
-    config['Files']['file_DYN_beg'] = file_DYN_beg
-    
-if file_DYN_end == None : 
-    if LMDZ : file_DYN_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restart.nc'
-    if ICO  : file_DYN_end = f'{RunDir}/ICO_{JobName}_{YearEnd}1231_restart.nc'
-    config['Files']['file_DYN_end'] = file_DYN_end
-    
-if file_SRF_beg == None :
-    file_SRF_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'
-    config['Files']['file_SRF_beg'] = file_SRF_beg
-if file_SRF_end == None :
-    file_SRF_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc'
-    config['Files']['file_SRF_end'] = file_SRF_end
-    
-if file_OCE_beg == None : 
-    file_OCE_beg = f'{RunDir}/OCE_{JobName}_{YearRes}1231_restart.nc'
-    config['Files']['file_OCE_beg'] = file_OCE_beg
-if file_OCE_end == None :
-    file_OCE_end = f'{RunDir}/OCE_{JobName}_{YearEnd}1231_restart.nc'
-    config['Files']['file_OCE_end'] = file_OCE_end
-if file_ICE_beg == None :
-    file_ICE_beg = f'{RunDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc'
-    config['Files']['file_ICE_beg'] = file_ICE_beg
-if file_ICE_end == None : 
-    file_ICE_end = f'{RunDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc'
-    config['Files']['file_ICE_end'] = file_ICE_end
-
-liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg]
-liste_end = [file_ATM_end, file_DYN_end, file_SRF_end]
-
-if Routing == 'SIMPLE' :
-    if file_RUN_beg == None : 
-        file_RUN_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc'
-        config['Files']['file_RUN_beg'] = file_RUN_beg
-    if file_RUN_end == None : 
-        file_RUN_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc'
-        config['Files']['file_RUN_end'] = file_RUN_end
-        
-    liste_beg.append ( file_RUN_beg )
-    liste_end.append ( file_RUN_end )
-    echo ( f'{file_RUN_beg = }' )
-    echo ( f'{file_RUN_end = }' )
-
-echo ( f'{file_ATM_beg = }' )
-echo ( f'{file_ATM_end = }' )
-echo ( f'{file_DYN_beg = }' )
-echo ( f'{file_DYN_end = }' )
-echo ( f'{file_SRF_beg = }' )
-echo ( f'{file_SRF_end = }' )
-echo ( f'{file_RUN_beg = }' )
-echo ( f'{file_RUN_end = }' )
-echo ( f'{file_OCE_beg = }' )
-echo ( f'{file_OCE_end = }' )
-echo ( f'{file_ICE_beg = }' )
-echo ( f'{file_ICE_end = }' )
-
-echo ('\nExtract restart files from tar : ATM, ICO and SRF')
-for resFile in liste_beg + liste_end :
-    if os.path.exists ( os.path.join (RunDir, resFile) ) :
-        echo ( f'file found : {resFile = }' )
-    else :
-        base_file = Path (file_name).stem # basename, and remove suffix
-        command =  f'cd {RunDir} ; tar xf {tar_restart_beg} {base_resFile}.nc'
-        echo ( command )
-        ierr = os.system ( command )
-        if ierr == 0 : echo ( f'tar done : {base_resFile}')
-        else         : raise Exception ( f'command failed : {command}' )
-
-echo ('\nOpening ATM SRF and ICO restart files')
-d_ATM_beg = xr.open_dataset ( os.path.join (RunDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze()
-d_ATM_end = xr.open_dataset ( os.path.join (RunDir, file_ATM_end), decode_times=False, decode_cf=True).squeeze()
-d_SRF_beg = xr.open_dataset ( os.path.join (RunDir, file_SRF_beg), decode_times=False, decode_cf=True).squeeze()
-d_SRF_end = xr.open_dataset ( os.path.join (RunDir, file_SRF_end), decode_times=False, decode_cf=True).squeeze()
-d_DYN_beg = xr.open_dataset ( os.path.join (RunDir, file_DYN_beg), decode_times=False, decode_cf=True).squeeze()
-d_DYN_end = xr.open_dataset ( os.path.join (RunDir, file_DYN_end), decode_times=False, decode_cf=True).squeeze()
-
-echo ('\nExtract and rebuild OCE and ICE restarts')
-def get_ndomain (zfile) :
-    #d_zfile = xr.open_dataset (zfile, decode_times=False).squeeze()
-    #ndomain_opa = d_zfile.attrs['DOMAIN_number_total']
-    #d_zfile.close ()
-    ndomain_opa = subprocess.getoutput ( f'ls {zfile}_*.nc | wc -l' ) #.format()
-    return int (ndomain_opa)
-
-def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, RunDirComp=RunDirOCE ) :
-    '''Extract restart file from tar. Rebuild ocean files if needed'''
-    echo ( f'file to extract : {file_name} ' )
-    if os.path.exists ( file_name ) :
-        echo ( f'-- File ready : {file_name}' )
-    else : 
-        echo ( f'-- Extracting {file_name}' )
-        base_resFile = Path (file_name).stem # basename, and remove suffix
-        # Try to extract the rebuilded file
-        command =  f'cd {RunDirComp} ; tar xf {tar_restart} {base_resFile}.nc'
-        echo ( command )
-        ierr = os.system ( command )
-        if ierr == 0 :
-            echo ( f'tar done : {base_resFile}')
-            command = f'cd {RunDirComp} ; mv {base_resFile}.nc {RunDir}'
-            ierr = os.system ( command )
-            if ierr == 0 : echo ( f'command done : {command}' )
-            else         : raise Exception ( f'command failed : {command}' )
-        else :
-            if not os.path.exists ( os.path.join (RunDir, f'{base_resFile}_0000.nc') ):
-                command =  f'cd {RunDirComp} ; tar xf {tar_restart_end} {base_resFile}_*.nc'
-                echo ( command )
-                ierr = os.system ( command )
-                if ierr == 0 : echo ( f'tar done : {file_OCE_beg}')
-                else         : raise Exception ( f'command failed : {command}' )
-                echo ('extract ndomain' )
-            ndomain_opa = get_ndomain ( os.path.join (RunDir, f'{base_resFile}') )
-            command = f'cd {RunDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {RunDir}'
-            echo ( command )
-            ierr = os.system ( command )
-            if ierr == 0 : echo ( f'Rebuild done : {base_resFile}.nc')
-            else         : raise Exception ( f'command failed : {command}' )
-
-extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirOCE )
-extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, RunDirComp=RunDirOCE )
-extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirICE )
-extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, RunDirComp=RunDirICE )
-
-echo ('Opening OCE and ICE restart files')
-if NEMO == 3.6 : 
-    d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()
-    d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()
-    d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()
-if NEMO == 4.0 or NEMO == 4.2 : 
-    d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-
 ## Write the full configuration
 config_out = open (FullIniFile, 'w')
@@ -456 +469 @@
 config_out.close ()
 
-for var in d_SRF_beg.variables :
-    d_SRF_beg[var] = d_SRF_beg[var].where (  d_SRF_beg[var]<1.e20, 0.)
-    d_SRF_end[var] = d_SRF_end[var].where (  d_SRF_end[var]<1.e20, 0.)
-
-if ICO :
-    d_RUN_beg = xr.open_dataset ( os.path.join (RunDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze()
-    d_RUN_end = xr.open_dataset ( os.path.join (RunDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze()
-
-def kg2Sv  (val, rho=OCE_rho_liq) :
-    '''From kg to Sverdrup'''
-    return val/dtime_sec*1.0e-6/rho
-
-def kg2myear (val, rho=OCE_rho_liq) :
-    '''From kg to m/year'''
-    return val/OCE_aire_tot/rho/NbYear
-
-def var2prt (var, small=False) :
-    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000.
-    else     :  return var , kg2Sv(var)      , kg2myear(var)
-
-def prtFlux (Desc, var, Form='F', small=False) :
-    if small :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
-    else :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
-    echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) )
-    return None
-
-
-
 # ATM grid with cell surfaces
+if LMDZ :
+    echo ('ATM grid with cell surfaces : LMDZ')
+    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
+    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
+    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire'] )    [0], cumulPoles=True, dim1D='cell' )
+    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
+    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
+    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
+    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
+    SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
+    SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
+    SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1D='cell', cumulPoles=True )
+    SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1D='cell', cumulPoles=True )
+    SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1D='cell' )
 if ICO :
     if ATM_HIS == 'latlon' :
-        jpja, jpia = d_ATM_his['aire'][0].shape   
-        file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' )
-        config['Files']['file_ATM_aire'] = file_ATM_aire
-        echo ( f'Aire sur grille reguliere : {file_ATM_aire = }' )
-        d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze()
-        ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True )
-        ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] )
-        ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] )
-        ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] )
-        ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] )
-        SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] )
-        #SRF_aire = ATM_aire * lmdz.geo2point (d_SRF_his ['Contfrac'] )
-        #SRF_aire = ATM_aire * ATM_fter
+        echo ( 'ATM areas and fractions on latlon grid' )
+        if 'lat_dom_out' in d_ATM_his.variables :
+            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
+            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
+        else :
+            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
+            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
+        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumulPoles=True, dim1D='cell' )
+        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
+        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
+        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
+        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
+
     if ATM_HIS == 'ico' :
-        echo ( f'Aire sur grille icosaedre : {file_ATM_aire = }' )
-        
-if LMDZ :
-    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=True )
-    ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] )
-    ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] )
-    ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] )
-    ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] )
-    #SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] )
-    SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'] )
-    
+        echo ( 'ATM areas and fractions on ICO grid' )
+        ATM_aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
+        ATM_lat  =  rprec (d_ATM_his ['lat']         )
+        ATM_lon  =  rprec (d_ATM_his ['lon']         )
+        ATM_fter =  rprec (d_ATM_his ['fract_ter'][0])
+        ATM_foce =  rprec (d_ATM_his ['fract_oce'][0])
+        ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0])
+        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0])
+
+    if SRF_HIS == 'latlon' :
+        echo ( 'SRF areas and fractions on latlon grid' )
+        if 'lat_domain_landpoints_out' in d_SRF_his  : 
+            SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
+            SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
+        else : 
+            if 'lat_domain_landpoints_out' in d_SRF_his :
+                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
+                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
+            else :
+                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
+                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
+                
+        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1D='cell', cumulPoles=True )
+        SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1D='cell', cumulPoles=True )
+        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1D='cell', cumulPoles=True )
+        SRF_aire = SRF_areafrac
+ 
+    if SRF_HIS == 'ico' :
+        echo ( 'SRF areas and fractions on latlon grid' )
+        SRF_lat       =  rprec (d_SRF_his ['lat']     )
+        SRF_lon       =  rprec (d_SRF_his ['lon']     )
+        SRF_areas     =  rprec (d_SRF_his ['Areas']   ) 
+        SRF_contfrac  =  rprec (d_SRF_his ['Contfrac'])
+        SRF_aire      =  SRF_areas * SRF_contfrac
 ATM_fsea      = ATM_foce + ATM_fsic
 ATM_flnd      = ATM_fter + ATM_flic
@@ -525 +543 @@
 ATM_aire_fsea = ATM_aire * ATM_fsea
 
-SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.)
-
-if ICO :
-    # Area on icosahedron grid
-    file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
-    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze()
-    d_DYN_aire = d_DYN_aire.rename ( {'cell':'cell_mesh'} )
-    DYN_aire   = d_DYN_aire['aire']
-
-    DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic']
-    DYN_flnd = 1.0 - DYN_fsea
-    
-if LMDZ :
-    # Area on lon/lat grid
-    DYN_aire = ATM_aire
-    DYN_fsea = ATM_fsea
-    DYN_flnd = ATM_flnd
-    DYN_fter = d_ATM_beg['FTER']
-    DYN_flic = d_ATM_beg['FTER']
-
-def ATM_stock_int (stock) :
-    '''Integrate (* surface) stock on atmosphere grid'''
-    ATM_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
-    return ATM_stock_int
+#SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.)
+
+# if ICO :
+#     if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
+#     config['Files']['file_DYN_aire'] = file_DYN_aire
+
+# if ICO :
+#     # Area on icosahedron grid
+#     d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
+           
+#     DYN_lat = d_DYN_aire['lat']
+#     DYN_lon = d_DYN_aire['lon']
+
+#     DYN_aire = d_DYN_aire['aire']
+#     DYN_fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
+#     DYN_flnd = 1.0 - DYN_fsea
+#     DYN_fter = d_ATM_beg['FTER']
+#     DYN_flic = d_ATM_beg['FLIC']
+#     DYN_aire_fter = DYN_aire * DYN_fter
+    
+# if LMDZ :
+#     # Area on lon/lat grid
+#     DYN_aire = ATM_aire
+#     DYN_fsea = ATM_fsea
+#     DYN_flnd = ATM_flnd
+#     DYN_fter = rprec (d_ATM_beg['FTER'])
+#     DYN_flic = rprec (d_ATM_beg['FLIC'])
+#     DYN_aire_fter = DYN_aire * DYN_fter
+
+# Functions computing integrals and sum
+# def ATM_stock_int (stock) :
+#     '''Integrate (* surface) stock on atmosphere grid'''
+#     ATM_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
+#     return ATM_stock_int
 
 def ATM_flux_int (flux) :
@@ -555 +583 @@
     return ATM_flux_int
 
-def SRF_stock_int (stock) :
-    '''Integrate (* surface) stock on land grid'''
-    ATM_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
-    return ATM_stock_int
+# def SRF_stock_int (stock) :
+#     '''Integrate (* surface) stock on land grid'''
+#     ATM_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
+#     return ATM_stock_int
 
 def SRF_flux_int (flux) :
@@ -565 +593 @@
     return SRF_flux_int
 
-def OCE_stock_int (stock) :
-    '''Integrate stock on ocean grid'''
-    OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) )
-    return OCE_stock_int
+def LIC_flux_int (flux) :
+    '''Integrate (* time * surface) flux on land ice grid'''
+    LIC_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() )
+    return LIC_flux_int
+
+# def OCE_stock_int (stock) :
+#     '''Integrate stock on ocean grid'''
+#     OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) )
+#     return OCE_stock_int
 
 def ONE_stock_int (stock) :
@@ -584 +617 @@
     OCE_flux_int = np.sum (  np.sort ( (flux * dtime_per_sec).to_masked_array().ravel()) )
     return OCE_flux_int
-    
-#if LMDZ :
-#    d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} )
 
 # Get mask and surfaces
@@ -607 +637 @@
 
 echo ( '\n====================================================================================' )
-echo ( '-- NEMO change in stores (for the records)' )
-#-- Note that the total number of days of the run should be diagnosed rather than assumed
-#-- Here the result is in Sv
-#
-#-- Change in ocean volume in freshwater equivalent
-
-OCE_ssh_beg = d_OCE_beg['sshn']
-OCE_ssh_end = d_OCE_end['sshn']
-OCE_sum_ssh_beg = OCE_stock_int ( OCE_ssh_beg )
-OCE_sum_ssh_end = OCE_stock_int ( OCE_ssh_end )
-
-OCE_mas_wat_beg = OCE_sum_ssh_beg * OCE_rho_liq
-OCE_mas_wat_end = OCE_sum_ssh_end * OCE_rho_liq
-
-echo ( '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) )
-dOCE_vol_liq = ( OCE_sum_ssh_end - OCE_sum_ssh_beg )
-dOCE_mas_liq = dOCE_vol_liq * OCE_rho_liq
-dOCE_mas_wat = dOCE_mas_liq
-
-echo ( 'dOCE vol    = {:12.3e} m^3'.format (dOCE_vol_liq) )
-echo ( 'dOCE ssh    = {:12.3e} m  '.format (dOCE_vol_liq/OCE_aire_tot) )
-echo ( 'dOCE mass   = {:12.3e} kg '.format (dOCE_mas_liq) )
-echo ( 'dOCE mass   = {:12.3e} Sv '.format (dOCE_mas_liq/dtime_sec*1E-6/OCE_rho_liq) )
-
-## Glace et neige
+echo ( f'-- ATM Fluxes  -- {Title} ' )
+
+if ATM_HIS == 'latlon' :
+    echo ( ' latlon case' )
+    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' )
+    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' )
+    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' )
+    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' )
+    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
+    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' )
+    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' )
+    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' )
+    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' )
+    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' )
+    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' )
+    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' )
+    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' )
+    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' )
+    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' )
+    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' )
+    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' )
+    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' )
+    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' )
+    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' )
+    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' )
+    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' )
+    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
+    echo ( f'End of LATLON case')
+    
+if ATM_HIS == 'ico' :
+    echo (' ico case')
+    ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
+    ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
+    ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
+    ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
+    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
+    ATM_fqcalving   = rprec (d_ATM_his ['fqcalving'])
+    ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  )
+    ATM_precip      = rprec (d_ATM_his ['precip']   )
+    ATM_snowf       = rprec (d_ATM_his ['snow']     )
+    ATM_evap        = rprec (d_ATM_his ['evap']     )
+    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
+    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
+    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
+    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
+    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
+    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
+    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
+    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
+    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
+    ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
+    ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
+    ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
+    ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
+    ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
+    ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
+    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
+    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
+    echo ( f'End of ico case ')
+
+
+    
+echo ( 'ATM wprecip_oce' )
+ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
+ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
+ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
+ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
+
+ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic
+ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic
+ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic
+ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic
+ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic
+ATM_wemp        = ATM_wevap - ATM_wprecip
+ATM_emp         = ATM_evap  - ATM_precip
+
+ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
+ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
+ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
+ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
+ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
+
+ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter
+ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce
+ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic
+ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic
+ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce
+
+if RUN_HIS == 'latlon' :
+    echo ( f'RUN costalflow Grille LATLON' )
+    if TestInterp :
+        echo ( f'RUN runoff TestInterp' )
+        RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' )
+        RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' )
+    else : 
+        echo ( f'RUN runoff' )
+        RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' )
+        RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' )
+        
+    RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' ) 
+    RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' )
+    RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' )
+    RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' ) 
+    RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' )
+
+if RUN_HIS == 'ico' :
+    echo ( f'RUN costalflow Grille ICO' )
+    RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])
+    RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )
+    RUN_runoff      =  rprec (d_RUN_his ['runoff']     )
+    RUN_drainage    =  rprec (d_RUN_his ['drainage']   )
+    RUN_riversret   =  rprec (d_RUN_his ['riversret']  )
+    
+    RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
+    RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
+
+## Correcting units of SECHIBA variables
+def mmd2SI ( Var ) :
+    '''Change unit from mm/d or m^3/s to kg/s if needed'''
+    if 'units' in VarT.attrs : 
+        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
+            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s'
+        if VarT.attrs['units'] == 'mm/d' :
+            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
+        if VarT.attrs['units'] in ['m^3', 'm3', ] :
+            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg'
+            
+for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
+    VarT = locals()['RUN_' + var]
+    mmd2SI (VarT)
+
+#for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
+#    VarT = locals()['SRF_' + var]
+#    mmd2SI (VarT)
+echo ( f'RUN input' )
+RUN_input  = RUN_runoff      + RUN_drainage
+RUN_output = RUN_coastalflow + RUN_riverflow
+
+echo ( f'ATM flw_wbilo' )
+ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
+ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      )
+ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    )
+ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      )
+ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      )
+ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       )
+
+ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
+ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
+ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
+ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
+ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
+ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
+ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
+
+LIC_flx_calving     = LIC_flux_int ( ATM_fqcalving  )
+LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    )
+
+ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
+ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
+ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
+ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
+
+LIC_flx_precip      = LIC_flux_int ( ATM_precip     )
+LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      )
+LIC_flx_evap        = LIC_flux_int ( ATM_evap       )
+LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  )
+
+ATM_flx_wrain_ter   = ATM_flux_int ( ATM_wrain_ter )
+ATM_flx_wrain_oce   = ATM_flux_int ( ATM_wrain_oce )
+ATM_flx_wrain_lic   = ATM_flux_int ( ATM_wrain_lic )
+ATM_flx_wrain_sic   = ATM_flux_int ( ATM_wrain_sic )
+ATM_flx_wrain_sea   = ATM_flux_int ( ATM_wrain_sea )
+
+ATM_flx_wsnow_ter   = ATM_flux_int ( ATM_wsnow_ter )
+ATM_flx_wsnow_oce   = ATM_flux_int ( ATM_wsnow_oce )
+ATM_flx_wsnow_lic   = ATM_flux_int ( ATM_wsnow_lic )
+ATM_flx_wsnow_sic   = ATM_flux_int ( ATM_wsnow_sic )
+ATM_flx_wsnow_sea   = ATM_flux_int ( ATM_wsnow_sea )
+
+ATM_flx_wevap_ter   = ATM_flux_int ( ATM_wevap_ter )
+ATM_flx_wevap_oce   = ATM_flux_int ( ATM_wevap_oce )
+ATM_flx_wevap_lic   = ATM_flux_int ( ATM_wevap_lic )
+ATM_flx_wevap_sic   = ATM_flux_int ( ATM_wevap_sic )
+ATM_flx_wevap_sea   = ATM_flux_int ( ATM_wevap_sea )
+ATM_flx_wprecip_lic = ATM_flux_int ( ATM_wprecip_lic )
+ATM_flx_wprecip_oce = ATM_flux_int ( ATM_wprecip_oce )
+ATM_flx_wprecip_sic = ATM_flux_int ( ATM_wprecip_sic )
+ATM_flx_wprecip_ter = ATM_flux_int ( ATM_wprecip_ter )
+ATM_flx_wprecip_sea = ATM_flux_int ( ATM_wprecip_sea )
+ATM_flx_wemp_lic    = ATM_flux_int ( ATM_wemp_lic )
+ATM_flx_wemp_oce    = ATM_flux_int ( ATM_wemp_oce )
+ATM_flx_wemp_sic    = ATM_flux_int ( ATM_wemp_sic )
+ATM_flx_wemp_ter    = ATM_flux_int ( ATM_wemp_ter )
+ATM_flx_wemp_sea    = ATM_flux_int ( ATM_wemp_sea )
+
+ATM_flx_emp          = ATM_flux_int ( ATM_emp )
+
+RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
+RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
+RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
+RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
+RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
+RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
+RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
+RUN_flx_input       = SRF_flux_int ( RUN_input      )
+RUN_flx_output      = ONE_flux_int ( RUN_output     )
+
+RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output)
+RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow)
+
+prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )          
+prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )          
+prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )          
+prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )          
+prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )          
+prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True)  
+prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True)  
+prtFlux ('calving              ', ATM_flx_calving       , 'f' ) 
+prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )       
+prtFlux ('precip               ', ATM_flx_precip        , 'f' )       
+prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )        
+prtFlux ('evap                 ', ATM_flx_evap          , 'f' )
+prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' )
+
+prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' )
+prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' )
+prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' )
+prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' )
+prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True )
+
+prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True )
+prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True )
+prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True )
+prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True )
+
+
+echo ( '\n====================================================================================' )
+echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
+prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
+prtFlux ('riverflow     ', RUN_flx_river      , 'f' )        
+prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' )  
+prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )    
+prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
+prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
+prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
+prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
+prtFlux ('river in      ', RUN_flx_input      , 'f' )   
+prtFlux ('river out     ', RUN_flx_output     , 'f' )   
+prtFlux ('river bil     ', RUN_flx_bil        , 'f' ) 
+    
+echo ( '\n====================================================================================' )
+echo ( f'-- OCE Fluxes  -- {Title} ' )
+
+# Read variable and computes integral over space and time
+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     = 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 :
-    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']
-    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']
-
-    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']
-    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']
-    
-    ICE_pnd_beg = 0.0 ; ICE_pnd_end = 0.0
-    ICE_fzl_beg = 0.0 ; ICE_fzl_end = 0.0
-
-    ICE_mas_wat_beg = OCE_stock_int ( (ICE_ice_beg*ICE_rho_ice + ICE_sno_beg*ICE_rho_sno) )
-    ICE_mas_wat_end = OCE_stock_int ( (ICE_ice_end*ICE_rho_ice + ICE_sno_end*ICE_rho_sno) )
-    
+    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   = 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_ice_beg = d_ICE_beg ['v_i']  ; ICE_ice_end = d_ICE_end ['v_i']
-    ICE_sno_beg = d_ICE_beg ['v_s']  ; ICE_sno_end = d_ICE_end ['v_s']
-    ICE_pnd_beg = d_ICE_beg ['v_ip'] ; ICE_pnd_end = d_ICE_end ['v_ip']
-    ICE_fzl_beg = d_ICE_beg ['v_il'] ; ICE_fzl_end = d_ICE_end ['v_il']
-    
-    ICE_mas_wat_beg = OCE_stock_int ( d_ICE_beg['snwice_mass'] )
-    ICE_mas_wat_end = OCE_stock_int ( d_ICE_end['snwice_mass'] )
-    
-    
-ICE_vol_ice_beg = OCE_stock_int ( ICE_ice_beg )
-ICE_vol_ice_end = OCE_stock_int ( ICE_ice_end )
-
-ICE_vol_sno_beg = OCE_stock_int ( ICE_sno_beg )
-ICE_vol_sno_end = OCE_stock_int ( ICE_sno_end )
-
-ICE_vol_pnd_beg = OCE_stock_int ( ICE_pnd_beg )
-ICE_vol_pnd_end = OCE_stock_int ( ICE_pnd_end )
-
-ICE_vol_fzl_beg = OCE_stock_int ( ICE_fzl_beg )
-ICE_vol_fzl_end = OCE_stock_int ( ICE_fzl_end )
-
-#-- Converting to freswater volume
-dICE_vol_ice = ICE_vol_ice_end - ICE_vol_ice_beg
-dICE_mas_ice = dICE_vol_ice * ICE_rho_ice
-
-dICE_vol_sno = ICE_vol_sno_end - ICE_vol_sno_beg
-dICE_mas_sno = dICE_vol_sno * ICE_rho_sno
-
-dICE_vol_pnd = ICE_vol_pnd_end - ICE_vol_pnd_beg
-dICE_mas_pnd = dICE_vol_pnd * ICE_rho_pnd
-
-dICE_vol_fzl= ICE_vol_fzl_end - ICE_vol_fzl_beg
-dICE_mas_fzl = dICE_vol_fzl * ICE_rho_pnd
-
+    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 :
-    dICE_mas_wat = dICE_mas_ice + dICE_mas_sno
-    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_ice + dICE_mas_sno
-
-if NEMO == 4.0 or NEMO == 4.2 :
-    dICE_mas_wat = ICE_mas_wat_end - ICE_mas_wat_beg
-    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat
-
-echo ( f'ICE_vol_ice_beg = {ICE_vol_ice_beg:12.6e} m^3 | ICE_vol_ice_end = {ICE_vol_ice_end:12.6e} m^3' )
-echo ( f'ICE_vol_sno_beg = {ICE_vol_sno_beg:12.6e} m^3 | ICE_vol_sno_end = {ICE_vol_sno_end:12.6e} m^3' )
-echo ( f'ICE_vol_pnd_beg = {ICE_vol_pnd_beg:12.6e} m^3 | ICE_vol_pnd_end = {ICE_vol_pnd_end:12.6e} m^3' )
-echo ( f'ICE_vol_fzl_beg = {ICE_vol_fzl_beg:12.6e} m^3 | ICE_vol_fzl_end = {ICE_vol_fzl_end:12.6e} m^3' )
-
-echo ( f'dICE_vol_ice   = {dICE_vol_ice:12.3e} m^3' )
-echo ( f'dICE_vol_sno   = {dICE_vol_sno:12.3e} m^3' )
-echo ( f'dICE_vol_pnd   = {dICE_vol_pnd:12.3e} m^3' )
-echo ( f'dICE_mas_ice   = {dICE_mas_ice:12.3e} m^3' )
-echo ( f'dICE_mas_sno   = {dICE_mas_sno:12.3e} m^3' )  
-echo ( f'dICE_mas_pnd   = {dICE_mas_pnd:12.3e} m^3' )  
-echo ( f'dICE_mas_fzl   = {dICE_mas_fzl:12.3e} m^3' )  
-echo ( f'dICE_mas_wat   = {dICE_mas_wat:12.3e} m^3' ) 
-
-
-SEA_mas_wat_beg = OCE_mas_wat_beg + ICE_mas_wat_beg
-SEA_mas_wat_end = OCE_mas_wat_end + ICE_mas_wat_end
-
-echo ( '\n------------------------------------------------------------')
-echo ( 'Variation du contenu en eau ocean + glace ' )
-echo ( 'dMass (ocean)   = {:12.6e} kg '.format(dSEA_mas_wat) )
-echo ( 'dVol  (ocean)   = {:12.3e} Sv '.format(dSEA_mas_wat/dtime_sec*1E-6/OCE_rho_liq) )
-echo ( 'dVol  (ocean)   = {:12.3e} m  '.format(dSEA_mas_wat*1E-3/OCE_aire_tot) )
-
-
-echo ( '\n------------------------------------------------------------------------------------' )
-echo ( '-- ATM changes in stores ' )
-
-#-- Change in precipitable water from the atmosphere daily and monthly files
-#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
-
-# ATM vertical grid
-ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
-ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
-klevp1 = ATM_Ahyb.shape[0]
-
-# Surface pressure
-if ICO : 
-    DYN_ps_beg = d_DYN_beg['ps']
-    DYN_ps_end = d_DYN_end['ps']
-    
-if LMDZ : 
-    DYN_ps_beg =  lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) )
-    DYN_ps_end =  lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) )
-    
-# 3D Pressure
-DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg
-DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end
-
-# Layer thickness
-DYN_sigma_beg = DYN_p_beg[0:-1]*0.
-DYN_sigma_end = DYN_p_end[0:-1]*0. 
-
-for k in np.arange (klevp1-1) :
-    DYN_sigma_beg[k,:] = (DYN_p_beg[k,:] - DYN_p_beg[k+1,:]) / Grav
-    DYN_sigma_end[k,:] = (DYN_p_end[k,:] - DYN_p_end[k+1,:]) / Grav
-
-DYN_sigma_beg = DYN_sigma_beg.rename ( {'klevp1':'sigs'} )
-DYN_sigma_end = DYN_sigma_end.rename ( {'klevp1':'sigs'} )
-
-##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases
-if LMDZ :
-    try : 
-        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ) )
-        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ) )
-    except :
-        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) )
-        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) )
-if ICO :
-    try :
-        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).rename ( {'lev':'sigs'} )
-        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).rename ( {'lev':'sigs'} )
-    except : 
-        DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) ).rename ( {'lev':'sigs'} )
-        DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) ).rename ( {'lev':'sigs'} )
-    
-DYN_mas_wat_beg = ATM_stock_int (DYN_sigma_beg * DYN_wat_beg)
-DYN_mas_wat_end = ATM_stock_int (DYN_sigma_end * DYN_wat_end)
-
-dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg
-
-echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' )
-echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
-echo ( 'dMass (atm)   = {:12.3e} kg '.format (dDYN_mas_wat) )
-echo ( 'dMass (atm)   = {:12.3e} Sv '.format (dDYN_mas_wat/dtime_sec*1.e-6/ATM_rho) )
-echo ( 'dMass (atm)   = {:12.3e} m  '.format (dDYN_mas_wat/ATM_aire_sea_tot/ATM_rho) )
-
-ATM_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']
-ATM_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']
-
-ATM_qs_beg  = d_ATM_beg['QS01']*d_ATM_beg['FTER']+d_ATM_beg['QS02']*d_ATM_beg['FLIC']+d_ATM_beg['QS03']*d_ATM_beg['FOCE']+d_ATM_beg['QS04']*d_ATM_beg['FSIC']
-ATM_qs_end  = d_ATM_end['QS01']*d_ATM_end['FTER']+d_ATM_end['QS02']*d_ATM_end['FLIC']+d_ATM_end['QS03']*d_ATM_end['FOCE']+d_ATM_end['QS04']*d_ATM_end['FSIC']
-
-ATM_qsol_beg = d_ATM_beg['QSOL']
-ATM_qsol_end = d_ATM_end['QSOL']
-
-ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
-ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER']
-ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
-ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC']
-ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
-ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE']
-ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
-ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC']  
-
-if ICO :
-   ATM_sno_beg   = ATM_sno_beg .rename ( {'points_physiques':'cell_mesh'} )
-   ATM_sno_end   = ATM_sno_end .rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs_beg    = ATM_qs_beg  .rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs_end    = ATM_qs_end  .rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qsol_beg  = ATM_qsol_beg.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qsol_end  = ATM_qsol_end.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs01_beg  = ATM_qs01_beg.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs01_end  = ATM_qs01_end.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs02_beg  = ATM_qs02_beg.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs02_end  = ATM_qs02_end.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs03_beg  = ATM_qs03_beg.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs03_end  = ATM_qs03_end.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs04_beg  = ATM_qs04_beg.rename ( {'points_physiques':'cell_mesh'} )
-   ATM_qs04_end  = ATM_qs04_end.rename ( {'points_physiques':'cell_mesh'} ) 
-
-echo ('Computing atmosphere integrals')
-ATM_mas_sno_beg   = ATM_stock_int ( ATM_sno_beg  )
-ATM_mas_sno_end   = ATM_stock_int ( ATM_sno_end  )
-ATM_mas_qs_beg    = ATM_stock_int ( ATM_qs_beg   )
-ATM_mas_qs_end    = ATM_stock_int ( ATM_qs_end   )
-ATM_mas_qsol_beg  = ATM_stock_int ( ATM_qsol_beg )
-ATM_mas_qsol_end  = ATM_stock_int ( ATM_qsol_end )
-ATM_mas_qs01_beg  = ATM_stock_int ( ATM_qs01_beg )
-ATM_mas_qs01_end  = ATM_stock_int ( ATM_qs01_end )
-ATM_mas_qs02_beg  = ATM_stock_int ( ATM_qs02_beg )
-ATM_mas_qs02_end  = ATM_stock_int ( ATM_qs02_end )
-ATM_mas_qs03_beg  = ATM_stock_int ( ATM_qs03_beg )
-ATM_mas_qs03_end  = ATM_stock_int ( ATM_qs03_end )
-ATM_mas_qs04_beg  = ATM_stock_int ( ATM_qs04_beg )
-ATM_mas_qs04_end  = ATM_stock_int ( ATM_qs04_end )
-
-dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg
-dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg
-dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg
-
-dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg
-dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg
-dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg
-dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg
-
-echo ( '\nVariation du contenu en neige atmosphere (calottes)' )
-echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
-echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_sno ) )
-echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1e-6/ICE_rho_ice) )
-echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_sno/ATM_aire_sea_tot/ATM_rho) )
-
-echo ( '\nVariation du contenu humidite du sol' )
-echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) )
-echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_qs ) )
-echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_qs/dtime_sec*1e-6/ATM_rho) )
-echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_qs/ATM_aire_sea_tot/ATM_rho) )
-
-echo ( '\nVariation du contenu en eau+neige atmosphere ' )
-echo ( 'dMass (eau + neige atm) = {:12.3e} kg '.format (  dDYN_mas_wat + dATM_mas_sno) )
-echo ( 'dMass (eau + neige atm) = {:12.3e} Sv '.format ( (dDYN_mas_wat + dATM_mas_sno)/dtime_sec*1E-6/ATM_rho) )
-echo ( 'dMass (eau + neige atm) = {:12.3e} m  '.format ( (dDYN_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot/ATM_rho) )
-
-echo ( '\n------------------------------------------------------------------------------------' )
-echo ( '-- SRF changes ' )
-
-if Routing == 'SIMPLE' :
-    RUN_mas_wat_beg = ONE_stock_int ( d_RUN_beg ['fast_reservoir'] +  d_RUN_beg ['slow_reservoir'] +  d_RUN_beg ['stream_reservoir'])
-    RUN_mas_wat_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] +  d_RUN_end ['slow_reservoir'] +  d_RUN_end ['stream_reservoir'])
-
-if Routing == 'SECHIBA' : 
-    RUN_mas_wat_beg = ONE_stock_int ( d_SRF_beg['fastres']  + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \
-                                    + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] )
-    RUN_mas_wat_end = ONE_stock_int ( d_SRF_end['fastres']  + d_SRF_end['slowres'] + d_SRF_end['streamres'] \
-                                    + d_SRF_end['floodres'] + d_SRF_end['lakeres'] + d_SRF_end['pondres'] )
-
-dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg
-    
-echo ( '\nWater content in routing ' )
-echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) )
-echo ( 'dMass (routing)  = {:12.3e} kg '.format(dRUN_mas_wat) )
-echo ( 'dMass (routing)  = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) )
-echo ( 'dMass (routing)  = {:12.3e} m  '.format(dRUN_mas_wat/OCE_aire_tot*1E-3) )
-
-print ('Reading SRF restart')
-tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; tot_watveg_beg  = tot_watveg_beg .where (tot_watveg_beg < 1E10, 0.)
-tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; tot_watsoil_beg = tot_watsoil_beg.where (tot_watsoil_beg< 1E10, 0.)
-snow_beg        = d_SRF_beg['snow_beg']        ; snow_beg        = snow_beg       .where (snow_beg       < 1E10, 0.)
-
-tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; tot_watveg_end  = tot_watveg_end .where (tot_watveg_end < 1E10, 0.)
-tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; tot_watsoil_end = tot_watsoil_end.where (tot_watsoil_end< 1E10, 0.)
-snow_end        = d_SRF_end['snow_beg']        ; snow_end        = snow_end       .where (snow_end       < 1E10, 0.)
-
-if LMDZ :
-    tot_watveg_beg  = lmdz.geo2point (tot_watveg_beg)
-    tot_watsoil_beg = lmdz.geo2point (tot_watsoil_beg)
-    snow_beg        = lmdz.geo2point (snow_beg)
-    tot_watveg_end  = lmdz.geo2point (tot_watveg_end)
-    tot_watsoil_end = lmdz.geo2point (tot_watsoil_end)
-    snow_end        = lmdz.geo2point (snow_end)
-
-SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg
-SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end
-
-#SRF_mas_wat_beg = d_SRF_beg['tot_watveg_beg']+d_SRF_beg['tot_watsoil_beg'] + d_SRF_beg['snow_beg']
-#SRF_mas_wat_end = d_SRF_end['tot_watveg_beg']+d_SRF_end['tot_watsoil_beg'] + d_SRF_end['snow_beg']
-
-#if LMDZ :
-#    tot_watveg_beg  = tot_watveg_beg .rename ( {'y':'points_phyiques'} )
-#    tot_watsoil_beg = tot_watsoil_beg.rename ( {'y':'points_phyiques'} )
-#    snow_beg        = snow_beg       .rename ( {'y':'points_phyiques'} )
-#    tot_watveg_end  = tot_watveg_end .rename ( {'y':'points_phyiques'} )
-#    tot_watsoil_end = tot_watsoil_end.rename ( {'y':'points_phyiques'} )
-#    snow_end        = snow_end       .rename ( {'y':'points_phyiques'} )
-#    SRF_wat_beg     = SRF_wat_beg    .rename ( {'y':'points_phyiques'} )
-#    SRF_wat_end     = SRF_wat_end    .rename ( {'y':'points_phyiques'} )
-if ICO :
-    tot_watveg_beg  = tot_watveg_beg .rename ( {'y':'cell_mesh'} )
-    tot_watsoil_beg = tot_watsoil_beg.rename ( {'y':'cell_mesh'} )
-    snow_beg        = snow_beg       .rename ( {'y':'cell_mesh'} )
-    tot_watveg_end  = tot_watveg_end .rename ( {'y':'cell_mesh'} )
-    tot_watsoil_end = tot_watsoil_end.rename ( {'y':'cell_mesh'} )
-    snow_end        = snow_end       .rename ( {'y':'cell_mesh'} )
-    SRF_wat_beg     = SRF_wat_beg    .rename ( {'y':'cell_mesh'} )
-    SRF_wat_end     = SRF_wat_end    .rename ( {'y':'cell_mesh'} )  
-
-print ('Computing integrals')
-
-print ( ' 1/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_watveg_beg   = ATM_stock_int ( tot_watveg_beg  )
-print ( ' 2/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_watsoil_beg  = ATM_stock_int ( tot_watsoil_beg )
-print ( ' 3/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_snow_beg     = ATM_stock_int ( snow_beg        )
-print ( ' 4/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_watveg_end   = ATM_stock_int ( tot_watveg_end  )
-print ( ' 5/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_watsoil_end  = ATM_stock_int ( tot_watsoil_end )
-print ( ' 6/6', end='' ) ; sys.stdout.flush ()
-SRF_mas_snow_end     = ATM_stock_int ( snow_end        )
-print (' -- ') ; sys.stdout.flush ()
-
-dSRF_mas_watveg   = SRF_mas_watveg_end   - SRF_mas_watveg_beg
-dSRF_mas_watsoil  = SRF_mas_watsoil_end  - SRF_mas_watsoil_beg
-dSRF_mas_snow     = SRF_mas_snow_end     - SRF_mas_snow_beg
-
-echo ('\nLes differents reservoirs')
-echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' )
-echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' )
-echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' )
-
-echo ( 'dMass (watveg)  = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watveg , dSRF_mas_watveg /dtime_sec*1E-9, dSRF_mas_watveg /OCE_aire_tot*1E-3) )
-echo ( 'dMass (watsoil) = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watsoil, dSRF_mas_watsoil/dtime_sec*1E-9, dSRF_mas_watsoil/OCE_aire_tot*1E-3) )
-echo ( 'dMass (sno)     = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_snow   , dSRF_mas_snow   /dtime_sec*1E-9, dSRF_mas_snow   /OCE_aire_tot*1E-3  ) )
-
-SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg 
-SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end
-dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg
-
-echo ( '\nWater content in surface ' )
-echo ( 'SRF_mas_wat_beg  = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) )
-echo ( 'dMass (water srf) = {:12.3e} kg '.format (dSRF_mas_wat) )
-echo ( 'dMass (water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-6/ATM_rho) )
-echo ( 'dMass (water srf) = {:12.3e} m  '.format (dSRF_mas_wat/ATM_aire_sea_tot/ATM_rho) )
-
-echo ( '\nWater content in  ATM + SRF + RUN ' )
-echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '.
-           format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg,
-                   DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) )
-echo ( 'dMass (water atm+srf+run) = {:12.6e} kg '.format ( dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat) )
-echo ( 'dMass (water atm+srf+run) = {:12.3e} Sv '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/dtime_sec*1E-6/ATM_rho) )
-echo ( 'dMass (water atm+srf+run) = {:12.3e} m  '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/ATM_aire_sea_tot/ATM_rho) )
-
-echo ( '\n------------------------------------------------------------------------------------' )
-echo ( '-- Change in all components' )
-echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg'.
-           format (SEA_mas_wat_beg + DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg,
-                   SEA_mas_wat_end + DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) )
-echo ( 'dMass (water CPL)         = {:12.3e} kg '.format ( dSEA_mas_wat + dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat) )
-echo ( 'dMass (water CPL)         = {:12.3e} Sv '.format ((dSEA_mas_wat + dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/dtime_sec*1E-9) )
-echo ( 'dMass (water CPL)         = {:12.3e} m  '.format ((dSEA_mas_wat + dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/OCE_aire_tot*1E-3) )
-
-
+    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
+    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 = 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    = 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.
+    OCE_calv_icb = 0. ; OCE_mas_calv_icb = 0.
+
+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
+
+prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True)
+prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True)
+prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True)
+prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True)
+prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True)
+prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True)
+prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e',     )
+prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True)
+prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True)
+prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e',     )  
+prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True)
+prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True)
+prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True)
+prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True)
+prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True)
+prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True)
+prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True)
+prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True)
+prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e'      )
+prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True)
+prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True)
+prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True)
+prtFlux ('  RAIN         ', OCE_mas_rain      , 'e'      )
+prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True)
+prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True)
+prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True)
+prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True)
+prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
+
+
+echo (' ')
+
+prtFlux ( 'wbilo sea  ', ATM_flux_int (ATM_wbilo_sea), 'e',   )
+prtFlux ( 'costalflow ', ONE_flux_int (RUN_coastalflow), 'e',   )
+prtFlux ( 'riverflow  ', RUN_flx_river  , 'e',   )
+prtFlux ( 'costalflow ', RUN_flx_coastal, 'e',   )
+prtFlux ( 'runoff     ', RUN_flx_river+RUN_flx_coastal, 'e',   )
+
+ATM_to_OCE   =  ATM_flux_int (ATM_wbilo_sea) - RUN_flx_river - RUN_flx_coastal - ATM_flx_calving
+#OCE_from_ATM = -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf
+OCE_from_ATM = OCE_mas_all
+
+prtFlux ( 'ATM_to_OCE  ', ATM_to_OCE  , 'e', True )
+prtFlux ( 'OCE_from_ATM', OCE_from_ATM, 'e', True )