diff --git a/MY_RUN/KTEST/002_3Drelief/003_python/plot_002_3DRelief.py b/MY_RUN/KTEST/002_3Drelief/003_python/plot_002_3DRelief.py
index 079e324052720e530282675cb68dbaad869336eb..6f57fae700bef1f17d9783d15172b05175cf3e73 100644
--- a/MY_RUN/KTEST/002_3Drelief/003_python/plot_002_3DRelief.py
+++ b/MY_RUN/KTEST/002_3Drelief/003_python/plot_002_3DRelief.py
@@ -10,7 +10,6 @@ mpl.use('Agg')
 from read_MNHfile import read_netcdf
 from Panel_Plot import PanelPlot
 from misc_functions import comp_altitude2DVar
-import cartopy.crs as ccrs
 import os
 
 os.system('rm -f tempgraph*')
@@ -21,8 +20,7 @@ path=""
 
 LnameFiles = ['REL3D.1.EXP01.002.nc']
 
-Dvar_input = {
-'f1':['ZS', 'UT', 'WT','ni_u','nj_u','level','ZTOP', 'ni','nj','level_w','time']}
+Dvar_input = {'f1':['ZS', 'UT', 'WT','ni_u','nj_u','level','ZTOP', 'ni','nj','level_w','time']}
 
 #  Read the variables in the files
 Dvar = {}
@@ -53,7 +51,6 @@ Lstep = [0.25, 0.1]
 Lstepticks = Lstep
 Lcolormap = ['gist_rainbow_r']*len(Lplot)
 Ltime = [Dvar['f1']['time']]*len(Lplot)
-Lprojection = [ccrs.PlateCarree()]*len(Lplot)
 fig1 = Panel1.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel,
                                 orog=orog, Lylim=Lylim, colorbar=True, Ltime=Ltime)
diff --git a/MY_RUN/KTEST/004_Reunion/007_python/plot_004_Reunion.py b/MY_RUN/KTEST/004_Reunion/007_python/plot_004_Reunion.py
index 2a607c828d7eff707b7bb5d75f8a68be83f5ff0a..ea59e343330e5dfc43a383bb976da141751a6705 100644
--- a/MY_RUN/KTEST/004_Reunion/007_python/plot_004_Reunion.py
+++ b/MY_RUN/KTEST/004_Reunion/007_python/plot_004_Reunion.py
@@ -6,15 +6,13 @@ Creation : 07/01/2021
 
 Last modifications
 """
-
 import matplotlib as mpl
 mpl.use('Agg')
 from read_MNHfile import read_netcdf
 from Panel_Plot import PanelPlot
-from misc_functions import comp_altitude2DVar, windvec_verti_proj, mean_operator
+from misc_functions import comp_altitude2DVar, mean_operator
 import cartopy.crs as ccrs
 import numpy as np
-import math
 import copy
 import os
 
@@ -22,8 +20,6 @@ os.system('rm -f tempgraph*')
 #
 #  User's parameter / Namelist
 #
-path=""
-
 LnameFiles = ['REUNI.1.00A20.004dia.nc', 'REUNI.1.00A20.004.nc']
 
 Dvar_input = {
@@ -32,7 +28,7 @@ Dvar_input = {
 
 #  Read the variables in the files
 Dvar = {}
-Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
+Dvar = read_netcdf(LnameFiles, Dvar_input, path="", removeHALO=True)
 
 ################################################################
 #########          PANEL 1  # Horizontal cross-section
@@ -40,7 +36,7 @@ Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
 Panel1 = PanelPlot(2,2, [20,20],'004_Reunion horizontal sections')
 
 Dvar['f1']['WIND'] = np.sqrt(Dvar['f1']['UT']**2 + Dvar['f1']['VT']**2)
-Lplot = [ Dvar['f1']['ZS'][:,:], Dvar['f1']['WIND'][0,:,:], Dvar['f1']['ALT_THETA'][0,:,:], Dvar['f1']['ALT_PRESSURE'][0,:,:]]
+Lplot = [ Dvar['f1']['ZS'][:,:], Dvar['f1']['WIND'][0,:,:], Dvar['f1']['ALT_THETA'][:,:], Dvar['f1']['ALT_PRESSURE'][:,:]]
 
 LaxeX = [Dvar['f1']['longitude']]*len(Lplot)
 LaxeY = [Dvar['f1']['latitude']]*len(Lplot)
@@ -79,8 +75,6 @@ Lscale = [400]*len(Lplot1)
 fig2 = Panel1.pvector(Lxx=LaxeX, Lyy=LaxeY, Llevel=Llvl, Lvar1=Lplot1, Lvar2=Lplot2, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
                       Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lproj=Lprojection, Lid_overlap=[2,6], ax=fig1.axes, Lscale=Lscale)
 
-
-
 ################################################################
 #########          PANEL 2  # Vertical cross-section
 ###############################################################
@@ -141,4 +135,4 @@ Lcolor=['lightgray']
 fig4 = Panel2.pvector(Lxx=LaxeX, Lyy=LaxeZ, Lvar1=Lplot1, Lvar2=Lplot2, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
                         Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[6], ax=fig3.axes, Lscale=Lscale, Lylim=Lylim, Lxlim=Lxlim, Lcolor=Lcolor)
 
-Panel2.save_graph(2,fig4)
+Panel2.save_graph(2,fig4)
\ No newline at end of file
diff --git a/MY_RUN/KTEST/005_ARM/003_python/plot_005_ARM.py b/MY_RUN/KTEST/005_ARM/003_python/plot_005_ARM.py
index 7427083f029f7a175f234f47b3bed95fa23a28a2..46764697097f879f0c99455ba9a268a071b76511 100644
--- a/MY_RUN/KTEST/005_ARM/003_python/plot_005_ARM.py
+++ b/MY_RUN/KTEST/005_ARM/003_python/plot_005_ARM.py
@@ -17,28 +17,29 @@ os.system('rm -f tempgraph*')
 #  User's parameter / Namelist
 #
 path=""
-
 LnameFiles = ['ARM__.1.CEN4T.000.nc' ]
+LG_LES = '/LES_budgets/Cartesian/Not_time_averaged/Not_normalized/cart/'
 
 Dvar_input = {
-'f1':['MEAN_TH','MEAN_U','MEAN_V','MEAN_RC','MEAN_RR',
-      'SBG_TKE','SBG_WTHL','SBG_WRT',
-      'THLUP_MF','RTUP_MF','RVUP_MF','RCUP_MF','RIUP_MF','WUP_MF',
-      'MAFLX_MF','DETR_MF','ENTR_MF','FRCUP_MF','THVUP_MF','WTHL_MF',
-      'WRT_MF','WTHV_MF','WU_MF','WV_MF',
+'f1':[(LG_LES,'MEAN_TH') , (LG_LES,'MEAN_U')  , (LG_LES,'MEAN_V') , (LG_LES,'MEAN_RC'), (LG_LES,'MEAN_RR'),
+      (LG_LES,'SBG_TKE') , (LG_LES,'SBG_WTHL'), (LG_LES,'SBG_WRT'),
+      (LG_LES,'THLUP_MF'), (LG_LES,'RTUP_MF') , (LG_LES,'RVUP_MF'), (LG_LES,'RCUP_MF'), (LG_LES,'RIUP_MF'), (LG_LES,'WUP_MF'),
+      (LG_LES,'MAFLX_MF'), (LG_LES,'DETR_MF') , (LG_LES,'ENTR_MF'), (LG_LES,'FRCUP_MF'), (LG_LES,'THVUP_MF'), (LG_LES,'WTHL_MF'),
+      (LG_LES,'WRT_MF')  , (LG_LES,'WTHV_MF') , (LG_LES,'WU_MF')  , (LG_LES,'WV_MF'),
       'level_les','time_les']
 }
 
 #  Read the variables in the files
 Dvar = {}
-Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=False)
+Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=False, get_data_only=True)
 
 ################################################################
 #########          PANEL 1
 ###############################################################
 Panel1 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.03, labelcolorbarpad = 13, colorbaraspect=40)
 
-Lplot = [Dvar['f1']['MEAN_TH'],Dvar['f1']['MEAN_U'], Dvar['f1']['MEAN_V'], Dvar['f1']['MEAN_RC'], Dvar['f1']['MEAN_RR'],Dvar['f1']['SBG_TKE'],]
+Lplot = [Dvar['f1'][(LG_LES,'MEAN_TH')],Dvar['f1'][(LG_LES,'MEAN_U')], Dvar['f1'][(LG_LES,'MEAN_V')], 
+         Dvar['f1'][(LG_LES,'MEAN_RC')], Dvar['f1'][(LG_LES,'MEAN_RR')],Dvar['f1'][(LG_LES,'SBG_TKE')]]
 LaxeX = [Dvar['f1']['time_les']/3600.]*len(Lplot)
 LaxeZ = [Dvar['f1']['level_les']]*len(Lplot)
 Ltitle = ['Mean potential temperature TH', 'Mean U', 'Mean V', 'Mean cloud mixing ratio RC', 'Mean precipitation RR', 'Subgrid TKE']
@@ -57,7 +58,6 @@ LaddWhite = [False, False, True, True,True, True]
 fig1 = Panel1.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 LaddWhite_cm=LaddWhite, Lylim=Lylim)
-
 Panel1.save_graph(1,fig1)
 
 ################################################################
@@ -65,7 +65,8 @@ Panel1.save_graph(1,fig1)
 ###############################################################
 Panel2 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.03, labelcolorbarpad = 13, colorbaraspect=40)
 
-Lplot = [Dvar['f1']['SBG_WTHL'], Dvar['f1']['SBG_WRT'], Dvar['f1']['THLUP_MF'], Dvar['f1']['RTUP_MF'], Dvar['f1']['RVUP_MF'], Dvar['f1']['RCUP_MF']]
+Lplot = [Dvar['f1'][(LG_LES,'SBG_WTHL')], Dvar['f1'][(LG_LES,'SBG_WRT')], Dvar['f1'][(LG_LES,'THLUP_MF')], 
+         Dvar['f1'][(LG_LES,'RTUP_MF')], Dvar['f1'][(LG_LES,'RVUP_MF')], Dvar['f1'][(LG_LES,'RCUP_MF')]]
 LaxeX = [Dvar['f1']['time_les']/3600.]*len(Lplot)
 LaxeZ = [Dvar['f1']['level_les']]*len(Lplot)
 Ltitle = ['Subgrid vertical liquid potential temp. flux', 'Subgrid vertical RT flux', 
@@ -85,7 +86,6 @@ LaddWhite = [False, False, False, False, False, False]
 fig2 = Panel2.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 LaddWhite_cm=LaddWhite, Lylim=Lylim)
-
 Panel2.save_graph(2,fig2)
 
 ################################################################
@@ -93,7 +93,8 @@ Panel2.save_graph(2,fig2)
 ###############################################################
 Panel3 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.03, labelcolorbarpad = 13, colorbaraspect=40)
 
-Lplot = [Dvar['f1']['RIUP_MF'], Dvar['f1']['WUP_MF'], Dvar['f1']['MAFLX_MF'], Dvar['f1']['DETR_MF'], Dvar['f1']['ENTR_MF'], Dvar['f1']['FRCUP_MF']]
+Lplot = [Dvar['f1'][(LG_LES,'RIUP_MF')], Dvar['f1'][(LG_LES,'WUP_MF')], Dvar['f1'][(LG_LES,'MAFLX_MF')], 
+         Dvar['f1'][(LG_LES,'DETR_MF')], Dvar['f1'][(LG_LES,'ENTR_MF')], Dvar['f1'][(LG_LES,'FRCUP_MF')]]
 LaxeX = [Dvar['f1']['time_les']/3600.]*len(Lplot)
 LaxeZ = [Dvar['f1']['level_les']]*len(Lplot)
 Ltitle = ['Updraft ice mixing ratio', 'Updraft vertical velocity', 
@@ -113,16 +114,15 @@ LaddWhite = [True]*len(Lplot)
 fig3 = Panel3.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 LaddWhite_cm=LaddWhite, Lylim=Lylim)
-
 Panel3.save_graph(3,fig3)
 
-
 ################################################################
 #########          PANEL 4
 ###############################################################
 Panel4 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.03, labelcolorbarpad = 13, colorbaraspect=40)
 
-Lplot = [Dvar['f1']['THVUP_MF'], Dvar['f1']['WTHL_MF'], Dvar['f1']['WRT_MF'], Dvar['f1']['WTHV_MF'], Dvar['f1']['WU_MF'], Dvar['f1']['WV_MF']]
+Lplot = [Dvar['f1'][(LG_LES,'THVUP_MF')], Dvar['f1'][(LG_LES,'WTHL_MF')], Dvar['f1'][(LG_LES,'WRT_MF')], 
+         Dvar['f1'][(LG_LES,'WTHV_MF')], Dvar['f1'][(LG_LES,'WU_MF')], Dvar['f1'][(LG_LES,'WV_MF')]]
 LaxeX = [Dvar['f1']['time_les']/3600.]*len(Lplot)
 LaxeZ = [Dvar['f1']['level_les']]*len(Lplot)
 Ltitle = ['Updraft virtual potential temperature', 'Subgrid WTHL flux from Mass-Flux scheme', 
@@ -142,5 +142,4 @@ LaddWhite = [False, False, True, False, False, False]
 fig4 = Panel4.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 LaddWhite_cm=LaddWhite, Lylim=Lylim)
-
-Panel4.save_graph(4,fig4)
+Panel4.save_graph(4,fig4)
\ No newline at end of file
diff --git a/MY_RUN/KTEST/007_16janvier/010_python/plot_007_16janvier.py b/MY_RUN/KTEST/007_16janvier/010_python/plot_007_16janvier.py
index 7e32b6b951b437dcba03357de097726a2bd5d985..9bc7193671db8a65ad0364bbd55d7a4b4ea65fb5 100644
--- a/MY_RUN/KTEST/007_16janvier/010_python/plot_007_16janvier.py
+++ b/MY_RUN/KTEST/007_16janvier/010_python/plot_007_16janvier.py
@@ -32,7 +32,7 @@ Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
 ################################################################
 #########          PANEL 1  
 ###############################################################
-Panel1 = PanelPlot(2,2, [20,20],'007_janvier domaine 1 16JAN.1.12B18.001dg.nc')
+Panel1 = PanelPlot(2,2, [20,20],'007_janvier domaine 1 16JAN.1.12B18.001dg.nc', minmaxpad=1.05)
 
 Lplot = [ Dvar['f1']['ZS'],Dvar['f1']['THT850HPA'], Dvar['f1']['MRV700HPA'],Dvar['f1']['ALT_PRESSURE']]
 lon = [Dvar['f1']['longitude']]*len(Lplot)
@@ -73,7 +73,7 @@ Panel1.save_graph(1,fig2)
 ################################################################
 #########          PANEL 2 
 ###############################################################
-Panel2 = PanelPlot(2,2, [20,20],'007_janvier domaine 2 16JAN.1.12B18.001dg.nc')
+Panel2 = PanelPlot(2,2, [20,20],'007_janvier domaine 2 16JAN.1.12B18.001dg.nc', minmaxpad=1.05)
 
 Lplot = [ Dvar['f2']['ZS'],Dvar['f2']['THT850HPA'], Dvar['f2']['MRV700HPA'],Dvar['f2']['ALT_PRESSURE']]
 lon = [Dvar['f2']['longitude']]*len(Lplot)
diff --git a/MY_RUN/KTEST/011_KW78CHEM/004_python/plot_011_KW78CHEM.py b/MY_RUN/KTEST/011_KW78CHEM/004_python/plot_011_KW78CHEM.py
index d271fe993d8733a11387ff34e6ebd728f36d6e3c..3f792a31363a78a8d5db2c9e5163afa1597afc82 100644
--- a/MY_RUN/KTEST/011_KW78CHEM/004_python/plot_011_KW78CHEM.py
+++ b/MY_RUN/KTEST/011_KW78CHEM/004_python/plot_011_KW78CHEM.py
@@ -1,6 +1,5 @@
 #!/usr/bin/env python3
 """
-
 @author: Quentin Rodier
 Creation : 07/01/2021
 
@@ -36,7 +35,7 @@ Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
 ################################################################
 #########          PANEL 1
 ###############################################################
-Panel1 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.01)
+Panel1 = PanelPlot(2,3, [25,14],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.01, lateralminmaxpad=0.9)
 
 Lplot = [ Dvar['f1']['INPRR'], Dvar['f1']['ACPRR'], Dvar['f1']['PABST'],Dvar['f2']['ALT_CLOUD'],Dvar['f2']['ALT_CLOUD'] ]
 
@@ -92,7 +91,7 @@ Panel1.save_graph(1,fig2)
 ###############################################################
 Panel2 = PanelPlot(2,2, [20,20],'', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.01, colorbaraspect=40, labelcolorbarpad = 13)
 
-# Interpoler COT','O3T à 3000 et 5000m avec une moyenne sur2 niveaux
+# Interpoler COT','O3T à 3000 et 5000m avec une moyenne sur 2 niveaux
 Dvar['f1']['COT3000m'] = (Dvar['f1']['COT'][6,:,:] + Dvar['f1']['COT'][5,:,:])/2.0
 Dvar['f1']['O3T3000m'] = (Dvar['f1']['O3T'][6,:,:] + Dvar['f1']['O3T'][5,:,:])/2.0
 Dvar['f1']['COT5000m'] = (Dvar['f1']['COT'][10,:,:] + Dvar['f1']['COT'][9,:,:])/2.0
@@ -131,7 +130,6 @@ Dvar['f1']['UM'] = tomass.MXM(Dvar['f1']['UT'])
 Dvar['f1']['VM'] = tomass.MYM(Dvar['f1']['VT'])
 Dvar['f1']['WM'] = tomass.MZM(Dvar['f1']['WT'])
 
-
 angle_sec1, RVT_sec1, axe_m1 = oblique_proj(Dvar['f1']['RVT'], Dvar['f1']['ni'], Dvar['f1']['nj'], Dvar['f1']['level'], i_beg, j_beg, i_end, j_end)
 WIND_proj = windvec_verti_proj(Dvar['f1']['UM'], Dvar['f1']['VM'], Dvar['f1']['level'], angle_sec1)
 angle_sec1, WIND_sec1, axe_m1 = oblique_proj(WIND_proj, Dvar['f1']['ni'], Dvar['f1']['nj'], Dvar['f1']['level'], i_beg, j_beg, i_end, j_end)
@@ -161,7 +159,6 @@ fig3 = Panel3.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lyl
                                 Lstep=Lstep, Lstepticks=Lstepticks, Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 Ltime=Ltime, Lpltype=Lpltype, LaddWhite_cm=LaddWhite)
 
-
 Lplot1 = [ WIND_sec1]
 Lplot2 = [ WT_sec1]
 Ltitle = ['Wind']
@@ -174,7 +171,6 @@ Lscale = [200]*len(Lplot)
 fig4 = Panel3.pvector(Lxx=LaxeX, Lyy=LaxeZ, Lvar1=Lplot1, Lvar2=Lplot2, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
                         Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[0], ax=fig3.axes, Lscale=Lscale)
 
-
 Lplot = [RRT_sec1]
 LaxeX = [axe_m1]
 LaxeZ = [Dvar['f1']['level']]
@@ -196,7 +192,6 @@ LaddWhite = [True]*len(Lplot)
 fig5 = Panel3.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lminval=Lminval, Lmaxval=Lmaxval, 
                                 Lstep=Lstep, Lstepticks=Lstepticks, LcolorLine=LcolorLine, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 Ltime=Ltime, Lpltype=Lpltype, LaddWhite_cm=LaddWhite, ax=fig4.axes,Lid_overlap=[2],colorbar=False)
-
 Panel3.save_graph(3,fig5)
 
 ################################################################
@@ -292,4 +287,4 @@ fig8 = Panel5.psectionV(Lxx=LaxeX, Lzz=LaxeZ, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lyl
                                 Lstep=Lstep, Lstepticks=Lstepticks, LcolorLine=LcolorLine, Lcbarlabel=Lcbarlabel, Lfacconv=Lfacconv, 
                                 Lpltype=Lpltype, ax=fig7.axes,Lid_overlap=[0, 2, 4, 6, 8, 10, 12, 14, 16],colorbar=False)
 
-Panel5.save_graph(5,fig8)
+Panel5.save_graph(5,fig8)
\ No newline at end of file
diff --git a/MY_RUN/KTEST/012_dust/005_python/plot_012_dust.py b/MY_RUN/KTEST/012_dust/005_python/plot_012_dust.py
index 7448f10d0b38061642aacbb934f40796745cc745..2116f12329b9163e64761721f006b0e9e06efd1f 100644
--- a/MY_RUN/KTEST/012_dust/005_python/plot_012_dust.py
+++ b/MY_RUN/KTEST/012_dust/005_python/plot_012_dust.py
@@ -9,7 +9,7 @@ import matplotlib as mpl
 mpl.use('Agg')
 from read_MNHfile import read_netcdf
 from Panel_Plot import PanelPlot
-from misc_functions import mean_operator, convert_date
+from misc_functions import mean_operator
 import cartopy.crs as ccrs
 import numpy as np
 import os
@@ -21,8 +21,7 @@ os.system('rm -f tempgraph*')
 path=""
 LnameFiles = ['DUST7.1.SEG02.004.nc']
 
-Dvar_input = {
-'f1':['ZS', 'UT','VT', 'WT','THT',
+Dvar_input = {'f1':['ZS', 'UT','VT', 'WT','THT',
       'DSTM03T','DSTM33T','DSTM02T','DSTM32T','DSTM01T','DSTM31T','F_DST001P1','F_DST002P1','F_DST003P1',
       'latitude','longitude','level',
       'INPRR','ACPRR','PABST','RCT','RVT','RRT','LSTHM']}
@@ -30,6 +29,7 @@ Dvar_input = {
 #  Read the variables in the files
 Dvar = {}
 Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
+
 ################################################################
 #########          PANEL 1
 ###############################################################
@@ -74,6 +74,7 @@ fig2 = Panel1.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Llevel=Llvl, Lvar2=Lplot2,
                         Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[0], ax=fig1.axes, Lscale=Lscale)
 
 Panel1.save_graph(1,fig2)
+
 ################################################################
 #########          PANEL 2
 ###############################################################
@@ -97,6 +98,7 @@ fig3 = Panel2.psectionH(lon=lon, lat=lat, Lvar=Lplot, Llevel=Llvl, Lxlab=Lxlab,
                                 Ltime=Ltime, LaddWhite_cm=LaddWhite, Lproj=Lprojection, ax=[])
 
 Panel2.save_graph(2,fig3)
+
 ################################################################
 #########          PANEL 3
 ###############################################################
@@ -119,4 +121,4 @@ fig4 = Panel3.psectionH(lon=lon, lat=lat, Lvar=Lplot, Lxlab=Lxlab, Lylab=Lylab,
                                 Lstep=[], Lstepticks=[], Lcolormap=Lcolormap, Lcbarlabel=Lcbarlabel,Lcarte=Lcarte,
                                 Ltime=Ltime, LaddWhite_cm=LaddWhite, Lproj=Lprojection, ax=[])
 
-Panel3.save_graph(3,fig4)
+Panel3.save_graph(3,fig4)
\ No newline at end of file
diff --git a/MY_RUN/KTEST/014_LIMA/004_python/plot_014_LIMA.py b/MY_RUN/KTEST/014_LIMA/004_python/plot_014_LIMA.py
index b64203496df724da0109d0798a84d4b66302de4e..61f335890fa2c86098b3b9e1aeb1c33eaa038cb2 100644
--- a/MY_RUN/KTEST/014_LIMA/004_python/plot_014_LIMA.py
+++ b/MY_RUN/KTEST/014_LIMA/004_python/plot_014_LIMA.py
@@ -20,25 +20,27 @@ os.system('rm -f tempgraph*')
 path=""
 LnameFiles = ['XPREF.1.SEG01.000.nc' ]
 
-Dvar_input = {'f1':[('RI','AVEF'), ('CICE','AVEF'), ('RS','AVEF'),('RG','AVEF'), ('CIFNFREE01','AVEF'),('CIFNNUCL01','AVEF') ]}
+Dvar_input = {'f1':[('/Budgets/RI','AVEF'), ('/Budgets/CICE','AVEF'), ('/Budgets/RS','AVEF'),
+                    ('/Budgets/RG','AVEF'), ('/Budgets/CIFNFREE01','AVEF'),('/Budgets/CIFNNUCL01','AVEF') ]}
 Dvar_input_coord_budget = {'f1':['cart_level', 'cart_ni']}
 Dvar_input_coord = {'f1':['ZS','ZTOP']}
 
 #  Read the variables in the files
 Dvar, Dvar_coord_budget, Dvar_coord = {}, {}, {}
-Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=False)
-Dvar_coord_budget = read_netcdf(LnameFiles, Dvar_input_coord_budget, path=path, removeHALO=False)
+Dvar = read_netcdf(LnameFiles, Dvar_input, path=path, removeHALO=True)
+Dvar_coord_budget = read_netcdf(LnameFiles, Dvar_input_coord_budget, path=path, removeHALO=True)
 Dvar_coord = read_netcdf(LnameFiles, Dvar_input_coord, path=path, removeHALO=True)
 
-Dvar['f1']['altitude'], Dvar['f1']['ni_2D'] = comp_altitude1DVar(Dvar['f1'][('CIFNFREE01','AVEF')], Dvar_coord['f1']['ZS'],Dvar_coord['f1']['ZTOP'], Dvar_coord_budget['f1']['cart_level'],Dvar_coord_budget['f1']['cart_ni'])
+Dvar['f1']['altitude'], Dvar['f1']['ni_2D'] = comp_altitude1DVar(Dvar['f1'][('/Budgets/CIFNNUCL01','AVEF')], Dvar_coord['f1']['ZS'],
+    Dvar_coord['f1']['ZTOP'], Dvar_coord_budget['f1']['cart_level'],Dvar_coord_budget['f1']['cart_ni'])
 
 ################################################################
 #########          PANEL 1
 ###############################################################
 Panel1 = PanelPlot(2,3, [25,14],'014_LIMA', titlepad=25, minmaxpad=1.04, timepad=-0.07, colorbarpad=0.03, labelcolorbarpad = 13, colorbaraspect=40)
 
-Lplot = [Dvar['f1'][('RI','AVEF')],Dvar['f1'][('CICE','AVEF')], Dvar['f1'][('RS','AVEF')],
-         Dvar['f1'][('RG','AVEF')],Dvar['f1'][('CIFNFREE01','AVEF')], Dvar['f1'][('CIFNNUCL01','AVEF')]]
+Lplot = [Dvar['f1'][('/Budgets/RI','AVEF')],Dvar['f1'][('/Budgets/CICE','AVEF')], Dvar['f1'][('/Budgets/RS','AVEF')],
+         Dvar['f1'][('/Budgets/RG','AVEF')],Dvar['f1'][('/Budgets/CIFNFREE01','AVEF')], Dvar['f1'][('/Budgets/CIFNNUCL01','AVEF')]]
 LaxeX = [Dvar['f1']['ni_2D']]*len(Lplot)
 LaxeZ = [Dvar['f1']['altitude']]*len(Lplot)
 Ltitle = ['Ice water content', 'Ice concentration', 'Snow water content',
diff --git a/src/LIB/Python/Panel_Plot.py b/src/LIB/Python/Panel_Plot.py
index 24134085e04ad816580a545bfb65de987ec3db8b..8e095faf07fd46b6dcd4f01df5dddca6b0b5759a 100644
--- a/src/LIB/Python/Panel_Plot.py
+++ b/src/LIB/Python/Panel_Plot.py
@@ -1,9 +1,12 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Created on Wed Feb 24 10:49:45 2021
+MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+MNH_LIC for details. version 1.
 
-@author: rodierq
+@author: 07/2021 Quentin Rodier
 """
 import matplotlib as mpl
 mpl.use('Agg')
@@ -11,6 +14,7 @@ import matplotlib.pyplot as plt
 from matplotlib import cm
 from matplotlib.colors import ListedColormap
 import numpy as np
+import cartopy
 import cartopy.feature as cfeature
 
 class PanelPlot():
@@ -56,8 +60,12 @@ class PanelPlot():
         # Grid lines and labels
       if 'PlateCarree' in str(projo):
         gl = ax.gridlines(crs=self.projo, draw_labels=True, linewidth=1, color='gray')
-        gl.top_labels = False
-        gl.right_labels = False
+        if float(cartopy.__version__[:4]) >= 0.18:
+          gl.top_labels = False
+          gl.right_labels = False
+        else:
+          gl.xlabels_top = False
+          gl.ylabels_right = False
         
         #  Coastlines
       if self.drawCoastLines and 'GeoAxes' in str(type(ax)):
@@ -188,7 +196,7 @@ class PanelPlot():
                  Lfacconv=[], ax=[], Lid_overlap=[], colorbar=True, orog=[], Lxlim=[], Lylim=[], Ltime=[], Lpltype=[], LaddWhite_cm=[]):
       """
         Horizontal cross section plot
-        Arguments :
+        Parameters :
             - Lxx    : List of x or y coordinate variable or time axis
             - Lzz    : List of z coordinates variable
             - Lvar   : List of variables to plot
@@ -253,7 +261,7 @@ class PanelPlot():
         
         #  Print time validity
         if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
-        
+
         # Number of contours level
         if not Lstep[i]: #  Default value of number of steps is 20
             Lstep[i] = (Lmaxval[i] - Lminval[i])/20  
@@ -315,7 +323,7 @@ class PanelPlot():
                  Lylim=[], Ltime=[], LaxisColor=[]):
       """
         XY (multiple)-lines plot
-        Arguments :
+        Parameters :
             - Lxx    : List of variables to plot or coordinates along the X axis #TODO : ajouter Lfacconv pour les deux axes. Impact tous les cas test avec lignes X/Y
             - Lyy    : List of variables to plot or coordinates along the Y axis
             - Lxlab  : List of x-axis label
@@ -405,7 +413,7 @@ class PanelPlot():
                  Lid_overlap=[], colorbar=True, Ltime=[], LaddWhite_cm=[], Lpltype=[], Lcbformatlabel=[]):
       """
         Horizontal cross section plot
-        Arguments :
+        Parameters :
             - lon    : longitude 2D array
             - lat    : latitude 2D array
             - Lvar   : List of variables to plot
@@ -554,7 +562,7 @@ class PanelPlot():
                 Lylim=[], Lxlim=[]):
       """
         Horizontal vectors lines
-        Arguments :
+        Parameters :
             - Lxx    : List of x or y coordinate variable (lat or ni or nm)
             - Lyy    : List of y coordinates variable (lon or level)
             - Lvar1   : List of wind-component along x/y or oblic axis (3D for hor. section, 2D for vertical section)
@@ -667,7 +675,7 @@ class PanelPlot():
                 Lylim=[], Lxlim=[]):
       """
         Wind stream lines
-        Arguments :
+        Parameters :
             - Lxx    : List of x or y coordinate variable (lat or ni or nm)
             - Lyy    : List of y coordinates variable (lon or level)
             - Lvar1   : List of wind-component along x/y or oblic axis (3D for hor. section, 2D for vertical section)
@@ -768,7 +776,7 @@ class PanelPlot():
                  Lylim=[], Ltime=[], LaxisColor=[]):
       """
         XY Histogram
-        Arguments :
+        Parameters :
             - Lbins    : List of bins
             - Lvar   : List of the value for each bin
             - Lxlab  : List of x-axis label
diff --git a/src/LIB/Python/misc_functions.py b/src/LIB/Python/misc_functions.py
index 6b2c0932742f0d47e6607275a890f990a50226ce..aa91aa7a311dc721c17f3f28cc7a8fd4f78df5b1 100644
--- a/src/LIB/Python/misc_functions.py
+++ b/src/LIB/Python/misc_functions.py
@@ -1,9 +1,12 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Created on Thu Feb 25 11:25:31 2021
+MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+MNH_LIC for details. version 1.
 
-@author: rodierq
+@author: 07/2021 Quentin Rodier
 """
 import copy
 from scipy.interpolate import RectBivariateSpline
@@ -11,8 +14,6 @@ import numpy as np
 import math
 
 def convert_date(datesince, time_in_sec):
-  print(type(datesince))
-  print(type(time_in_sec))
   return str(time_in_sec) + datesince[:33]
 
 class mean_operator():
@@ -37,103 +38,172 @@ class mean_operator():
             out[k,:,:] = (var[k,:,:] + var[k+1,:,:])*0.5
         return out
 
-def windvec_verti_proj(u, v, lvl, angle):
-  """
-     Compute the projected horizontal wind vector on an axis with a given angle w.r.t. the x/ni axes (West-East)
-     Arguments :
-            - u : U-wind component
-            - v : V-wind component
-            - angle (radian) of the new axe w.r.t the x/ni axes (West-East). angle = 0 for (z,x) sections, angle=pi/2 for (z,y) sections
-     Returns :
-            - a 3D wind component projected on the axe to be used with Panel_Plot.vector as Lvar1
-  """ 
-  out = copy.deepcopy(u)
-  for k in range(len(lvl)):
-      out[k,:,:] = u[k,:,:]*math.cos(angle) + v[k,:,:]*math.sin(angle)
-  return out
+def windvec_verti_proj(u, v, level, angle):
+    """Compute the projected horizontal wind vector on an axis with a given angle w.r.t. the x/ni axes (West-East)
+    
+    Parameters
+    ----------
+    u : array 3D
+        U-wind component
+    
+    v : array 3D
+        V-wind component
+
+    level : array 1D
+        level dimension array
+        
+    angle : float
+        angle (radian) of the new axe w.r.t the x/ni axes (West-East). angle = 0 for (z,x) sections, angle=pi/2 for (z,y) sections
+    
+    Returns
+    -------
+    
+    projected_wind : array 3D
+        a 3D wind component projected on the axe to be used with Panel_Plot.pvector as Lvar1
+    """ 
+    projected_wind = copy.deepcopy(u)
+    for k in range(len(level)):
+        projected_wind[k,:,:] = u[k,:,:]*math.cos(angle) + v[k,:,:]*math.sin(angle)
+    return projected_wind
 
 def oblique_proj(var, ni, nj, lvl, i_beg, j_beg, i_end, j_end):
-  """
-     Compute an oblique projection of a 3D variable w.r.t. its axes
-     Arguments :
-            - var          : 3D var to project (example THT)
-            - ni           : 1D x-axe of the 3D dimension
-            - nj           : 1D y-axe of the 3D dimension
-            - level        : 1D level variable (level or level_w)
-            - i_beg, j_beg : coordinate of the begin point of the new axe
-            - i_end, j_end : coordinate of the end point of the new axe
-     Returns :
-            - a 2D (z,m) variable projected on the oblique axe
-            - a 1D m new axe (distance from the beggining point)
-            - the angle (radian) of the new axe w.r.t the x/ni axes (West-East)
-  """
-  dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) #  Distance de la section oblique  m
-  out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
-  axe_m = np.zeros(int(dist_seg)+1) #Axe des abscisses qui sera tracé selon la coupe
-  axe_m_coord = [] #Coordonnées x,y des points qui composent l'axe
-  axe_m_coord.append( (ni[i_beg],nj[j_beg]) ) #Le premier point est celui donné par l'utilisateur
-  for m in range(int(dist_seg)): #Discrétisation selon distance de la coupe / int(distance_de_la_coupe)
-    axe_m_coord.append( (axe_m_coord[0][0] + (ni[i_end]-ni[i_beg])/(int(dist_seg))*(m+1), 
-                         axe_m_coord[0][1] + (nj[j_end]-nj[j_beg])/(int(dist_seg))*(m+1) ))
-    axe_m[m+1] = np.sqrt((ni[i_beg]-axe_m_coord[m+1][0])**2 + (nj[j_beg]-axe_m_coord[m+1][1])**2)
-  
-  for k in range(len(lvl)):
-    a=RectBivariateSpline(ni, nj,var[k,:,:],kx=1,ky=1) #Interpolation par niveau à l'ordre 1 pour éviter des valeurs négatives de champs strictement > 0
-    for m in range(int(dist_seg)+1):
-      out_var[k,m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1]) # La fonction ev de RectBivariate retourne la valeur la plus proche du point considéré
-      
-  angle_proj = math.acos((ni[i_end]-ni[i_beg])/axe_m[-1])
-  return angle_proj, out_var, axe_m
+    """Compute an oblique projection of a 3D variable w.r.t. its axes
+    
+    Parameters
+    ----------
+    var : array 3D
+        the 3D variable to project (e.g. THT)
+    
+    ni : array 1D
+        1D x-axis of the 3D dimension
+    
+    nj : array 1D
+        1D y-axis of the 3D dimension
+    
+    level : array 1D
+        1D z-axe of the 3D dimension
+    
+    i_beg, j_beg : int
+        coordinate of the begin point of the new axe
+    
+    i_end, j_end : int
+        coordinate of the end point of the new axe       
+    
+     Returns
+     -------
+    angle_proj : float
+        the angle (radian) of the new axe w.r.t the x/ni axes (West-East)
+         
+    out_var : array 2D
+        a 2D (z,m) variable projected on the oblique axe
+    
+    axe_m : array 1D
+        a 1D m new axe (distance from the beggining point)
+    
+    """
+    dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) #  Distance de la section oblique  m
+    out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
+    axe_m = np.zeros(int(dist_seg)+1) #Axe des abscisses qui sera tracé selon la coupe
+    axe_m_coord = [] #Coordonnées x,y des points qui composent l'axe
+    axe_m_coord.append( (ni[i_beg],nj[j_beg]) ) #Le premier point est celui donné par l'utilisateur
+    for m in range(int(dist_seg)): #Discrétisation selon distance de la coupe / int(distance_de_la_coupe)
+        axe_m_coord.append( (axe_m_coord[0][0] + (ni[i_end]-ni[i_beg])/(int(dist_seg))*(m+1), 
+                             axe_m_coord[0][1] + (nj[j_end]-nj[j_beg])/(int(dist_seg))*(m+1) ))
+        axe_m[m+1] = np.sqrt((ni[i_beg]-axe_m_coord[m+1][0])**2 + (nj[j_beg]-axe_m_coord[m+1][1])**2)
+    
+    for k in range(len(lvl)):
+        a=RectBivariateSpline(ni, nj,var[k,:,:],kx=1,ky=1) #Interpolation par niveau à l'ordre 1 pour éviter des valeurs négatives de champs strictement > 0
+        for m in range(int(dist_seg)+1):
+            out_var[k,m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1]) # La fonction ev de RectBivariate retourne la valeur la plus proche du point considéré
+    
+    angle_proj = math.acos((ni[i_end]-ni[i_beg])/axe_m[-1])
+    return angle_proj, out_var, axe_m
 
 def comp_altitude1DVar(oneVar2D, orography, ztop, level, n_xory):
-  """
-     Compute and returns an altitude and x or y grid mesh variable in 2D following the topography in 1D
-     To be used with 2D simulations
-     Arguments :
-            - oneVar2D  : a random netCDF 2D var (example UT, THT)
-            - orography : 1D orography (ZS)
-            - ztop      : scalar of the top height of the model (ZTOP)
-            - level     : 1D level variable (level or level_w)
-            - n_xory:   : 1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
-     Returns :
-            - a 2D altitude variable with topography taken into account
-            - a 2D directionnal variable
-  """
-  n_xory_2D = copy.deepcopy(oneVar2D)
-  altitude  =  copy.deepcopy(oneVar2D)
-  for i in range(len(level)):
-    n_xory_2D[i,:] = n_xory
-  for j in range(len(n_xory)):
+    """Compute and returns an altitude and x or y grid mesh variable in 2D following the topography in 1D
+    To be used with 2D simulations
+    
+    Parameters
+    ----------
+    oneVar2D : array 2D
+        a 2D array (e.g. UT, THT)
+        
+    orography : array 1D
+        1D orography (ZS)
+        
+    ztop : real
+        scalar of the top height of the model (ZTOP)
+        
+    level : array 1D
+        1D level variable (level or level_w)
+        
+    n_xory : array 1D  
+        1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
+        
+    Returns
+    -------
+    altitude
+        a 2D altitude variable with topography taken into account
+        
+    n_xory_2D
+        a 2D directionnal variable duplicated from n_xory
+    """
+    n_xory_2D = copy.deepcopy(oneVar2D)
+    altitude  =  copy.deepcopy(oneVar2D)
+    
     for k in range(len(level)):
-        altitude[k,j] = orography[j] + level[k]*((ztop-orography[j])/ztop)
-  return altitude, n_xory_2D
+        n_xory_2D[k,:] = n_xory
+        for j in range(len(n_xory)):
+            for k in range(len(level)):
+                altitude[k,j] = orography[j] + level[k]*((ztop-orography[j])/ztop)
+    return altitude, n_xory_2D
 
 def comp_altitude2DVar(oneVar3D, orography, ztop, level, n_y, n_x):
-  """
-     Compute and returns an altitude and x or y grid mesh variable in 3D following the topography in 2D
-     To be used with 3D simulations in cartesian coordinates
-     Arguments :
-            - oneVar3D  : a random netCDF 3D var (example UT, THT)
-            - orography : 2D orography (ZS)
-            - ztop      : scalar of the top height of the model (ZTOP)
-            - level     : 1D level variable (level or level_w)
-            - n_xory:   : 1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
-     Returns :
-            - a 3D altitude variable with topography taken into account
-            - a 3D directionnal variable
-  """
-  n_x3D = copy.deepcopy(oneVar3D)
-  n_y3D = copy.deepcopy(oneVar3D)
-  altitude  =  copy.deepcopy(oneVar3D)
-  for i in range(len(level)):
-    n_y3D[i,:] = n_y
-    n_x3D[i,:] = n_x
-  for i in range(oneVar3D.shape[2]):
-    for j in range(oneVar3D.shape[1]):
-      if ztop==0:
-        altitude[:,i,j] = level[:]
-      else:
-        for k in range(len(level)):
-            altitude[k,j,i] = orography[j,i] + level[k]*((ztop-orography[j,i])/ztop)
-
-  return altitude, n_x3D, n_y3D
+    """Compute and returns an altitude and x or y grid mesh variable in 3D following the topography in 2D
+    To be used with 3D simulations
+    
+    Parameters
+    ----------
+    oneVar3D : array 3D
+        a 3D array (e.g. UT, THT)
+        
+    orography : array 2D
+        2D orography (ZS)
+        
+    ztop : real
+        scalar of the top height of the model (ZTOP)
+        
+    level : array 1D
+        1D level variable (level or level_w)
+        
+    n_x : array 1D  
+        1D directionnal grid variable along i (ni_u, or ni_v)
+        
+    n_y : array 1D  
+        1D directionnal grid variable along j (nj_u, or nj_v)
+        
+    Returns
+    -------
+    altitude
+        a 3D altitude variable with topography taken into account
+        
+    n_x3D
+        a 3D directionnal variable duplicated from n_x
+        
+    n_y3D
+        a 3D directionnal variable duplicated from n_y
+    """
+    n_x3D = copy.deepcopy(oneVar3D)
+    n_y3D = copy.deepcopy(oneVar3D)
+    altitude  =  copy.deepcopy(oneVar3D)
+    for i in range(len(level)):
+        n_y3D[i,:] = n_y
+        n_x3D[i,:] = n_x
+    for i in range(oneVar3D.shape[2]):
+        for j in range(oneVar3D.shape[1]):
+            if ztop==0:
+                altitude[:,i,j] = level[:]
+            else:
+                for k in range(len(level)):
+                    altitude[k,j,i] = orography[j,i] + level[k]*((ztop-orography[j,i])/ztop)
+    return altitude, n_x3D, n_y3D
\ No newline at end of file
diff --git a/src/LIB/Python/read_MNHfile.py b/src/LIB/Python/read_MNHfile.py
index 7ad0bd539810b5413023b7339b11e479f3a7d548..e9a4ce4397d923c009d93d93cb2864d7061dd8d9 100644
--- a/src/LIB/Python/read_MNHfile.py
+++ b/src/LIB/Python/read_MNHfile.py
@@ -1,287 +1,344 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Created on Mon Feb 22 10:29:13 2021
+MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+MNH_LIC for details. version 1.
 
-@author: rodierq
+@author: 07/2021 Quentin Rodier
 """
-
 import netCDF4 as nc
 import numpy as np
 
-def read_netcdf(LnameFiles, Dvar_input, path='.', removeHALO=True):
-  Dvar = {}
-  for i,nameFiles in enumerate(LnameFiles):
-    f_nb = 'f' + str(i+1)
-    print('Reading file ' + f_nb)
-    print(path + nameFiles)
-    if '000' in nameFiles[-6:-3]: #time series file (diachronic)
+def read_netcdf(LnameFiles, Dvar_input, path='.', get_data_only=True, del_empty_dim=True, removeHALO=True):
+    """Read a netCDF4 Meso-NH file
+    For each file, call functions to read diachronic or synchronous file
+    
+    Parameters
+    ----------
+    LnameFiles : list of str
+        list of Meso-NH netCDF4 file (diachronic or synchronous)
+    
+    Dvar_input : Dict{'fileNumber' : 'var_name',('group_name','var_name')}
+        where
+        'fileNumber' is a str corresponding to 'f' + the file number in LnameFiles (by order)
+        'var_name' is the exact str of the netCDF4 variable name
+        ('group_name','var_name') is the exact tuple of the (sub-)groups name and the netCDF4 variable name
+        e.g. : {'f1':['ZS', 'WT','ni', 'level'],
+                'f2':[('/LES_budgets/Cartesian/Not_time_averaged/Not_normalized/cart/',MEAN_TH'),('/Budgets/RI','AVEF')]
+                }
+    
+    path : str
+        unique path of the files
+    
+    get_data_only : bool, default: True
+        if True,  the function returns Dvar as masked_array (only data)
+        if False, the function returns Dvar as netCDF4._netCDF4.Variable
+    
+    del_empty_dim : bool, default: True
+        if get_data_only=True and del_empty_dim=True, returns Dvar as an array without dimensions with size 1 and 0
+        e.g. : an array of dimensions (time_budget, cart_level, cart_nj, cart_ni) with shape (180,1,50,1) is returned (180,50)
+        
+    removeHALO : bool, default: True
+        if True, remove first and last (NHALO=1) point [1:-1] if get_data_only=True on each 
+        level, level_w, ni, ni_u, ni_v, nj, nj_u, nj_v dimensions
+           
+    Returns
+    -------
+    Dvar : Dict 
+        Dvar[ifile]['var_name']                if the group contains only one variable
+        Dvar[ifile][('group_name','var_name')] if the group contains more than one variable
+    """
+    Dvar = {}
+    for i,nameFiles in enumerate(LnameFiles):
+        f_nb = 'f' + str(i+1)
+        print('Reading file ' + f_nb)
+        print(path + nameFiles)
         theFile = nc.Dataset(path + nameFiles,'r')
-        if theFile['MASDEV'][0] <= 54:
-            read_TIMESfiles_54(theFile, f_nb, Dvar_input, Dvar)
-        else : # 55 groups variables
-            read_TIMESfiles_55(theFile, f_nb, Dvar_input, Dvar, removeHALO)
+        Dvar[f_nb] = {}
+        if '000' in nameFiles[-6:-3]: 
+            if theFile['MASDEV'][0] <= 54:
+                raise TypeError('The python lib is available for MNH >= 5.5')
+            else:
+                Dvar[f_nb] = read_TIMESfiles_55(theFile, Dvar_input[f_nb], Dvar[f_nb], get_data_only, del_empty_dim, removeHALO)
+        else:
+            Dvar[f_nb]= read_BACKUPfile(theFile, Dvar_input[f_nb], Dvar[f_nb], get_data_only, del_empty_dim, removeHALO)
         theFile.close()
-    else:
-        read_BACKUPfile(nameFiles, f_nb, Dvar_input, Dvar, path=path, removeHALO=removeHALO)
-  return Dvar #Return the dic of [files][variables]
-
-def read_BACKUPfile(nameF, ifile, Dvar_input, Dvar_output, path='.', removeHALO=True):
-  theFile = nc.Dataset(path + nameF,'r')
-  Dvar_output[ifile] = {} #initialize dic for each files 
-       
-  #  Reading date since beginning of the model run
-  Dvar_output[ifile]['time'] =  theFile.variables['time'][0]
-  Dvar_output[ifile]['date'] = nc.num2date(Dvar_output[ifile]['time'],units=theFile.variables['time'].units, calendar = theFile.variables['time'].calendar)
+    return Dvar
 
-  for var in Dvar_input[ifile]: #For each files
-    #  Read variables
-    n_dim = theFile.variables[var].ndim
-    name_dim = theFile.variables[var].dimensions
+def read_var(theFile, Dvar, var_name, get_data_only=True, del_empty_dim=True, removeHALO=True):
+    """Read a netCDF4 variable
+    
+    Parameters
+    ----------
+    theFile : netCDF4._netCDF4.Dataset
+        a Meso-NH diachronic netCDF4 file
+        
+    var_name : str
+        a Meso-NH netCDF4 variable name
+        
+    get_data_only : bool, default: True
+        if True,  the function returns Dvar as masked_array (only data)
+        if False, the function returns Dvar as netCDF4._netCDF4.Variable
     
-    if (n_dim ==0) or (n_dim == 1 and 'time' in name_dim):  #Scalaires or Variable time
-      Dvar_output[ifile][var] = theFile.variables[var][0].data
-    else:      
+    del_empty_dim : bool, default: True
+        if get_data_only=True and del_empty_dim=True, returns Dvar as an array without dimensions with size 1 and 0
+        e.g. : an array of dimensions (time_budget, cart_level, cart_nj, cart_ni) with shape (180,1,50,1) is returned (180,50)
+        
+    removeHALO : bool, default: True
+        if True, remove first and last (NHALO=1) point [1:-1] if get_data_only=True on each 
+        level, level_w, ni, ni_u, ni_v, nj, nj_u, nj_v dimensions
+           
+    Returns
+    -------
+    Dvar : Dict 
+        Dvar['var_name']                if the group contains only one variable
+        Dvar[('group_name','var_name')] if the group contains more than one variable
+    """
+    try:
+        var_dim = theFile.variables[var_name].ndim
+        var_dim_name = theFile.variables[var_name].dimensions
+    except:
+        raise KeyError("Group and variable name not found in the file, check the group/variable name with ncdump -h YourMNHFile.000.nc. You asked for variable : " + var_name)
+    
+    if not get_data_only:
+        Dvar[var_name] = theFile.variables[var_name]
+    else: 
+        if var_dim == 0:
+            Dvar[var_name] = theFile.variables[var_name][0].data
+        elif var_dim == 1:
+            Dvar[var_name] = theFile.variables[var_name][:]
+        elif var_dim == 2:
+            Dvar[var_name] = theFile.variables[var_name][:,:]
+        elif var_dim == 3:
+            Dvar[var_name] = theFile.variables[var_name][:,:,:]
+        elif var_dim == 4:
+            Dvar[var_name] = theFile.variables[var_name][:,:,:,:]
+        elif var_dim == 5:
+            Dvar[var_name] = theFile.variables[var_name][:,:,:,:,:]
+        elif var_dim == 6:
+            Dvar[var_name] = theFile.variables[var_name][:,:,:,:,:,:]
+        elif var_dim == 7:
+            Dvar[var_name] = theFile.variables[var_name][:,:,:,:,:,:,:]
         if removeHALO:
-          if n_dim == 1:
-            Dvar_output[ifile][var] = theFile.variables[var][1:-1] #Variables 1D
-          elif n_dim == 2:
-            if theFile.variables[var].shape[0] == 1 and 'size' in name_dim[1]: #Variables 2D with the second dimension not a coordinate (--> list of strings : chemicals)
-                Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
-            elif theFile.variables[var].shape[0] == 1: #Variables 2D with first dim = 0
-                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1] #Variables 2D
-            else:
-              Dvar_output[ifile][var] = theFile.variables[var][1:-1,1:-1] #Variables 2D
-          elif n_dim == 5: #Variables time, sizeXX, level, nj, ni (ex: chemical budgets in 3D)
-              Dvar_output[ifile][var] = theFile.variables[var][0, :, 1:-1,:1:-1,1:-1]
-          elif n_dim == 4 and 'time' in name_dim and ('level' in name_dim or 'level_w' in name_dim): # time,z,y,x
-              if theFile.variables[var].shape[1] == 1: #Variables 4D time,z,y,x with time=0 z=0
-                Dvar_output[ifile][var] = theFile.variables[var][0,0,1:-1,1:-1] #Variables 2D y/x
-              elif theFile.variables[var].shape[2] == 1: #Variable 2D  (0,zz,0,xx)
-                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,0,1:-1] #Variables 2D z/y
-              elif theFile.variables[var].shape[3] == 1: #Variable 2D  (0,zz,yy,0)
-                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1,0] #Variables 2D z/x
-              ## ATTENTION VARIABLE 1D codé en 4D non faite
-              else: #Variable 3D simple
-                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1,1:-1] #Variables time + 3D     
-          elif n_dim == 4 and 'time' in name_dim and 'level' not in name_dim and 'level_w' not in name_dim: # time,nb_something,y,x
-               Dvar_output[ifile][var] = theFile.variables[var][0,:,1:-1,1:-1] #Variables 2D y/x
-          elif n_dim == 3 and 'time' in name_dim: # time, y, x 
-            Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1]
-          else :
-            Dvar_output[ifile][var] = theFile.variables[var][1:-1,1:-1,1:-1]  #Variables 3D
-        else:
-          if n_dim == 1:
-            Dvar_output[ifile][var] = theFile.variables[var][:] #Variables 1D  
-          elif n_dim == 2:
-            if theFile.variables[var].shape[0] == 1 and 'size' in name_dim[1]: #Variables 2D with the second dimension not a coordinate (--> list of strings : chemicals)
-              Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
-            elif theFile.variables[var].shape[0] == 1: #Variables 2D with first dim = 0
-              Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
-            else:
-              Dvar_output[ifile][var] = theFile.variables[var][:,:] #Variables 2D
-          elif n_dim == 5: #Variables time, sizeXX, level, nj, ni (ex: chemical budgets in 3D)
-              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,:,:]
-          elif n_dim == 4: # time,z,y,x
-            if theFile.variables[var].shape[1] == 1: #Variables 4D time,z,y,x with time=0 z=0
-              Dvar_output[ifile][var] = theFile.variables[var][0,0,:,:] #Variables 2D y/x
-            elif theFile.variables[var].shape[2] == 1: #Variable 2D  (0,zz,0,xx)
-              Dvar_output[ifile][var] = theFile.variables[var][0,:,0,:] #Variables 2D z/y
-            elif theFile.variables[var].shape[3] == 1: #Variable 2D  (0,zz,yy,0)
-              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,0] #Variables 2D z/x
-            ## ATTENTION VARIABLE 1D codé en 4D non faite
-            else: #Variable 3D simple
-              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,:] #Variables time + 3D
-          elif n_dim ==3 and name_dim in var.dimensions: # time, y, x
-            Dvar_output[ifile][var] = theFile.variables[var][0,:,:]
-          else:
-            Dvar_output[ifile][var] = theFile.variables[var][:,:,:]  #Variables 3D
-        #  For all variables except scalars, change Fill_Value to NaN
-        Dvar_output[ifile][var]= np.where(Dvar_output[ifile][var] != -99999.0, Dvar_output[ifile][var], np.nan)
-        Dvar_output[ifile][var]= np.where(Dvar_output[ifile][var] != 999.0, Dvar_output[ifile][var], np.nan)
-
-  theFile.close()
-  return Dvar_output #Return the dic of [files][variables]
-
-def read_TIMESfiles_54(theFile, ifile, Dvar_input, Dvar_output):
-    Dvar_output[ifile] = {} #initialize dic for each files 
+            for i in range(8):
+                try:
+                    if var_dim_name[i]=='level' or var_dim_name[i]=='level_w' or \
+                    var_dim_name[i]=='ni' or var_dim_name[i]=='ni_u' or var_dim_name[i]=='ni_v' or \
+                    var_dim_name[i]=='nj' or var_dim_name[i]=='nj_u' or var_dim_name[i]=='nj_v':
+                        if var_dim != 0:
+                            Dvar[var_name] = removetheHALO(i+1, Dvar[var_name])
+                except:
+                    break
+        if del_empty_dim:
+            Ldimtosqueeze=[]
+            var_shape = theFile.variables[var_name].shape
+            for i in range(8):
+                try:
+                    if var_shape[i]==1: Ldimtosqueeze.append(i)
+                except IndexError:
+                    break
+            Ldimtosqueeze=tuple(Ldimtosqueeze)
+            Dvar[var_name] = np.squeeze(Dvar[var_name], axis=Ldimtosqueeze) 
+        
+    return Dvar
 
-    #  Level variable is automatically read without the Halo
-    Dvar_output[ifile]['level'] = theFile.variables['level'][1:-1]
+def read_from_group(theFile, Dvar, group_name, var_name, get_data_only=True, del_empty_dim=True,removeHALO=True):
+    """Read a variable from a netCDF4 group 
     
-    #  Time variable is automatically read (time since begging of the run) from the 1st variable of the asked variable to read
-    suffix, name_first_var = remove_PROC(Dvar_input[ifile][0])
-    try: #  It is possible that there is only one value (one time) in the .000 file, as such time series are not suitable and the following line can't be executed. The time variable is then not written
-        increment = theFile.variables[name_first_var+'___DATIM'][1,-1] - theFile.variables[name_first_var+'___DATIM'][0,-1] #-1 is the last entry of the date (current UTC time in seconds)
-        length_time = theFile.variables[name_first_var+'___DATIM'].shape[0]
-        Dvar_output[ifile]['time'] = np.arange(increment,increment*(length_time+1),increment)
+    Parameters
+    ----------
+    theFile : netCDF4._netCDF4.Dataset
+        a Meso-NH diachronic netCDF4 file
+    
+    group_name : str
+        a Meso-NH netCDF4 group name
+        
+    var_name : str
+        a Meso-NH netCDF4 variable name
+        
+    get_data_only : bool, default: True
+        if True,  the function returns Dvar as masked_array (only data)
+        if False, the function returns Dvar as netCDF4._netCDF4.Variable
+    
+    del_empty_dim : bool, default: True
+        if get_data_only=True and del_empty_dim=True, returns Dvar as an array without dimensions with size 1 and 0
+        e.g. : an array of dimensions (time_budget, cart_level, cart_nj, cart_ni) with shape (180,1,50,1) is returned (180,50)
+        
+    removeHALO : bool, default: True
+        if True, remove first and last (NHALO=1) point [1:-1] if get_data_only=True on each 
+        level, level_w, ni, ni_u, ni_v, nj, nj_u, nj_v dimensions
+        
+    Returns
+    -------
+    Dvar : Dict 
+        Dvar['var_name']                if the group contains only one variable
+        Dvar[('group_name','var_name')] if the group contains more than one variable
+    """
+    try:
+        var_dim = theFile[group_name].variables[var_name].ndim
+        var_dim_name = theFile[group_name].variables[var_name].dimensions
     except:
-        pass
+        raise KeyError("Group and variable name not found in the file, check the group/variable name with ncdump -h YourMNHFile.000.nc. You asked for group/variable : " + group_name + var_name)
     
-    for var in Dvar_input[ifile]: #For each files
-        suffix, var_name = remove_PROC(var)
-        n_dim = theFile.variables[var].ndim
-        name_dim = theFile.variables[var].dimensions
+    if not get_data_only:
+        Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name]
+    else:
+        if var_dim == 0:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][0].data
+        if var_dim == 1:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:]
+        elif var_dim == 2:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:]
+        elif var_dim == 3:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:,:]
+        elif var_dim == 4:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:,:,:]
+        elif var_dim == 5:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:,:,:,:]
+        elif var_dim == 6:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:,:,:,:,:]
+        elif var_dim == 7:
+            Dvar[(group_name,var_name)] = theFile[group_name].variables[var_name][:,:,:,:,:,:,:]
+        if removeHALO:
+            for i in range(8):
+                try:
+                    if var_dim_name[i]=='level' or var_dim_name[i]=='level_w' or \
+                    var_dim_name[i]=='ni' or var_dim_name[i]=='ni_u' or var_dim_name[i]=='ni_v' or \
+                    var_dim_name[i]=='nj' or var_dim_name[i]=='nj_u' or var_dim_name[i]=='nj_v':
+                        if var_dim != 0:
+                            Dvar[var_name] = removetheHALO(i+1, Dvar[var_name])
+                except:
+                    break
+        if del_empty_dim:
+            Ldimtosqueeze=[]
+            var_shape = Dvar[(group_name,var_name)].shape
+            for i in range(8):
+                try:
+                    if var_shape[i]==1: Ldimtosqueeze.append(i)
+                except IndexError:
+                    break
+            Ldimtosqueeze=tuple(Ldimtosqueeze)                
+            Dvar[(group_name,var_name)] = np.squeeze(Dvar[(group_name,var_name)], axis=Ldimtosqueeze)
+            
+        # LES budget needs to be transposed to use psection functions without specifying .T each time
+        if 'LES_budgets' in group_name: 
+            Dvar[(group_name,var_name)] = Dvar[(group_name,var_name)].T   
+    return Dvar
 
-        #  First, test if the variable is a dimension/coordinate variable
-        if (n_dim ==0):  #  Scalaires variable
-             Dvar_output[ifile][var] = theFile.variables[var][0].data
-             pass
-        elif n_dim == 1:
-            Dvar_output[ifile][var_name] = theFile.variables[var][1:-1]  #  By default, the Halo is always removed because is not in the other variables in any .000 variable
-            pass
-        elif n_dim == 2:
-            Lsize1 = list_size1(n_dim, name_dim)
-            if Lsize1 == [True, False]: Dvar_output[ifile][var_name] = theFile.variables[var][0,1:-1] 
-            pass
+def read_BACKUPfile(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_dim=True, removeHALO=True):
+    """Read variables from Meso-NH MASDEV >= 5.5.0 synchronous file
+    For all variables in Dvar_input of one file, call functions to read the variable of the group+variable
+    
+    Parameters
+    ----------
+    theFile : netCDF4._netCDF4.Dataset
+        a Meso-NH diachronic netCDF4 file
         
-        Lsize1 = list_size1(n_dim, name_dim)
-        if Lsize1 == [True, False, False, True, True]: Dvar_output[ifile][var_name] = theFile.variables[var][0,:,:,0,0].T # Need to be transposed here
-        if Lsize1 == [True, True, False, True, False]: Dvar_output[ifile][var_name] = theFile.variables[var][0,0,:,0,:]
-
-    return Dvar_output #Return the dic of [files][variables]
+    Dvar_input : Dict{'var_name',('group_name','var_name')}
+        with
+        'var_name' is the exact str of the netCDF4 variable name
+        ('group_name','var_name') is the exact tuple of the (sub-)groups name and the netCDF4 variable name
+        e.g. : {'f1':['ZS', 'WT','ni', 'level'],
+                'f2':[('/LES_budgets/Cartesian/Not_time_averaged/Not_normalized/cart/',MEAN_TH'),('/Budgets/RI','AVEF')]
+                }
+    
+    get_data_only: bool, default: True
+        if True,  the function returns Dvar as masked_array (only data)
+        if False, the function returns Dvar as netCDF4._netCDF4.Variable
+    
+    del_empty_dim: bool, default: True
+        if get_data_only=True and del_empty_dim=True, returns Dvar as masked_array without dimensions with size 1 and 0
+        e.g. : an array of dimensions (time_budget, cart_level, cart_nj, cart_ni) with shape (180,1,50,1) is returned (180,50)
+    
+    Returns
+    -------
+    Dvar : Dict 
+    Dvar['var_name']                if the group contains only one variable
+    Dvar[('group_name','var_name')] if the group contains more than one variable
+    """     
+    #  Reading date since beginning of the model run
+    Dvar['time'] = theFile.variables['time'][0]
+    Dvar['date'] = nc.num2date(Dvar['time'],units=theFile.variables['time'].units, calendar = theFile.variables['time'].calendar)
+           
+    for var in Dvar_input:
+        if type(var) == tuple:
+            Dvar = read_from_group(theFile, Dvar, var[0], var[1], get_data_only, del_empty_dim, removeHALO)
+        else:
+            Dvar = read_var(theFile, Dvar, var, get_data_only, del_empty_dim, removeHALO)
+            
+            #  For all variables except scalars, change Fill_Value to NaN
+            Dvar[var]= np.where(Dvar[var] != -99999.0, Dvar[var], np.nan)
+            Dvar[var]= np.where(Dvar[var] != 999.0, Dvar[var], np.nan)
+    return Dvar
 
-def read_TIMESfiles_55(theFile, ifile, Dvar_input, Dvar_output, removeHALO=False):
-    """
-        Read variables from MNH MASDEV >= 5.5.0 
-        Parameters :
-            - Dvar_input : dictionnary of {file : var}. var can be either 
-                - a string = the variable name 
-                - or a tuple of ('group_name','var_name')
-                If the variable desired is in a group_name and the group_name is not specified, it is assumed group_name = variable_name
-            except for specific variable as (cart, neb, clear, cs1, cs2, cs3) type
-        Return :
-        Dvar_output : dictionnary of Dvar_output[ifile][variables or tuple (group,variables) if the user specified a tuple]
-    """
-    Dvar_output[ifile] = {} #initialize dic for each files 
-    def read_var(theFile, Dvar_output, var):
-        suffix, var_name = remove_PROC(var)
-        try: #  NetCDF4 Variables
-            n_dim = theFile.variables[var_name].ndim
-            #  First, test if the variable is a dimension/coordinate variable
-            if (n_dim ==0):  #  Scalaires variable
-                Dvar_output[var_name] = theFile.variables[var_name][0].data
-            else:
-                if(removeHALO):
-                    if n_dim == 1:
-                        Dvar_output[var_name] = theFile.variables[var_name][1:-1]
-                    elif n_dim == 2:
-                        Dvar_output[var_name] = theFile.variables[var_name][1:-1,1:-1] 
-                    else: 
-                        raise NameError("Lecture des variables de dimension sup a 2 pas encore implementees pour fichiers .000")
-                else:
-                    if n_dim == 1:
-                        Dvar_output[var_name] = theFile.variables[var_name][:]
-                    elif n_dim == 2:
-                        Dvar_output[var_name] = theFile.variables[var_name][:,:] 
-                    else: 
-                        raise NameError("Lecture des variables de dimension sup a 2 pas encore implementees pour fichiers .000")
-        except KeyError: # NetCDF4 Group not specified by the user
-            if '(cart)' in var_name or '(neb)' in var_name or '(clear)' in var_name or '(cs1)' in var_name or '(cs2)' in var_name or '(cs3)' in var_name or 'AVEF' in suffix or 'INIF' in suffix or 'ENDF' in suffix:
-            # If users specify the complete variable name with averaging type
-              group_name = get_group_from_varname(var_name)
-            else:
-              group_name = var_name
-            read_from_group(theFile, Dvar_output, group_name, var)
-        return Dvar_output
+def read_TIMESfiles_55(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_dim=True, removeHALO=True):
+    """Read variables from Meso-NH MASDEV >= 5.5.0 diachronic file
+    For all variables in Dvar_input of one file, call functions to read the variable of the group+variable
 
-    def read_from_group(theFile, Dvar_output, group_name, var):
-        """
-        Read variables from MNH MASDEV >= 5.5.0 
-        Parameters :
-            - var : the variable name
-            - group_name : the group name            
-        Return :
-        Dvar_output : dictionnary of :
-            - Dvar_output[ifile]['var_name']                if the group contains only one variable
-            - Dvar_output[ifile][('group_name','var_name']  if the group contains more than one variable
-        """
-        if '___' in var:
-            suffix, var_name = remove_PROC(var)
-        else:
-            suffix = var
-            var_name = var
-        if group_name == 'TSERIES' or group_name =='AVION': #always 1D
-            Dvar_output[(group_name,var)] = theFile.groups[group_name].variables[var][:]
-        elif group_name == 'ZTSERIES' or group_name =='AVIONZ': #always 2D 
-            Dvar_output[(group_name,var)] = theFile.groups[group_name].variables[var][:,:].T
-        elif 'XTSERIES' in group_name: #always 2D
-            Dvar_output[(group_name,var)] = theFile.groups[group_name].variables[var][:,:].T
-        elif theFile.groups[group_name].type == 'TLES' : #  LES type
-            try: #By default, most variables read are 2D cart and user does not specify it in the variable name
-              whites = ' '*(17 - len('(cart)') - len(var_name))
-              Dvar_output[var] = theFile.groups[var].variables[var + whites + '(cart)'][:,:].T
-            except:
-              try: #Variable 3D sv,time_les, level_les
-                  Dvar_output[var] = theFile.groups[group_name].variables[var][:,:,:]
-              except:
-                try: #Variable 2D with type of variable specified (cart, neb, clear, cs1, cs2, cs3) 
-                  Dvar_output[var] = theFile.groups[group_name].variables[var][:,:].T
-                except ValueError: #  Variable 1D
-                  Dvar_output[var] = theFile.groups[group_name].variables[var][:]
-        elif theFile.groups[group_name].type == 'CART':  #  Budget CART type
-            shapeVar = theFile.groups[group_name].variables[suffix].shape
-            Ltosqueeze=[] #  Build a tuple with the number of the axis which are 0 dimensions to be removed by np.squeeze
-            if shapeVar[0]==1: Ltosqueeze.append(0)
-            if shapeVar[1]==1: Ltosqueeze.append(1)
-            if shapeVar[2]==1: Ltosqueeze.append(2)
-            if shapeVar[3]==1: Ltosqueeze.append(3)
-            Ltosqueeze=tuple(Ltosqueeze)
-            if len(theFile.groups[group_name].variables.keys()) > 1: # If more than one variable in the group
-              Dvar_output[(group_name,var)] = np.squeeze(theFile.groups[group_name].variables[suffix][:,:,:,:], axis=Ltosqueeze) 
-            else:
-              Dvar_output[group_name] = np.squeeze(theFile.groups[group_name].variables[suffix][:,:,:,:], axis=Ltosqueeze) 
-        elif theFile.groups[group_name].type == 'MASK':  #  Budget MASK type
-            shapeVar = theFile.groups[group_name].variables[suffix].shape
-            Ltosqueeze=[] #  Build a tuple with the number of the axis which are 0 dimensions to be removed by np.squeeze
-            if shapeVar[0]==1: Ltosqueeze.append(0)
-            if shapeVar[1]==1: Ltosqueeze.append(1)
-            if shapeVar[2]==1: Ltosqueeze.append(2)
-            Ltosqueeze=tuple(Ltosqueeze)
-            if len(theFile.groups[group_name].variables.keys()) > 1: # If more than one variable in the group
-              Dvar_output[(group_name,var)] = np.squeeze(theFile.groups[group_name].variables[suffix][:,:,:], axis=Ltosqueeze).T
-            else:
-              Dvar_output[group_name] = np.squeeze(theFile.groups[group_name].variables[suffix][:,:,:], axis=Ltosqueeze).T
-        else:
-            raise NameError("Type de groups variables not implemented in read_MNHfile.py")
-        return Dvar_output
-    for var in Dvar_input[ifile]: #For each var
+    Parameters
+    ----------
+    theFile : netCDF4._netCDF4.Dataset
+        a Meso-NH diachronic netCDF4 file
+        
+    Dvar_input : Dict{'var_name',('group_name','var_name')}
+        with
+        'var_name' is the exact str of the netCDF4 variable name
+        ('group_name','var_name') is the exact tuple of the (sub-)groups name and the netCDF4 variable name
+        e.g. : {'f1':['ZS', 'WT','ni', 'level'],
+                'f2':[('/LES_budgets/Cartesian/Not_time_averaged/Not_normalized/cart/',MEAN_TH'),('/Budgets/RI','AVEF')]
+                }
+    
+    get_data_only: bool, default: True
+        if True,  the function returns Dvar as masked_array (only data)
+        if False, the function returns Dvar as netCDF4._netCDF4.Variable
+    
+    del_empty_dim: bool, default: True
+        if get_data_only=True and del_empty_dim=True, returns Dvar as masked_array without dimensions with size 1 and 0
+        e.g. : an array of dimensions (time_budget, cart_level, cart_nj, cart_ni) with shape (180,1,50,1) is returned (180,50)
+    
+    Returns
+    -------
+    Dvar : Dict 
+    Dvar[ifile]['var_name']                if the group contains only one variable
+    Dvar[ifile][('group_name','var_name')] if the group contains more than one variable
+    """  
+    for var in Dvar_input:
+        print(var)
         if type(var) == tuple:
-            Dvar_output[ifile] = read_from_group(theFile, Dvar_output[ifile], var[0], var[1])
+            Dvar = read_from_group(theFile, Dvar, var[0], var[1], get_data_only, del_empty_dim, removeHALO)
         else:
-            Dvar_output[ifile] = read_var(theFile, Dvar_output[ifile], var)
-    
-    return Dvar_output #Return the dic of [files][variables]
+            Dvar = read_var(theFile, Dvar, var, get_data_only, del_empty_dim, removeHALO)
+    return Dvar
 
-def list_size1(n_dim, named_dim):
-    Lsize1 = []
-    for i in range(n_dim):
-        if 'size1' == named_dim[i]:
-            Lsize1.append(True)
-        else:
-            Lsize1.append(False)
-    return Lsize1
+def removetheHALO(idim, var):
+    """Remove a NHALO=1 point [1:-1] at a given dimension idim of a variable var
     
-def remove_PROC(var):
-    if '___PROC' in var:
-        var_name = var[:-8]
-        suffix = "" # No need of suffix for MNHVERSION < 550 (suffix is for NetCDF4 group)
-    elif  '___ENDF' in var or '___INIF' in var or '___AVEF' in var:
-        var_name = var[:-7]
-        suffix = var[-4:]
-    else:
-        var_name = var
-        suffix = ''
-    return suffix, var_name
-
-def get_group_from_varname(var):
-    group_name=''
-    if '___ENDF' in var or '___INIF' in var or '___AVEF' in var: #Variable CART
-        suff, group_name = remove_PROC(var)
-    else: #Variable with whites in names ex: 'MEAN_TH     (cart)'
-        for i in range(len(var)):
-          if var[i] != ' ':
-            group_name+=var[i]
-          else: # As soon as the caracter is a blank, the group variable is set
-            break
-    return group_name 
+    Parameters
+    ----------
+    idim: int
+        the dimension over which remove the first and last point
+    
+    var: array
+        a Meso-NH netCDF4 variable name
+    
+    Returns
+    -------
+    var : array 
+    """  
+    if idim == 1:
+        var = var[1:-1]
+    elif idim == 2:
+        var = var[:,1:-1]
+    elif idim == 3:
+        var = var[:,:,1:-1]
+    elif idim == 4:
+        var = var[:,:,:,1:-1]
+    elif idim == 5:
+        var = var[:,:,:,:,1:-1]
+    elif idim == 6:
+        var = var[:,:,:,:,:,1:-1]
+    elif idim == 7:
+        var = var[:,:,:,:,:,:,1:-1]
+    return var
\ No newline at end of file