Index: /TOOLS/WATER_BUDGET/ATM_waterbudget.py
===================================================================
--- /TOOLS/WATER_BUDGET/ATM_waterbudget.py	(revision 6517)
+++ /TOOLS/WATER_BUDGET/ATM_waterbudget.py	(revision 6518)
@@ -12,10 +12,12 @@
 ##
 ##
-## SVN information
-#  $Author$
-#  $Date$
-#  $Revision$
-#  $Id$
-#  $HeadURL$
+# SVN Information
+SVN = {
+    'Author'  : "$Author$",
+    'Date'    : "$Date$",
+    'Revision': "$Revision$",
+    'Id'      : "$Id: ATM_waterbudget.py 6508 2023-06-13 10:58:38Z omamce $",
+    'HeadURL' : "$HeadUrl: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/WATER_BUDGET/ATM_waterbudget.py $"
+    }
 
 ###
@@ -34,19 +36,28 @@
 
 ## 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' ; RUN_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
+TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=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_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None ; ContinueOnError=False ; ErrorCount=0
-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
-SortIco = False
+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_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
 
 # Read command line arguments
@@ -74,5 +85,5 @@
             exec  ( f'{VarName} = wu.setNone ({VarName})' )
             exec  ( f'wu.{VarName} = {VarName}' )
-            print ( f'{VarName:25} set to : {locals()[VarName]:}' )
+            print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
             #exec ( f'del {VarName}' )
 
@@ -81,7 +92,7 @@
 ##-- 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 ( '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 or SI3
+if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3 or SI3
 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)
@@ -89,12 +100,13 @@
 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 
+if wu.unDefined ( 'YearLength' )   : YearLength  = 365.25 * 86400.    #-- 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}
-
-config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco}
+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 }
 
 ## --------------------------
@@ -106,5 +118,5 @@
 
 if libIGCM :
-    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
+    config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
 
 ##-- Import specific module
@@ -161,5 +173,6 @@
 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 }
+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 files
@@ -213,15 +226,23 @@
      
 echo ('\nOpen history files' )
-if wu.unDefined ( 'file_ATM_his' ) : 
-    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 wu.unDefined ( 'file_SRF_his' ) :
-    file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
+    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
 
@@ -449,7 +470,7 @@
 
 ## Write the full configuration
-config_out = open (FullIniFile, 'w')
-config.write ( config_out )
-config_out.close ()
+#config_out = open (FullIniFile, 'w')
+#config.write ( config_out )
+#config_out.close ()
 
 if ICO :
@@ -483,47 +504,79 @@
 # ATM grid with cell surfaces
 if LMDZ :
-    #ATM_lat  = lmdz.geo2point ( d_ATM_his ['lat']         , dim1D='cell' )
-    #ATM_lon  = lmdz.geo2point ( d_ATM_his ['lon']         , dim1D='cell' )
-    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire']     [0], cumulPoles=True, dim1D='cell' )
-    ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell' )
-    ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell' )
-    ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell' )
-    ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell' )
-    #SRF_lat  = lmdz.geo2point ( d_SRF_his ['lat']      [0], dim1D='cell' )
-    #SRF_lon  = lmdz.geo2point ( d_SRF_his ['lon']      [0], dim1D='cell' )
-    
-    #??SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell' )
-    SRF_aire =  ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'], dim1D='cell' )
+    ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' )
+    ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' )
+    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire']     [0], cumulPoles=True, dim1D='cell_latlon' )
+    ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' )
+    ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' )
+    ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' )
+    ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' )
+    SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' )
+    SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' )
+    SRF_aire      = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell_latlonl', cumulPoles=True )
+    SRF_areas     = lmdz.geo2point ( d_SRF_his ['Areas'],  dim1D='cell_latlonl', cumulPoles=True )
+    SRF_contfrac  = lmdz.geo2point ( d_SRF_his ['Contfrac'], dim1D='cell_latlonl' )
+
+    #SRF_res_aire = SRF_aire
     
 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
+        jpja, jpia = d_ATM_his['aire'][0].shape
+        if wu.unDefined ( 'file_ATM_aire' ) : 
+            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_lat  = d_ATM_his ['lat']
-        #ATM_lon  = d_ATM_his ['lon']
-        ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True, dim1D='cell' )
-        ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell' )
-        ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell' )
-        ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell' )
-        ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell' )
-        #SRF_lat  = d_SRF_his ['lat']
-        #SRF_lon  = d_SRF_his ['lon']
-        SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell' )
-        
+        try :
+            ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' )
+            ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' )
+        except :
+            ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat_dom_out']+0*d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' )
+            ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat_dom_out']+  d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' )
+        ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True, dim1D='cell_latlon' )
+        ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' )
+        ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' )
+        ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' )
+        ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' )
+
+    if SRF_HIS == 'latlon' :
+        try :
+            SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' )
+            SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' )
+        except :
+            try : 
+                SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_dom_out']+0*d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' )
+                SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_dom_out']+  d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' )
+            except : 
+                SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_domain_landpoints_out']+0*d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' )
+                SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_domain_landpoints_out']+  d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' )
+        SRF_areas    = d_SRF_his ['Areas'].values
+        SRF_areafrac = d_SRF_his ['AreaFrac'].values
+        SRF_areas    = xr.DataArray ( SRF_areas   , coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
+        SRF_areafrac = xr.DataArray ( SRF_areafrac, coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
+        SRF_contfrac = SRF_areafrac / SRF_areas
+        #SRF_aire      = lmdz.geo2point (SRF_areas * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
+        #SRF_aire      = lmdz.geo2point (SRF_areas, dim1D='cell_latlon', cumulPoles=True )
+        #SRF_aire      = lmdz.geo2point (d_ATM_aire ['aire'].squeeze().values * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
+        SRF_areas     = lmdz.geo2point ( SRF_areas   , dim1D='cell_latlon', cumulPoles=True )
+        SRF_areafrac  = lmdz.geo2point ( SRF_areafrac, dim1D='cell_latlon', cumulPoles=True )
+        SRF_contfrac  = lmdz.geo2point ( SRF_contfrac, dim1D='cell_latlon' )
+        SRF_aire      = SRF_areafrac
+
     if ATM_HIS == 'ico' :
         echo ( f'Grille icosaedre' )
         ATM_aire =  d_ATM_his ['aire']     [0][ATM_his_keysort].squeeze()
-        #ATM_lat  =  d_ATM_his ['lat']         [ATM_his_keysort]
-        #ATM_lon  =  d_ATM_his ['lon']         [ATM_his_keysort]
+        ATM_lat  =  d_ATM_his ['lat']         [ATM_his_keysort]
+        ATM_lon  =  d_ATM_his ['lon']         [ATM_his_keysort]
         ATM_fter =  d_ATM_his ['fract_ter'][0][ATM_his_keysort]
         ATM_foce =  d_ATM_his ['fract_oce'][0][ATM_his_keysort]
         ATM_fsic =  d_ATM_his ['fract_sic'][0][ATM_his_keysort]
         ATM_flic =  d_ATM_his ['fract_lic'][0][ATM_his_keysort]
-        #SRF_lat  =  d_SRF_his ['lat']         [SRF_his_keysort]
-        #SRF_lon  =  d_SRF_his ['lon']         [SRF_his_keysort]
-        SRF_aire =  d_SRF_his ['Areas'][SRF_his_keysort] * d_SRF_his ['Contfrac'][SRF_his_keysort]
+        
+    if SRF_HIS == 'ico' : 
+        SRF_lat       =  d_SRF_his ['lat'][SRF_his_keysort]
+        SRF_lon       =  d_SRF_his ['lon'][SRF_his_keysort]
+        SRF_aire      =  d_SRF_his ['Areas'][SRF_his_keysort] * d_SRF_his ['Contfrac'][SRF_his_keysort]
+        SRF_areas     =  d_SRF_his ['Areas'][SRF_his_keysort] 
+        SRF_contfrac  =  d_SRF_his ['Contfrac'][SRF_his_keysort]
 
 ATM_fsea      = ATM_foce + ATM_fsic
@@ -537,8 +590,13 @@
 
 SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0.)
+
+## Write the full configuration
+config_out = open (FullIniFile, 'w')
+config.write ( config_out )
+config_out.close ()
             
 if ICO :
     # Area on icosahedron grid
-    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze()
+    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
 
     if SortIco : 
@@ -552,5 +610,5 @@
 
     DYN_aire = d_DYN_aire['aire'][DYN_aire_keysort]
-    DYN_fsea = d_DYN_aire ['fract_oce'][DYN_aire_keysort] + d_DYN_aire['fract_sic'][DYN_aire_keysort]
+    DYN_fsea = d_DYN_aire['fract_oce'][DYN_aire_keysort] + d_DYN_aire['fract_sic'][DYN_aire_keysort]
     DYN_flnd = 1.0 - DYN_fsea
     DYN_fter = d_ATM_beg['FTER'][ATM_beg_keysort]
@@ -568,8 +626,8 @@
 
 # Functions computing integrals and sum
-def ATM_stock_int (stock) :
+def DYN_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
+    DYN_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
+    return DYN_stock_int
 
 def ATM_flux_int (flux) :
@@ -580,4 +638,5 @@
 def SRF_stock_int (stock) :
     '''Integrate (* surface) stock on land grid'''
+    #SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) # Marche avec ICO et lon/lat na
     SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
     return SRF_stock_int
@@ -598,14 +657,22 @@
     return ONE_flux_int
 
-
-#if LMDZ :
-#    d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} )
-
 ATM_aire_sea     = ATM_aire * ATM_fsea
 ATM_aire_tot     = ONE_stock_int (ATM_aire)
 ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea)
 
-
-echo ( 'Area of atmosphere/(4pi R^2) : {:12.5f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )
+DYN_aire_tot     =  ONE_stock_int ( DYN_aire )
+SRF_aire_tot     =  ONE_stock_int ( SRF_aire )
+
+echo ('')
+echo ( f'ATM DYN : Area of atmosphere : {DYN_aire_tot:18.9e}' )
+echo ( f'ATM HIS : Area of atmosphere : {ATM_aire_tot:18.9e}' )
+echo ( f'ATM SRF : Area of atmosphere : {SRF_aire_tot:18.9e}' )
+echo ('')
+echo ( 'ATM DYN : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) )
+echo ( 'ATM HIS : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )
+echo ( 'ATM SRF : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) )
+echo ('')
+echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' )
+
 
 if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) :
@@ -627,6 +694,6 @@
     DYN_psol_end = d_DYN_end['ps'][DYN_end_keysort]
 if LMDZ : 
-    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' )
-    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' )
+    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' )
+    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' )
     
 echo ( '3D Pressure at the interface layers (not scalar points)' )
@@ -635,12 +702,10 @@
 
 klevp1 = ATM_Bhyb.shape[-1]
-cell        = DYN_psol_beg.shape[-1]
-klev = klevp1 - 1
+cell   = DYN_psol_beg.shape[-1]
+klev   = klevp1 - 1
 
 echo ( 'Layer thickness (pressure)' )
 DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
 DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
-#DYN_dpres_beg = DYN_dpres_beg
-#DYN_dpres_end = DYN_dpres_end
     
 for k in np.arange (klevp1-1) :
@@ -669,6 +734,6 @@
     
 echo ( 'Compute water content : vertical and horizontal integral' )
-DYN_mas_wat_beg = ATM_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav
-DYN_mas_wat_end = ATM_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav
+DYN_mas_wat_beg = DYN_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav
+DYN_mas_wat_end = DYN_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav
 
 echo ( 'Variation of water content' )
@@ -683,9 +748,13 @@
 prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True )
 
-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_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']
@@ -707,23 +776,4 @@
 ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC']  
 
-# if LMDZ :
-#     ATM_sno_beg       = lmdz.geo2point ( ATM_sno_beg    , dim1D='cell' )
-#     ATM_sno_end       = lmdz.geo2point ( ATM_sno_end    , dim1D='cell' )
-#     ATM_qs_beg        = lmdz.geo2point ( ATM_qs_beg     , dim1D='cell' )
-#     ATM_qs_end        = lmdz.geo2point ( ATM_qs_end     , dim1D='cell' )
-#     ATM_qsol_beg      = lmdz.geo2point ( ATM_qsol_beg   , dim1D='cell' )
-#     ATM_qsol_end      = lmdz.geo2point ( ATM_qsol_end   , dim1D='cell' )
-#     ATM_qs01_beg      = lmdz.geo2point ( ATM_qs01_beg   , dim1D='cell' )
-#     ATM_qs01_end      = lmdz.geo2point ( ATM_qs01_end   , dim1D='cell' )
-#     ATM_qs02_beg      = lmdz.geo_point ( ATM_qs02_beg   , dim1D='cell' )
-#     ATM_qs02_end      = lmdz.geo2point ( ATM_qs02_end   , dim1D='cell' )
-#     ATM_qs03_beg      = lmdz.geo2point ( ATM_qs03_beg   , dim1D='cell' )
-#     ATM_qs03_end      = lmdz.geo2point ( ATM_qs03_end   , dim1D='cell' )
-#     ATM_qs04_beg      = lmdz.geo2point ( ATM_qs04_beg   , dim1D='cell' )
-#     ATM_qs04_end      = lmdz.geo2point ( ATM_qs04_end   , dim1D='cell' )
-#     LIC_sno_beg       = lmdz.geo2point ( LIC_sno_beg    , dim1D='cell' )
-#     LIC_sno_end       = lmdz.geo2point ( LIC_sno_end    , dim1D='cell' )
-#     LIC_runlic0_beg   = lmdz.geo2point ( LIC_runlic0_beg, dim1D='cell' )
-#     LIC_runlic0_end   = lmdz.geo2point ( LIC_runlic0_end, dim1D='cell' )
 if ICO :
      ATM_sno_beg       = ATM_sno_beg    [ATM_beg_keysort]
@@ -745,30 +795,28 @@
      LIC_runlic0_beg   = LIC_runlic0_beg[ATM_beg_keysort]
      LIC_runlic0_end   = LIC_runlic0_end[ATM_end_keysort]
-
-
    
 LIC_qs_beg = ATM_qs02_beg
 LIC_qs_end = ATM_qs02_end
 
-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 )
-
-LIC_mas_sno_beg     = ATM_stock_int ( LIC_sno_beg     )
-LIC_mas_sno_end     = ATM_stock_int ( LIC_sno_end     )
-LIC_mas_runlic0_beg = ATM_stock_int ( LIC_runlic0_beg )
-LIC_mas_runlic0_end = ATM_stock_int ( LIC_runlic0_end )
+ATM_mas_sno_beg     = DYN_stock_int ( ATM_sno_beg  )
+ATM_mas_sno_end     = DYN_stock_int ( ATM_sno_end  )
+
+ATM_mas_qs_beg      = DYN_stock_int ( ATM_qs_beg   )
+ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   )
+ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg )
+ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end )
+ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg )
+ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end )
+ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg )
+ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end )
+ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg )
+ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end )
+ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg )
+ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end )
+
+LIC_mas_sno_beg     = DYN_stock_int ( LIC_sno_beg     )
+LIC_mas_sno_end     = DYN_stock_int ( LIC_sno_end     )
+LIC_mas_runlic0_beg = DYN_stock_int ( LIC_runlic0_beg )
+LIC_mas_runlic0_end = DYN_stock_int ( LIC_runlic0_end )
 
 LIC_mas_qs_beg  = ATM_mas_qs02_beg
@@ -826,20 +874,22 @@
 
 if Routing == 'SECHIBA' :
-        RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
-        RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
-        RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
-        RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
-        RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
-        RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
-        
-        RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
-        RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
-        RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
-        RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
-        RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
-        RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
-
-RUN_mas_wat_beg = RUN_mas_wat_fast_beg + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg
-RUN_mas_wat_end = RUN_mas_wat_fast_end + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end
+    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
+    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
+    RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
+    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
+    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
+    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
+    
+    RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
+    RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
+    RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
+    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
+    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
+    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
+
+RUN_mas_wat_beg = RUN_mas_wat_fast_beg  + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + \
+                  RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg
+RUN_mas_wat_end = RUN_mas_wat_fast_end  + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + \
+                  RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end
 
 dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg
@@ -874,5 +924,5 @@
 echo ( f'\nWater content in routing  -- {Title} ' )
 echo ( '------------------------------------------------------------------------------------' )
-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 ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' )
 prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   )
 
@@ -890,12 +940,12 @@
 
 if LMDZ :
-    SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell')
-    SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell')
-    SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell')
-    SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell')
-    SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell')
-    SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell')
-    SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell')
-    SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell')
+    SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell_latlon')
+    SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell_latlon')
+    SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell_latlon')
+    SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell_latlon')
+    SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell_latlon')
+    SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell_latlon')
+    SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell_latlon')
+    SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell_latlon')
   
 if ICO :
@@ -942,5 +992,4 @@
 dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg
 
-
 echo ( '------------------------------------------------------------------------------------' )
 echo ( f'\nSurface reservoirs  -- {Title} ')
@@ -950,8 +999,8 @@
 echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' )
 
-prtFlux ('dMass (watveg) ', dSRF_mas_watveg , 'e' , True)
-prtFlux ('dMass (watsoil)', dSRF_mas_watsoil, 'e' , True)
-prtFlux ('dMass (snow)   ', dSRF_mas_snow   , 'e' , True)
-prtFlux ('dMass (lake)   ', dSRF_mas_lake   , 'e' , True)
+prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True )
+prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True )
+prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True )
+prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True )
 
 SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg 
@@ -976,19 +1025,19 @@
 print ( 'Step 1', end='' )
 if ATM_HIS == 'latlon' :
-    ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell' )
-    ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic'], dim1D='cell' )
-    ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter'], dim1D='cell' )
-    ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic'], dim1D='cell' )
-    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell' )
-    ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving'], dim1D='cell' )
-    ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']  , dim1D='cell' )
-    ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']   , dim1D='cell' )
-    ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']     , dim1D='cell' )
-    ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']     , dim1D='cell' )
-    ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell' )
-    ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell' )
-    ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell' )
-    ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell' )
-    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell' )
+    ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell_latlon' )
+    ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic'], dim1D='cell_latlon' )
+    ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter'], dim1D='cell_latlon' )
+    ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic'], dim1D='cell_latlon' )
+    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' )
+    ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving'], dim1D='cell_latlon' )
+    ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']  , dim1D='cell_latlon' )
+    ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']   , dim1D='cell_latlon' )
+    ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']     , dim1D='cell_latlon' )
+    ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']     , dim1D='cell_latlon' )
+    ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell_latlon' )
+    ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell_latlon' )
+    ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' )
+    ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' )
+    ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' )
     
 print ( ' - Step 2', end='' )
@@ -1024,12 +1073,12 @@
 print ( ' - Step 4', end='' )
 if RUN_HIS == 'latlon' : 
-    RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell' ) 
-    RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']  , dim1D='cell' )
-    RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']     , dim1D='cell' )
-    RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']   , dim1D='cell' )
-    RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']  , dim1D='cell' )
-    
-    RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'], dim1D='cell' ) 
-    RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']  , dim1D='cell' )
+    RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell_latlon' ) 
+    RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']  , dim1D='cell_latlon' )
+    RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']     , dim1D='cell_latlon' )
+    RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']   , dim1D='cell_latlon' )
+    RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']  , dim1D='cell_latlon' )
+    
+    RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'], dim1D='cell_latlon' ) 
+    RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']  , dim1D='cell_latlon' )
 
 print ( ' - Step 5', end='' )
@@ -1046,9 +1095,9 @@
 print ( ' - Step 6', end='' )
 if SRF_HIS == 'latlon' :
-    SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  , dim1D='cell')
-    SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  , dim1D='cell')
-    SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] , dim1D='cell')
-    SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] , dim1D='cell')
-    SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell' )
+    SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  , dim1D='cell_latlon')
+    SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  , dim1D='cell_latlon')
+    SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] , dim1D='cell_latlon')
+    SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] , dim1D='cell_latlon')
+    SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell_latlon' )
 
 print ( '- Step 7', end='' )
@@ -1140,5 +1189,9 @@
 prtFlux ('precip        ', ATM_flx_precip     , 'f' )       
 prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )        
-prtFlux ('evap          ', ATM_flx_evap       , 'f' )         
+prtFlux ('evap          ', ATM_flx_evap       , 'f' )
+prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )   
+
+echo ( '\n====================================================================================' )
+echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
 prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
 prtFlux ('riverflow     ', RUN_flx_river      , 'f' )        
@@ -1152,5 +1205,4 @@
 prtFlux ('river out     ', RUN_flx_output     , 'f' )   
 prtFlux ('river bil     ', RUN_flx_bil        , 'f' )   
-prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )   
 
 ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river
@@ -1277,7 +1329,7 @@
 echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) )
 
-# PRINT *,"routing water Balance ;  before : ",Â water_balance_before," ; after : ",water_balance_after,  &a
-# 1150	            " ; delta : ",Â 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%"
-
 echo ( ' ' )
 echo ( f'{Title = }' )
+
+echo ( 'SVN Information' )
+for clef in SVN.keys () : echo ( SVN[clef] )
Index: /TOOLS/WATER_BUDGET/WaterUtils.py
===================================================================
--- /TOOLS/WATER_BUDGET/WaterUtils.py	(revision 6517)
+++ /TOOLS/WATER_BUDGET/WaterUtils.py	(revision 6518)
@@ -15,9 +15,9 @@
 ##
 ## SVN information
-#  $Author: omamce $
-#  $Date: 2022-12-08 10:24:05 +0100 (Thu, 08 Dec 2022) $
-#  $Revision: 6277 $
+#  $Author$
+#  $Date$
+#  $Revision$
 #  $Id: ATM_waterbudget.py 6277 2022-12-08 09:24:05Z omamce $
-#  $HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/WATER_BUDGET/ATM_waterbudget.py $
+#  $HeadURL$
 
 ## Import system modules
