diff --git a/MY_RUN/KTEST/003_KW78/004_python/plot_003_KW78.py b/MY_RUN/KTEST/003_KW78/004_python/plot_003_KW78.py
index 9ec2501f9495c586925e8dc649351a883f6db275..139b4a13f35d2d4fdd92d5e94b46716cf7f170ca 100644
--- a/MY_RUN/KTEST/003_KW78/004_python/plot_003_KW78.py
+++ b/MY_RUN/KTEST/003_KW78/004_python/plot_003_KW78.py
@@ -69,7 +69,7 @@ Ltitle = ['Wind at K=2', 'Wind at 3000m', 'Wind at 5000m']
 Lxlab = ['x (m)']*len(Lplot)
 Lylab = ['y (m)']*len(Lplot)
 Llegendval = [10,10,10]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.002]*len(Lplot)
 Lcolor = ['black']*len(Lplot)
@@ -79,7 +79,7 @@ lat = [Dvar['f1']['nj_u'], Dvar['f2']['nj'],  Dvar['f2']['nj'] ]
 Lscale = [200]*len(Lplot)
 fig2 = Panel1.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Lvar2=Lplot2, Lcarte=[500,23500,500,23500], Llevel=Llvl, 
                       Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
-                      Lcolor=Lcolor, Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[4,6,8], ax=fig1.axes, Lscale=Lscale)
+                      Lcolor=Lcolor, Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[4,6,8], ax=fig1.axes, Lscale=Lscale)
 #  Oblique projection
 i_beg, j_beg = (3,0)
 i_end, j_end = (22,21)
@@ -133,13 +133,13 @@ Lplot1 = [ WIND_sec1]
 Lplot2 = [ WT_sec1]
 Ltitle = ['Wind']
 Llegendval = [25]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.004]*len(Lplot)
 Lscale = [200]*len(Lplot)
 
 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=[0], ax=fig3.axes, Lscale=Lscale)
+                        Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[0], ax=fig3.axes, Lscale=Lscale)
 
 Lplot = [RRT_sec1]
 LaxeX = [axe_m1]
@@ -195,13 +195,13 @@ Lplot1 = [ Dvar['f1']['VM'][:,:,13]]
 Lplot2 = [ Dvar['f1']['WM'][:,:,13]]
 Ltitle = ['Wind']
 Llegendval = [25]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.004]*len(Lplot)
 Lscale = [200]*len(Lplot)
 
 fig7 = 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=fig6.axes, Lscale=Lscale)
+                        Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[0], ax=fig6.axes, Lscale=Lscale)
 
 
 Lplot = [Dvar['f1']['RRT'][:,:,13]]
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 ea59e343330e5dfc43a383bb976da141751a6705..5de8f4b897fab5ecd28427472f553f42100b6f99 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
@@ -23,7 +23,8 @@ os.system('rm -f tempgraph*')
 LnameFiles = ['REUNI.1.00A20.004dia.nc', 'REUNI.1.00A20.004.nc']
 
 Dvar_input = {
-'f1':['ZS', 'UT', 'VT', 'WT', 'THT', 'ALT_PRESSURE','ALT_U','ALT_V','ALT_THETA','level','ZTOP', 'longitude','latitude','level_w','time'],
+'f1':['ZS', 'UT', 'VT', 'WT', 'THT', 'ALT_PRESSURE','ALT_U','ALT_V','ALT_THETA','level','ZTOP', 'longitude','latitude','level_w','time',
+    'RN', 'H','LE','GFLUX','HU2M','T2M','W10M','CD','CH','Z0','CE','TS','Z0H'],
 'f2':['LSTHM', 'LSVM']}
 
 #  Read the variables in the files
@@ -65,7 +66,7 @@ Ltitle = ['wind vectors at K=2', 'wind vectors at z = 1500m ']
 Lxlab = ['longitude']*len(Lplot1)
 Lylab = ['latitude']*len(Lplot1)
 Llegendval = [25,25]
-Lcbarlabel = ['(m/s)']*len(Lplot1)
+Llegendlabel = ['(m/s)']*len(Lplot1)
 Larrowstep = [4]*len(Lplot1)
 Lwidth = [0.003]*len(Lplot1)
 Lcolor = ['black']*len(Lplot1)
@@ -73,7 +74,7 @@ Lprojection = [ccrs.PlateCarree()]*len(Lplot1)
 Llvl = [0]*len(Lplot1)
 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)
+                      Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lproj=Lprojection, Lid_overlap=[2,6], ax=fig1.axes, Lscale=Lscale)
 
 ################################################################
 #########          PANEL 2  # Vertical cross-section
@@ -122,7 +123,7 @@ Lplot1 = [ Dvar['f1']['VM'][:,:,i_slice]]
 Lplot2 = [ Dvar['f1']['WM'][:,:,i_slice]]
 Ltitle = ['Wind']
 Llegendval = [15]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Lxlab = ['longitude']*len(Lplot)
 Lylab = ['altitude (m)']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
@@ -133,6 +134,6 @@ Lxlim = [(-21.3,-20.9)]*len(Lplot)
 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)
+                        Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[6], ax=fig3.axes, Lscale=Lscale, Lylim=Lylim, Lxlim=Lxlim, Lcolor=Lcolor)
 
-Panel2.save_graph(2,fig4)
\ No newline at end of file
+Panel2.save_graph(2,fig4)
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 03bec46158136aacebbe3bb40e9d216a9c16aa85..f5292605c7d9fe3c358ad25b49fc999224725bf9 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
@@ -58,7 +58,7 @@ Ltitle = ['Wind at 850hPa', 'Wind at 700hPa', 'Wind at 9000m']
 Lxlab = ['longitude']*len(Lplot1)
 Lylab = ['latitude']*len(Lplot1)
 Llegendval = [20,20,40]
-Lcbarlabel = ['(m/s)']*len(Lplot1)
+Llegendlabel = ['(m/s)']*len(Lplot1)
 Larrowstep = [2]*len(Lplot1)
 Lwidth = [0.002]*len(Lplot1)
 Lcolor = ['black']*len(Lplot1)
@@ -66,7 +66,7 @@ Lprojection = [ccrs.PlateCarree()]*len(Lplot1)
 Llvl = [0]*len(Lplot1)
 fig2 = Panel1.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Lvar2=Lplot2, Lcarte=[], Llevel=Llvl, Lxlab=Lxlab, Lylab=Lylab, 
                       Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, Lproj=Lprojection,
-                      Lcolor=Lcolor, Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[2,4,6], ax=fig1.axes)
+                      Lcolor=Lcolor, Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[2,4,6], ax=fig1.axes)
 
 Panel1.save_graph(1,fig2)
 
@@ -99,7 +99,7 @@ Ltitle = ['Wind at 850hPa', 'Wind at 700hPa', 'Wind at 9000m']
 Llegendval = [20,20,40]
 Lxlab = ['longitude']*len(Lplot1)
 Lylab = ['latitude']*len(Lplot1)
-Lcbarlabel = ['(m/s)']*len(Lplot1)
+Llegendlabel = ['(m/s)']*len(Lplot1)
 Larrowstep = [2]*len(Lplot1)
 Lwidth = [0.002]*len(Lplot1)
 Lcolor = ['black']*len(Lplot1)
@@ -107,6 +107,6 @@ Lprojection = [ccrs.PlateCarree()]*len(Lplot1)
 Llvl = [0]*len(Lplot1)
 fig2 = Panel2.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Lvar2=Lplot2, Lcarte=[], Llevel=Llvl, Lxlab=Lxlab, Lylab=Lylab, 
                       Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, Lproj=Lprojection,
-                      Lcolor=Lcolor, Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[2,4,6], ax=fig1.axes)
+                      Lcolor=Lcolor, Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[2,4,6], ax=fig1.axes)
 
 Panel2.save_graph(2,fig2)
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 3f792a31363a78a8d5db2c9e5163afa1597afc82..9666ee0c3614aea18c5a55add922bbde82015df1 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
@@ -68,7 +68,7 @@ Ltitle = ['Wind at K=2', 'Wind at 3000m', 'Wind at 5000m']
 Lxlab = ['x (m)']*len(Lplot)
 Lylab = ['y (m)']*len(Lplot)
 Llegendval = [10,10,10]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.002]*len(Lplot)
 Lcolor = ['black']*len(Lplot)
@@ -78,7 +78,7 @@ lat = [Dvar['f1']['nj_u'], Dvar['f2']['nj'],  Dvar['f2']['nj'] ]
 Lscale = [200]*len(Lplot)
 fig2 = Panel1.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Lvar2=Lplot2, Lcarte=[500,23500,500,23500], Llevel=Llvl, 
                       Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
-                      Lcolor=Lcolor, Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[4,6,8], ax=fig1.axes, Lscale=Lscale)
+                      Lcolor=Lcolor, Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[4,6,8], ax=fig1.axes, Lscale=Lscale)
 #  Oblique projection
 i_beg, j_beg = (2,0)
 i_end, j_end = (21,22)
@@ -163,13 +163,13 @@ Lplot1 = [ WIND_sec1]
 Lplot2 = [ WT_sec1]
 Ltitle = ['Wind']
 Llegendval = [25]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.004]*len(Lplot)
 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)
+                        Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[0], ax=fig3.axes, Lscale=Lscale)
 
 Lplot = [RRT_sec1]
 LaxeX = [axe_m1]
@@ -287,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)
\ No newline at end of file
+Panel5.save_graph(5,fig8)
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 2116f12329b9163e64761721f006b0e9e06efd1f..34819dabe4071852324a264955e5668dc147fe89 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
@@ -64,14 +64,14 @@ Lplot1 = [ Dvar['f1']['UM']]
 Lplot2 = [ Dvar['f1']['VM']]
 Ltitle = [' vectors at 1st level']
 Llegendval = [20]
-Lcbarlabel = ['m/s']*len(Lplot)
+Llegendlabel = ['m/s']*len(Lplot)
 Larrowstep = [1]*len(Lplot)
 Lwidth = [0.002]*len(Lplot)
 Lcolor = ['black']*len(Lplot)
 Lscale = [200]*len(Lplot)
 
 fig2 = Panel1.pvector(Lxx=lon, Lyy=lat, Lvar1=Lplot1, Llevel=Llvl, Lvar2=Lplot2, Lxlab=Lxlab, Lylab=Lylab, Ltitle=Ltitle, Lwidth=Lwidth, Larrowstep=Larrowstep, 
-                        Llegendval=Llegendval, Lcbarlabel=Lcbarlabel, Lid_overlap=[0], ax=fig1.axes, Lscale=Lscale)
+                        Llegendval=Llegendval, Llegendlabel=Llegendlabel, Lid_overlap=[0], ax=fig1.axes, Lscale=Lscale)
 
 Panel1.save_graph(1,fig2)
 
@@ -121,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)
\ No newline at end of file
+Panel3.save_graph(3,fig4)
diff --git a/src/LIB/Python/Panel_Plot.py b/src/LIB/Python/Panel_Plot.py
index 328d26264ad1f209de65a0ecdc36f563e07dde70..2072d1cf9225c27f14b79849e8b11e5e49885ce8 100644
--- a/src/LIB/Python/Panel_Plot.py
+++ b/src/LIB/Python/Panel_Plot.py
@@ -9,7 +9,7 @@ MNH_LIC for details. version 1.
 @author: 07/2021 Quentin Rodier
 """
 import matplotlib as mpl
-mpl.use('Agg')
+#mpl.use('Agg')
 import matplotlib.pyplot as plt
 from matplotlib import cm
 from matplotlib.colors import ListedColormap
@@ -21,7 +21,7 @@ class PanelPlot():
     
     def __init__(self, nb_l, nb_c, Lfigsize, bigtitle, titlepad=40, minmaxpad=1.03, timepad=-0.06, lateralminmaxpad=0.86, 
                  labelcolorbarpad=6.0, colorbaraspect=20, colorbarpad=0.04, tickspad=0.8,
-                 minmaxTextSize=10, bigtitleSize=13, titleSize=12,
+                 minmaxTextSize=10, bigtitleSize=13, titleSize=12, legendSize=10,
                  xlabelSize=11, ylabelSize=11, timeSize=11, cbTicksLabelSize=11, cbTitleSize=11, xyTicksLabelSize=10, figBoxLinewidth=1,
                  xyTicksWidth=1, xyTicksLength=6):
 
@@ -45,6 +45,7 @@ class PanelPlot():
         self.titleSize = titleSize               #  Graph title fontsize
         self.xlabelSize = xlabelSize             #  X-label fontsize
         self.ylabelSize = ylabelSize             #  Y-label fontsize
+        self.legendSize = legendSize             #  X/Y plot legend fontsize
         self.timeSize = timeSize                 #  Time attribute of the graphs fontsize
         self.cbTicksLabelSize = cbTicksLabelSize #  Colorbar ticks label fontsize
         self.cbTitleSize = cbTitleSize           #  Colorbar title fontsize
@@ -57,7 +58,7 @@ class PanelPlot():
         #  Initialization of the panel plots
         self.fig = plt.figure(figsize=(self.Lfigsize[0],self.Lfigsize[1]))
         self.fig.set_dpi(125)
-        self.fig.suptitle(self.bigtitle,fontsize=16)
+        self.fig.suptitle(self.bigtitle,fontsize=bigtitleSize)
 
     def save_graph(self, iplt, fig, fig_name='tempgraph'):
       """
@@ -221,7 +222,7 @@ class PanelPlot():
                  Lstep=[], Lstepticks=[], Lcolormap=[], Lcbarlabel=[], LcolorLine=[],
                  Lfacconv=[], ax=[], Lid_overlap=[], colorbar=True, orog=[], Lxlim=[], Lylim=[], Ltime=[], Lpltype=[], LaddWhite_cm=[], LwhiteTop=[]):
       """
-        Horizontal cross section plot
+        Vertical cross section plot
         Parameters :
             - Lxx    : List of x or y coordinate variable or time axis
             - Lzz    : List of z coordinates variable
@@ -246,6 +247,7 @@ class PanelPlot():
             - colorbar   : show colorbar or not
             - LaddWhite_cm : List of boolean to add white color to a colormap at the last bottom (low value) tick colorbar
             - LwhiteTop    : List of boolean to add the white color at the first top (high value). If false, the white is added at the bottom if Laddwhite_cm=T
+            - orog         : Orography variable
       """
       self.ax = ax
       firstCall = (len(self.ax) == 0)
@@ -412,13 +414,16 @@ class PanelPlot():
         #  Legend
         #TODO : Handling legend with overlap two axis lines in the same box. For now, placement is by hand
         if not id_overlap: 
-          self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.95))
+          self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.95),fontsize=self.legendSize)
         else:
-          self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.90))
+          self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.90),fontsize=self.legendSize)
 
         #  Title
         if Ltitle: self.set_Title(self.ax, iax, Ltitle[i], id_overlap,Lxlab[i], Lylab[i])
     
+        #  Ticks label 
+        self.ax[iax].tick_params(axis='both', labelsize=self.xyTicksLabelSize, width=self.xyTicksWidth, length=self.xyTicksLength, pad=self.tickspad)
+        
         #  X/Y Axis label
         if id_overlap: 
           self.ax[iax].xaxis.tick_top()
@@ -577,7 +582,9 @@ class PanelPlot():
     
         #  X/Y Axis
         self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
-        
+        #  Ticks label 
+        self.ax[iax].tick_params(axis='both', labelsize=self.xyTicksLabelSize, width=self.xyTicksWidth, length=self.xyTicksLength, pad=self.tickspad)
+       
         #  Color label on contour-line
         if Lpltype[i]=='c': #  Contour
           if 'GeoAxes' in str(type(self.ax[self.i])): # cartopy does not like the levels arguments in clabel, known issue
@@ -594,11 +601,11 @@ class PanelPlot():
       return self.fig
     
     def pvector(self, Lxx=[], Lyy=[], Lvar1=[], Lvar2=[], Lcarte=[], Llevel=[], Lxlab=[], Lylab=[], 
-                Ltitle=[], Lwidth=[], Larrowstep=[], Lcolor=[], Llegendval=[], Lcbarlabel=[], 
+                Ltitle=[], Lwidth=[], Larrowstep=[], Lcolor=[], Llegendval=[], Llegendlabel=[], 
                 Lproj=[], Lfacconv=[], ax=[], coastLines=True, Lid_overlap=[], Ltime=[], Lscale=[],
                 Lylim=[], Lxlim=[]):
       """
-        Horizontal vectors lines
+        Vectors
         Parameters :
             - Lxx    : List of x or y coordinate variable (lat or ni or nm)
             - Lyy    : List of y coordinates variable (lon or level)
@@ -617,7 +624,7 @@ class PanelPlot():
             - Larrowstep : List of sub-sample (frequency) if too much arrows
             - Lcolor : List of colors for the arrows (default: black)
             - Llegendval : List of value for the legend of the default arrow
-            - Lcbarlabel : List of labels for the legend of the default arrow
+            - Llegendlabel : List of labels for the legend of the default arrow
             - Lproj      : List of ccrs cartopy projection
             - Lfacconv   : List of factors for unit conversion of each variables
             - coastLines : Boolean to plot coast lines and grid lines
@@ -703,7 +710,7 @@ class PanelPlot():
         if Lproj: self.draw_Backmap(coastLines, self.ax[iax], Lproj[i])
     
         # Arrow legend key
-        qk = self.ax[iax].quiverkey(cf, 1.0, -0.05, Llegendval[i], str(Llegendval[i]) + Lcbarlabel[i], labelpos='E', color='black')
+        qk = self.ax[iax].quiverkey(cf, 1.0, -0.05, Llegendval[i], str(Llegendval[i]) + Llegendlabel[i], labelpos='E', color='black')
            
       return self.fig
 
diff --git a/src/LIB/Python/read_MNHfile.py b/src/LIB/Python/read_MNHfile.py
index fba716062817184b910fdab3d1bcf2d7bd414082..38072becffe3b1d8c066c6f370683b4988cf7ccf 100644
--- a/src/LIB/Python/read_MNHfile.py
+++ b/src/LIB/Python/read_MNHfile.py
@@ -104,10 +104,10 @@ def read_netcdf(LnameFiles, Dvar_input, path='.', get_data_only=True, del_empty_
                 Dvar[keyFiles] = read_TIMESfiles_55(theFile, Dvar_input[keyFiles], Dvar[keyFiles], get_data_only, del_empty_dim, removeHALO)
         else:
             Dvar[keyFiles]= read_BACKUPfile(theFile, Dvar_input[keyFiles], Dvar[keyFiles], get_data_only, del_empty_dim, removeHALO)
-        theFile.close()
+        #theFile.close()
     return Dvar
 
-def read_var(theFile, Dvar, var_name, get_data_only=True, del_empty_dim=True, removeHALO=True):
+def read_var(theFile, Dvar, var_name, get_data_only, del_empty_dim=True, removeHALO=True):
     """Read a netCDF4 variable
     
     Parameters
@@ -141,7 +141,6 @@ def read_var(theFile, Dvar, var_name, get_data_only=True, del_empty_dim=True, re
         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: 
@@ -180,11 +179,10 @@ def read_var(theFile, Dvar, var_name, get_data_only=True, del_empty_dim=True, re
                 except IndexError:
                     break
             Ldimtosqueeze=tuple(Ldimtosqueeze)
-            Dvar[var_name] = np.squeeze(Dvar[var_name], axis=Ldimtosqueeze) 
-        
+            Dvar[var_name] = np.squeeze(Dvar[var_name], axis=Ldimtosqueeze)
     return Dvar
 
-def read_from_group(theFile, Dvar, group_name, var_name, get_data_only=True, del_empty_dim=True,removeHALO=True):
+def read_from_group(theFile, Dvar, group_name, var_name, get_data_only, del_empty_dim=True,removeHALO=True):
     """Read a variable from a netCDF4 group 
     
     Parameters
@@ -267,7 +265,7 @@ def read_from_group(theFile, Dvar, group_name, var_name, get_data_only=True, del
             Dvar[(group_name,var_name)] = Dvar[(group_name,var_name)].T   
     return Dvar
 
-def read_BACKUPfile(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_dim=True, removeHALO=True):
+def read_BACKUPfile(theFile, Dvar_input, Dvar, get_data_only, 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
     
@@ -301,19 +299,18 @@ def read_BACKUPfile(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_dim
     #  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)
+            if get_data_only:
+                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, Dvar_input, Dvar, get_data_only=True, del_empty_dim=True, removeHALO=True):
+def read_TIMESfiles_55(theFile, Dvar_input, Dvar, get_data_only, 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
 
@@ -345,7 +342,6 @@ def read_TIMESfiles_55(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_
     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 = read_from_group(theFile, Dvar, var[0], var[1], get_data_only, del_empty_dim, removeHALO)
         else:
@@ -354,6 +350,8 @@ def read_TIMESfiles_55(theFile, Dvar_input, Dvar, get_data_only=True, del_empty_
 
 def removetheHALO(idim, var):
     """Remove a NHALO=1 point [1:-1] at a given dimension idim of a variable var
+    TODO: with NHALO=3 with WENO5
+    TODO: with the use of get_data_only=False (NetCDF4 output)
     
     Parameters
     ----------
diff --git a/src/LIB/Python/readme b/src/LIB/Python/readme
new file mode 100644
index 0000000000000000000000000000000000000000..ddb601ec1ec0ed8b757ed310dda73811edb0840d
--- /dev/null
+++ b/src/LIB/Python/readme
@@ -0,0 +1 @@
+Tutorial given the 17/03/2022 to use the library can be found at https://github.com/QuentinRodier/MNHPy/blob/main/examples/mnhpy.ipynb
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 20c9a078db5e550dbf3bb03cd2fcc6b13ddc89a5..b2df24bd9b4d00d0e37550e347af3dcef9fc3abf 100644
--- a/src/MNH/compute_bl89_ml.f90
+++ b/src/MNH/compute_bl89_ml.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2006-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2006-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.
@@ -57,6 +57,7 @@ END MODULE MODI_COMPUTE_BL89_ML
 !!     R.Honnert Oct 2016 : Update with AROME
 !!     Q.Rodier  01/2019 : support RM17 mixing length as in bl89.f90 
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!  P. Wautelet 17/12/2021: bugfix: KK instead of JKK
 !!
 !-------------------------------------------------------------------------------
 !
@@ -162,9 +163,9 @@ IF (OUPORDN.EQV..TRUE.) THEN
      ! Lenght travelled by parcel to nullify energy
       ZLWORK2(J1D)=        ( - PG_O_THVREF(J1D) *                     &
             (  ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) )                              &
-            - XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+            - XRM17*PSHEAR(J1D,KK)*sqrt(abs(PTKEM_DEP(J1D))) &
           + SQRT (ABS(                                                       &
-            (XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) +  &
+            (XRM17*PSHEAR(J1D,KK)*sqrt(abs(PTKEM_DEP(J1D))) +  &
              PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2  &
             + 2. * ZINTE(J1D) * PG_O_THVREF(J1D)                        &
                  * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) ))    ) /             &
diff --git a/src/MNH/compute_entr_detr.f90 b/src/MNH/compute_entr_detr.f90
index 80a9d68db57ef58f31574c47e144887ab441ce7a..004f706a183ee10d37a85b4fcd79fa78d801209f 100644
--- a/src/MNH/compute_entr_detr.f90
+++ b/src/MNH/compute_entr_detr.f90
@@ -122,6 +122,7 @@ END MODULE MODI_COMPUTE_ENTR_DETR
 !!      R.Honnert Oct 2016 : Update with AROME
 !  P. Wautelet 08/02/2019: bugfix: compute ZEPSI_CLOUD only once and only when it is needed
 !  P. Wautelet 10/02/2021: bugfix: initialized PPART_DRY everywhere
+!  M. Mandement 24/01/2022:bugfix: init of theta_l in mixtures was too low (0.1K) 
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -232,7 +233,7 @@ INTEGER :: JI,JLOOP
   ZFRAC_ICE(:)=PFRAC_ICE(:) ! to not modify fraction of ice
  
   ZPRE(:)=PPRE_MINUS_HALF(:)
-  ZMIXTHL(:)=0.1
+  ZMIXTHL(:)=300.0
   ZMIXRT(:)=0.1
 
   !Initialize PPART_DRY everywhere to prevent access to non-initialized values
diff --git a/src/MNH/isocom.f b/src/MNH/isocom.f
index 451c1b485c342af59e48c962566e12ed682dc631..42481f83a3f489686c70fa8a696f9ffe8966ff04 100644
--- a/src/MNH/isocom.f
+++ b/src/MNH/isocom.f
@@ -3,6 +3,13 @@ CMNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 CMNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 CMNH_LIC for details. version 1.
 C=======================================================================
+C Modifications:
+C  P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+C  P. Wautelet 22/01/2019: replace obsolete SNGL intrinsics by REAL intrinsics
+C  P. Wautelet 19/04/2019: use kind(0.0d0) instead of kind=8
+C  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+C  P. Wautelet 17/12/2021: add missing definitions of parameter ONE (in POLY3 and POLY3B)
+C=======================================================================
 C
 C *** ISORROPIA CODE
 C *** SUBROUTINE ISOROPIA
@@ -123,11 +130,6 @@ C
 C *** COPYRIGHT 1996-2000, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY
 C *** WRITTEN BY ATHANASIOS NENES
 C
-C Modifications:
-C  P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
-C  P. Wautelet 22/01/2019: replace obsolete SNGL intrinsics by REAL intrinsics
-C  P. Wautelet 19/04/2019: use kind(0.0d0) instead of kind=8
-C  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 C=======================================================================
 C
       SUBROUTINE ISOROPIA (WI, RHI, TEMPI,  CNTRL, 
@@ -3882,7 +3884,7 @@ C
       SUBROUTINE POLY3B (A1, A2, A3, RTLW, RTHI, ROOT, ISLV)
 C
       IMPLICIT REAL(kind(0.0d0))           (A-H, O-Z)
-      PARAMETER (ZERO=0.D0, EPS=1D-15, MAXIT=100, NDIV=5)
+      PARAMETER (ZERO=0.D0, ONE=1.D0, EPS=1D-15, MAXIT=100, NDIV=5)
 C
       FUNC(X) = X**3.d0 + A1*X**2.0 + A2*X + A3
 C
diff --git a/src/MNH/mode_RBK90_Integrator.f90 b/src/MNH/mode_RBK90_Integrator.f90
index ae11720e3cd8fdc1ce7ab60d602dd20a5e4ce0d5..536726def5d7ba4cca6db9fab1f67c5e867db171 100644
--- a/src/MNH/mode_RBK90_Integrator.f90
+++ b/src/MNH/mode_RBK90_Integrator.f90
@@ -1,10 +1,11 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 ! Modifications:
-!  Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+!  P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+!  P. Wautelet 17/12/2021: remove unused time variables (and remove use of non initialised time values)
 !-----------------------------------------------------------------
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! 
@@ -1316,9 +1317,6 @@ USE MODI_CH_FCN
     REAL :: T, Y(N)
 !~~~> Output variables
     REAL :: Ydot(N)
-!~~~> Local variables
-    REAL :: Told
-    REAL :: TIME
 !
 INTEGER, INTENT(IN) :: KMI      ! model number
 INTEGER, INTENT(IN) :: KVECNPT
@@ -1330,8 +1328,6 @@ INTEGER, INTENT(IN) :: KVECNPT
     INTEGER       :: ISPECIES      ! Ancillary variable
 
 
-    Told = TIME
-    TIME = T
 !JPP   CALL Update_SUN()
 !JPP   CALL Update_RCONST()
 !JPP   CALL Fun( Y, FIX, RCONST, Ydot )
@@ -1352,7 +1348,6 @@ INTEGER, INTENT(IN) :: KVECNPT
       Ydot(JLSHIFT+1:JLSHIFT+ISPECIES) = ZFCN(JL,1:ISPECIES)
       JLSHIFT = JLSHIFT + ISPECIES
     END DO
-    TIME = Told
 
 END SUBROUTINE FunTemplate
 
@@ -1389,8 +1384,6 @@ INTEGER, INTENT(IN) :: KMI      ! model number
 INTEGER, INTENT(IN) :: KVECNPT
 !
 !~~~> Local variables
-    REAL :: Told
-    REAL :: TIME
 !#ifdef FULL_ALGEBRA    
 !##    INTEGER :: i, j
 !#endif   
@@ -1399,8 +1392,6 @@ INTEGER, INTENT(IN) :: KVECNPT
     INTEGER       :: JL, JLL, JLSHIFT   ! Loop indexes
     INTEGER       :: ISPECIES      ! Ancillary variable
 
-    Told = TIME
-    TIME = T
 !JPP    CALL Update_SUN()
 !JPP    CALL Update_RCONST()
 !#ifdef FULL_ALGEBRA    
@@ -1441,7 +1432,6 @@ INTEGER, INTENT(IN) :: KVECNPT
         JLSHIFT = JLSHIFT + NSPARSEDIM_NAQ
       END IF
     END DO
-    TIME = Told
 
 END SUBROUTINE JacTemplate
 
diff --git a/src/MNH/mode_RBK90_linearalgebra.f90 b/src/MNH/mode_RBK90_linearalgebra.f90
index 8fd23aa139e89b8c906915ccda98c36e80e102e0..de787d50d3e04f7b43d59230b4f92f7d525bf730 100644
--- a/src/MNH/mode_RBK90_linearalgebra.f90
+++ b/src/MNH/mode_RBK90_linearalgebra.f90
@@ -1,10 +1,11 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 ! Modifications:
 !  P. Wautelet 22/02/2019: DOUBLE COMPLEX -> COMPLEX(kind(0.0d0)) to respect Fortran standard
+!  P. Wautelet 17/12/2021: convert arithmetic if to classic if (deleted from Fortran 2018 standard)
 !-----------------------------------------------------------------
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! 
@@ -1261,46 +1262,49 @@ END SUBROUTINE KppSolveTR
                                  
       WDOT = 0.0D0 
       IF (N .LE. 0) RETURN 
-      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
+      IF (incX .EQ. incY) THEN
 !                                                                       
 !     Code for unequal or nonpositive increments.                       
 !                                                                       
-    5 IX = 1 
-      IY = 1 
-      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
-      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
-      DO i = 1,N 
-        WDOT = WDOT + DX(IX)*DY(IY) 
-        IX = IX + incX 
-        IY = IY + incY 
-      END DO 
-      RETURN 
+        IF ( incX < 1 ) THEN
+          IX = 1
+          IY = 1
+          IF (incX .LT. 0) IX = (-N+1)*incX + 1
+          IF (incY .LT. 0) IY = (-N+1)*incY + 1
+          DO i = 1,N
+            WDOT = WDOT + DX(IX)*DY(IY)
+            IX = IX + incX
+            IY = IY + incY
+          END DO
+        ELSE IF ( incX == 1 ) THEN
 !                                                                       
 !     Code for both increments equal to 1.                              
 !                                                                       
 !     Clean-up loop so remaining vector length is a multiple of 5.      
 !                                                                       
-   20 M = MOD(N,5) 
-      IF (M .EQ. 0) GO TO 40 
-      DO i = 1,M 
-         WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
-      IF (N .LT. 5) RETURN 
-   40 MP1 = M + 1 
-      DO i = MP1,N,5 
-          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
-                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
-      END DO 
-      RETURN 
+          M = MOD(N,5)
+          IF (M .EQ. 0) GO TO 40
+          DO i = 1,M
+            WDOT = WDOT + DX(i)*DY(i)
+          END DO
+          IF (N .LT. 5) RETURN
+   40     MP1 = M + 1
+          DO i = MP1,N,5
+            WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
+                     DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
+          END DO
+        ELSE
 !                                                                       
 !     Code for equal, positive, non-unit increments.                    
 !                                                                       
-   60 NS = N*incX 
-      DO i = 1,NS,incX 
-        WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
+          NS = N*incX
+          DO i = 1,NS,incX
+            WDOT = WDOT + DX(i)*DY(i)
+          END DO
+        END IF
+      END IF
 
-      END FUNCTION WDOT                                          
+      END FUNCTION WDOT
 
 
 !--------------------------------------------------------------
diff --git a/src/MNH/nband_model.fx90 b/src/MNH/nband_model.fx90
index e0aface0e240a114c19048fbf5ddf99cb6c969f5..6005b838759cf188856ece7a479f72945cf033a2 100644
--- a/src/MNH/nband_model.fx90
+++ b/src/MNH/nband_model.fx90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1996-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1996-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.
@@ -13,6 +13,7 @@
 *                                named nband_model.f90 and compiled with -Fixed
 *           J.Escobar (1/12/2017) bug => intialized all ZV=0.0 in spectr
 *           P. Wautelet 21/11/2019: replace several CONTINUE (workaround of problems with gfortran OpenACC)
+*           P. Wautelet 17/12/2021: comment ZBSUI variable (not used and was not initialized)
 *
       SUBROUTINE NBMVEC
      I ( KIDIA  ,KFDIA ,KLON,KLEV,KGL,KCABS,KNG1,KUABS
@@ -768,7 +769,7 @@ C
       DO 255 JL=KIDIA,KFDIA
       ZBLEV(JL,KLEV+1)=ZRES(JL)
       ZBINT(JL,KLEV+1)=ZBINT(JL,KLEV+1)+ZRES(JL)
-      ZBSUI(JL)=ZBSUI(JL)+ZBSUR(JL)
+C      ZBSUI(JL)=ZBSUI(JL)+ZBSUR(JL)
  255  CONTINUE
 C      IF (NIMP.EQ.0) THEN
 C      JL=KIDIA
diff --git a/src/MNH/nhoa_coupln.f90 b/src/MNH/nhoa_coupln.f90
index f8ea30639b021065de23652e151303524f0644f4..072dfbeb3ca7aa64c133a8d0b2855e4a3d1fdea8 100644
--- a/src/MNH/nhoa_coupln.f90
+++ b/src/MNH/nhoa_coupln.f90
@@ -49,6 +49,7 @@ SUBROUTINE NHOA_COUPL_n(KDAD,PTSTEP,KMI,KTCOUNT,KKU)
 !!   JL Redelsperger 03/2021   Version 0 
 !!    MODIFICATIONS
 !!    -------------
+!  P. Wautelet 17/12/2021: bugfix: set IIU and IJU values
 !!-----------------------------------------------------------------------------
 !
 !*      0.   DECLARATIONS
@@ -98,17 +99,20 @@ CHARACTER(LEN=4)                    :: ZINIT_TYPE
 !---Coupled OA MesoNH----------------------------------------------------------------------------
 !*      0.   INITIALISATION
 !            --------------
+IIU = SIZE(XRHODJ,1)
+IJU = SIZE(XRHODJ,2)
+
 ! allocate flux local array
-ALLOCATE(ZCOUPTFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPUFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPTFL(IIU,IJU))
+ALLOCATE(ZCOUPUFL(IIU,IJU))
+ALLOCATE(ZCOUPVFL(IIU,IJU))
 ! allocate sfc variable local array
-ALLOCATE(ZCOUPUA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPTA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPUO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPTO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPUA(IIU,IJU))
+ALLOCATE(ZCOUPVA(IIU,IJU))
+ALLOCATE(ZCOUPTA(IIU,IJU))
+ALLOCATE(ZCOUPUO(IIU,IJU))
+ALLOCATE(ZCOUPVO(IIU,IJU))
+ALLOCATE(ZCOUPTO(IIU,IJU))
 ! values in ocean sfc
 IKE=KKU-JPVEXT
 ZCOUPUO(:,:)= XUT(:,:,IKE)
diff --git a/src/MNH/recycling.f90 b/src/MNH/recycling.f90
index 9734eebc00e2f5bc16e8e923bfdc04cba2f120e7..18e0956cd05cb710e26718962a054e801d79636b 100644
--- a/src/MNH/recycling.f90
+++ b/src/MNH/recycling.f90
@@ -65,35 +65,17 @@ END MODULE MODI_RECYCLING
 !**** 0. DECLARATIONS
 !     ---------------
 !
-! module
-USE MODE_POS
-USE MODE_ll
-USE MODE_IO
-!USE MODI_SHUMAN
-!
-USE MODD_PARAMETERS
-USE MODD_CONF
-!
-USE MODD_CST
-!
-USE MODD_DIM_n
-USE MODD_CONF
-USE MODD_CONF_n
-USE MODD_GRID
-USE MODD_GRID_n
-USE MODD_METRICS_n
-USE MODD_TIME
-USE MODD_TIME_n
-USE MODD_DYN_n
-USE MODD_FIELD_n
-USE MODD_CURVCOR_n
-USE MODD_REF
-!
-USE MODD_VAR_ll,          ONLY: IP, NPROC
+USE MODD_CONF,           ONLY: CCONF
+USE MODD_FIELD_n,        ONLY: XTHT, XUT, XVT, XWT
+USE MODD_GRID_n,         ONLY: XXHAT, XYHAT
+USE MODD_METRICS_n,      ONLY: XDZZ
+USE MODD_PARAMETERS,     ONLY: JPHEXT
 USE MODD_RECYCL_PARAM_n
+
+USE MODE_MSG
+
 USE MODI_RECYCL_FLUC
-USE MODD_LUNIT_n,     ONLY : TLUOUT
-!
+
 IMPLICIT NONE
 !
 !------------------------------------------------------------------------------
@@ -106,9 +88,6 @@ REAL, DIMENSION(:,:)     ,INTENT(INOUT) :: PFLUCTUNE,PFLUCTVTE,PFLUCTVNS,PFLUCTU
 !------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
-INTEGER                                      :: IIU,IJU,IKU,JCOUNT,ICOUNT,ILUOUT
-INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE,IIP
-INTEGER :: IIBG,IIEG,IJBG,IJEG,IIMAX,IJMAX
 INTEGER :: PMINW,PMINE,PMINN,PMINS
 INTEGER :: JIDIST,JJDIST
 REAL    :: Z_DELTX,Z_DELTY
@@ -116,24 +95,14 @@ REAL    :: Z_DELTX,Z_DELTY
 !------------------------------------------------------------------------------
 !
 !       0.3  allocation
-CALL GET_DIM_EXT_ll('B',IIU,IJU)
-IKU=NKMAX+2*JPVEXT
 PMINW=0
 PMINN=0
 PMINS=0
 PMINE=0
 
-CALL GET_OR_ll('B',IIBG,IJBG)
-IIBG = IIBG+IIB-1
-IJBG = IJBG+IJB-1
-CALL GET_GLOBALDIMS_ll( IIMAX,IJMAX)
-IIEG=IIBG+IIE-IIB
-IJEG=IJBG+IJE-IJB
 Z_DELTX = XXHAT(2)-XXHAT(1)
 Z_DELTY = XYHAT(2)-XYHAT(1)
 
-
-ILUOUT = TLUOUT%NLU
 !------------------------------------------------------------------------------
 !       
 !**** 1. Recycling distance calculation 
diff --git a/src/MNH/subl_blowsnow.f90 b/src/MNH/subl_blowsnow.f90
index ad0cf278fa86b1eeedaf7c0a2da1344d8d86024a..96c6b9be4432ef3ee07b4f7d79368c232ca113d7 100644
--- a/src/MNH/subl_blowsnow.f90
+++ b/src/MNH/subl_blowsnow.f90
@@ -696,6 +696,7 @@ ZFCRI = MAX(ZFCRI1,ZFCRI2)
 !*      2    Calculate variances of the horizontal and vertical velocity components
 !
 ZS0   = ZFCRI*PZZ/PVMOD
+STOP 'Bug in TURB_FLUC: ZUSTAR used but not set'
 ZSIGU = 4.77 *ZUSTAR**2/ (1+33*ZS0)**0.66666
 ZSIGV = 2.76 *ZUSTAR**2/ (1+9.5*ZS0)**0.66666
 ZSIGW = 1.31 *ZUSTAR**2/ (1+3.12*ZS0)**0.66666
diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90
index 639cd19a0b7724e63c64d472537254f37d74b73e..0b73309297d87e82a1ce438e0fb1f7bd5cc38197 100644
--- a/src/MNH/turb.f90
+++ b/src/MNH/turb.f90
@@ -618,7 +618,7 @@ ZRVORD= XRV / XRD
 !
 !
 !Copy data into ZTHLM and ZRM only if needed
-IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. ORMC01) THEN
+IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. HTURBLEN=='ADAP' .OR. ORMC01) THEN
   ZTHLM(:,:,:) = PTHLT(:,:,:)
   ZRM(:,:,:,:) = PRT(:,:,:,:)
 END IF
diff --git a/src/MNH/write_desfmn.f90 b/src/MNH/write_desfmn.f90
index fb24c9bae47e3b7140e3a266768cffe9a5216841..8be686e6450155250e7d03c8712d67f659048a94 100644
--- a/src/MNH/write_desfmn.f90
+++ b/src/MNH/write_desfmn.f90
@@ -434,7 +434,6 @@ IF(LBU_RRH) WRITE(UNIT=ILUSEG,NML=NAM_BU_RRH)
 IF(LBU_RSV) WRITE(UNIT=ILUSEG,NML=NAM_BU_RSV)
 IF(LLES_MEAN .OR. LLES_RESOLVED .OR. LLES_SUBGRID .OR. LLES_UPDRAFT  &
 .OR. LLES_DOWNDRAFT .OR. LLES_SPECTRA) WRITE(UNIT=ILUSEG,NML=NAM_LES)
-WRITE(UNIT=ILUSEG,NML=NAM_BLANKn)
 IF(LFORCING .OR. LTRANS) WRITE(UNIT=ILUSEG,NML=NAM_FRC)
 IF(CCLOUD(1:3) == 'ICE')  WRITE(UNIT=ILUSEG,NML=NAM_PARAM_ICE)
 IF(CCLOUD == 'C2R2' .OR. CCLOUD == 'C3R5' .OR. CCLOUD == 'KHKO') &
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 08471b0cfbfcef022df6f7a06d7cd23c2cea0c31..0bb1d984b2ca0d5c8a88ff9c8f53b2255dafbb69 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -2666,9 +2666,10 @@ IF ( LMEAN_POVO ) THEN
   IWORK1(:,:)=0
   ZWORK21(:,:)=0.
   IF (XMEAN_POVO(1)>XMEAN_POVO(2)) THEN
-    XMEAN_POVO(1) = ZX0D
-    XMEAN_POVO(2) = XMEAN_POVO(1)
-    ZX0D          = XMEAN_POVO(2)
+    !Invert values (smallest must be first)
+    ZX0D = XMEAN_POVO(1)
+    XMEAN_POVO(1) = XMEAN_POVO(2)
+    XMEAN_POVO(2) = ZX0D
   END IF
   DO JK=IKB,IKE
     WHERE((XPABST(:,:,JK)>XMEAN_POVO(1)).AND.(XPABST(:,:,JK)<XMEAN_POVO(2)))
diff --git a/src/Makefile.MESONH.mk b/src/Makefile.MESONH.mk
index 72a5827a7f72fce43c57c373f3a00eadac2a000d..7e99233a2fc6af6ae1e525a5da2cf25a390f4e88 100644
--- a/src/Makefile.MESONH.mk
+++ b/src/Makefile.MESONH.mk
@@ -648,8 +648,8 @@ endif
 #
 ifeq "$(VER_CDF)" "CDFCTI"
 CDF_PATH?=/usr
-INC_NETCDF     = -I${CDF_PATH}/include
-LIB_NETCDF     = -L${CDF_PATH}/lib64 -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 -lsz -lz
+INC_NETCDF     ?= $(shell nf-config --fflags)
+LIB_NETCDF     ?= $(shell nf-config --flibs)
 INC            += $(INC_NETCDF)
 LIBS           += $(LIB_NETCDF)
 endif
diff --git a/src/SURFEX/init_megann.F90 b/src/SURFEX/init_megann.F90
index ee3b097f8a6b1a5ee6374c4458225f0cd1d16d81..cfb05aecfc8ad0e6d332681ec41d1da73c983a61 100644
--- a/src/SURFEX/init_megann.F90
+++ b/src/SURFEX/init_megann.F90
@@ -26,6 +26,7 @@ SUBROUTINE INIT_MEGAN_n(IO, S, K, NP, MSF, MGN, PLAT, HSV, PMEGAN_FIELDS)
 !!    Original: 25/10/14
 !!    Modified: 06/2017, J. Pianezze, adaptation for SurfEx v8.0
 !!    Modified: 06/2018, P. Tulet,  add PFT and LAI
+!!    Modified: 11/2021, P. Tulet,  update PFT with Ecoclimap-SG
 !!
 !!
 !!    EXTERNAL
@@ -40,10 +41,7 @@ USE MODD_MEGAN_n, ONLY : MEGAN_t
 USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
 USE MODD_ISBA_n, ONLY : ISBA_S_t, ISBA_P_t, ISBA_K_t, ISBA_NP_t
 !
-USE MODD_DATA_COVER_PAR, ONLY : NVT_C4, NVT_TRBE, NVT_TRBD, NVT_TEBE, &
-                    NVT_TEBD, NVT_TENE, NVT_BOBD, NVT_BONE, NVT_BOND, &
-                    NVT_BOGR, NVT_SHRB, NVT_GRAS, NVT_TROG, NVT_C3,   &
-                    NVT_NO, NVT_ROCK, NVT_SNOW, NVT_IRR, NVT_PARK
+USE MODD_DATA_COVER_PAR
 !
 USE MODD_SURF_PAR,   ONLY : XUNDEF
 USE MODD_DATA_COVER, ONLY : XDATA_LAI
@@ -94,6 +92,7 @@ REAL,DIMENSION(SIZE(K%XCLAY,1)) :: ZLAI
 !
 ALLOCATE(MGN%XPFT   (N_MGN_PFT,SIZE(K%XCLAY,1)))
 ALLOCATE(MGN%XEF    (N_MGN_SPC,SIZE(K%XCLAY,1)))
+ALLOCATE(MGN%XLAI   (SIZE(K%XCLAY,1)))
 ALLOCATE(MGN%NSLTYP (SIZE(K%XCLAY,1)))
 ALLOCATE(MGN%XBIOFLX(SIZE(K%XCLAY,1)))
 MGN%XBIOFLX(:) = 0.
@@ -313,134 +312,91 @@ ENDDO
 !
 ! 1 Needleleaf evergreen temperate tree
 ! -------------------------------------
-! utilisation de la classe NVT_CONI pour 30 < LAT < 60 
-WHERE ((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.))
-  MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
-END WHERE
-WHERE ((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.))
-  MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
-END WHERE
+! utilisation de la classe NVT_TENE 
+MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
 !
 ! 2 Needleleaf evergreen boreal tree
 ! -------------------------------------
-!utilisation de la classe NVT_CONI  pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
-  MGN%XPFT(2,:) = S%XVEGTYPE(:,NVT_BONE)
-END WHERE
+!utilisation de la classe NVT_BONE 
+MGN%XPFT(2,:) = S%XVEGTYPE(:,NVT_BONE)
 !
 !3 Needleleaf deciduous boreal tree
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
-  MGN%XPFT(3,:) = S%XVEGTYPE(:,NVT_BOND)
-END WHERE
+!utilisation de la classe NVT_BOND 
+MGN%XPFT(3,:) = S%XVEGTYPE(:,NVT_BOND)
 !
 !4 Broadleaf evergreen tropical tree
 ! -------------------------------------
-!utilisation de la classe NVT_EVER pour -30 < LAT < 30
-! et une hauteur d'arbre supérieur à 3 m
-WHERE (((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.)).AND.&
-       (ZH_TREE(:,IP_TRBE) .GE. 3.).AND.(ZH_TREE(:,IP_TRBE) .NE. XUNDEF))
+!utilisation de la classe NVT_TRBE
 MGN%XPFT(4,:) = S%XVEGTYPE(:,NVT_TRBE)
-END WHERE
 !
 !5 Broadleaf evergreen temperate tree
 ! -------------------------------------
-! utilisation de la classe NVT_EVER pour 30 < LAT < 60 
-! et une hauteur d'arbre supérieur à 3 m. 
-WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).AND.&
-       (ZH_TREE(:,IP_TEBE) .GE. 3.).AND.(ZH_TREE(:,IP_TEBE) .NE. XUNDEF))
 MGN%XPFT(5,:) = S%XVEGTYPE(:,NVT_TEBE)
-END WHERE
-WHERE (((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)).AND.&
-       (ZH_TREE(:,IP_TEBE) .GE. 3.).AND.(ZH_TREE(:,IP_TEBE) .NE. XUNDEF))
-MGN%XPFT(5,:) = S%XVEGTYPE(:,NVT_TEBE)
-END WHERE
 !
 !6 Broadleaf deciduous tropical tree
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour -30 < LAT < 30
-WHERE ((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.))
 MGN%XPFT(6,:) = S%XVEGTYPE(:,NVT_TRBD)
-END WHERE
 !
 !7 Broadleaf deciduous temperate tree
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour 30 < LAT < 60
-! en utilisant une hauteur d'arbre supérieur à 3 m
-WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).AND.&
-       (ZH_TREE(:,IP_TEBD) .GE. 3.).AND.(ZH_TREE(:,IP_TEBD) .NE. XUNDEF))
 MGN%XPFT(7,:) = S%XVEGTYPE(:,NVT_TEBD)
-END WHERE
-WHERE (((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)).AND.&
-       (ZH_TREE(:,IP_TEBD) .GE. 3.).AND.(ZH_TREE(:,IP_TEBD) .NE. XUNDEF))
-MGN%XPFT(7,:) = S%XVEGTYPE(:,NVT_TEBD)
-END WHERE
 !
 !8 Broadleaf deciduous boreal tree
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour LAT > 60
-WHERE (((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)).AND.&
-       (ZH_TREE(:,IP_BOBD) .GE. 3.).AND.(ZH_TREE(:,IP_BOBD) .NE. XUNDEF))
 MGN%XPFT(8,:) = S%XVEGTYPE(:,NVT_BOBD)
-END WHERE
 !
 !9 Broadleaf evergreen shrub
 ! -------------------------------------
-!utilisation de la classe NVT_EVER pour une hauteur d'arbre inférieure à 3 m
-WHERE (ZH_TREE(:,IP_SHRB) .LT. 3.)
+!utilisation de la classe NVT_SHBR pour -30 < LAT < 30
+WHERE (((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.)))
 MGN%XPFT(9,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(9,:) = 0.
 END WHERE
 !
 !10 Broadleaf deciduous temperate shrub
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour une hauteur d'arbre inférieure à 3 m
-! et  pour 30 < LAT < 60
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
-       ((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)))
-MGN%XPFT(10,:) = S%XVEGTYPE(:,NVT_SHRB)
-END WHERE
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
+!utilisation de la classe NVT_SHBR pour 30 < LAT < 60
+WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).OR.&
        ((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)))
 MGN%XPFT(10,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(10,:) = 0.
 END WHERE
 !
 !11 Broadleaf deciduous boreal_shrub
 ! -------------------------------------
-!utilisation de la classe NVT_TREE pour une hauteur d'arbre inférieure à 3 m
-! et  pour  LAT > 60
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
-       ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)))
+!utilisation de la classe NVT_SHBR pour LAT > 60
+WHERE (((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)))
 MGN%XPFT(11,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(11,:) = 0.
 END WHERE
 !
 !12 C3 arctic grass
 ! -------------------------------------
-!utilisation de la classe NVT_GRAS + NVT_PARK pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
-MGN%XPFT(12,:) = S%XVEGTYPE(:,NVT_GRAS) + S%XVEGTYPE(:,NVT_PARK)
-ELSEWHERE
+MGN%XPFT(12,:) = S%XVEGTYPE(:,NVT_BOGR) 
 !
 !13 C3 non-arctic grass
 ! -------------------------------------
-!utilisation de la classe NVT_GRAS + NVT_PARK ailleur
-MGN%XPFT(13,:) = S%XVEGTYPE(:,NVT_GRAS) +  S%XVEGTYPE(:,NVT_PARK)
-END WHERE
+MGN%XPFT(13,:) = S%XVEGTYPE(:,NVT_GRAS)
 !
 !14 C4 grass
 ! -------------------------------------
-! utilisation de la classe NVT_TROG 
 MGN%XPFT(14,:) = S%XVEGTYPE(:,NVT_TROG) 
 !
 !15 Corn
 ! -------------------------------------
-! utilisation de la classe NVT_C4
 MGN%XPFT(15,:) = S%XVEGTYPE(:,NVT_C4) 
 !
 !16 Wheat
 ! -------------------------------------
-! utilisation de la classe NVT_C3
+IF (NVT_C3W .NE. 0 ) THEN ! use ecoclimap_sg
+MGN%XPFT(16,:) = S%XVEGTYPE(:,NVT_C3W) + S%XVEGTYPE(:,NVT_C3S)
+ELSE  ! use ecaclimap2.0
 MGN%XPFT(16,:) = S%XVEGTYPE(:,NVT_C3)
+END IF
 !
 ! Emission factor
 MGN%XEF(:,:) = 0.
@@ -499,7 +455,6 @@ DO JSV=1, MSF%NMEGAN_NBR
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFBIDER")  MGN%XEF(18,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFSTRESS") MGN%XEF(19,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFOTHER")  MGN%XEF(20,:) = PMEGAN_FIELDS(:,JSV)
-!  IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "LAI")      PLAI(:,1)     = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT1")     MGN%XPFT(1,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT2")     MGN%XPFT(2,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT3")     MGN%XPFT(3,:) = PMEGAN_FIELDS(:,JSV)
@@ -516,6 +471,7 @@ DO JSV=1, MSF%NMEGAN_NBR
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT14")    MGN%XPFT(14,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT15")    MGN%XPFT(15,:) = PMEGAN_FIELDS(:,JSV)
   IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT16")    MGN%XPFT(16,:) = PMEGAN_FIELDS(:,JSV)
+  IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "LAI")      MGN%XLAI(:)    = PMEGAN_FIELDS(:,JSV)
 END DO
 
 #endif
diff --git a/src/SURFEX/modn_isban.F90 b/src/SURFEX/modn_isban.F90
index 2ea5770de36f09326b5d2219eb5302efe3a8cf43..8ea9b57710e4a013a2375ad9e0975d24e016b5a1 100644
--- a/src/SURFEX/modn_isban.F90
+++ b/src/SURFEX/modn_isban.F90
@@ -108,7 +108,7 @@ LOGICAL  :: LSNOWDRIFT_SUBLIM
 LOGICAL  :: LSNOW_ABS_ZENITH
  CHARACTER(3) :: CSNOWMETAMO
  CHARACTER(3) :: CSNOWRAD
- CHARACTER(LEN=6)  :: CCH_DRY_DEP, CPARAMBVOC
+ CHARACTER(LEN=6)  :: CCH_DRY_DEP, CPARAMBVOC=''
  CHARACTER(LEN=28) :: CCHEM_SURF_FILE
 !
 NAMELIST/NAM_ISBAn/CC1DRY,CSCOND,CSOILFRZ,CDIFSFCOND,CSNOWRES,CCPSURF, &
diff --git a/src/configure b/src/configure
index c0845c0d5f5e52acd0d7cca1d31c9316e03c1c79..bf9ce3128cccb7730631c1ec4bb936abd15aeaf5 100755
--- a/src/configure
+++ b/src/configure
@@ -19,7 +19,7 @@ export VERSION_CDFCXX=${VERSION_CDFCXX:-"4.3.1"}
 export VERSION_CDFF=${VERSION_CDFF:-"4.5.3"}
 export VERSION_GRIBAPI=${VERSION_GRIBAPI:-"1.26.0-Source"}
 export VERSION_ECCODES=${VERSION_ECCODES:-"2.18.0"}
-export ECCODES_DEFINITION_PATH=${ECCODES_DEFINITION_PATH:${SRC_MESONH}/src/LIB/eccodes-${VERSION_ECCODES}"/definitions/"}
+export ECCODES_DEFINITION_PATH=${ECCODES_DEFINITION_PATH:-${SRC_MESONH}/src/LIB/eccodes-${VERSION_ECCODES}"/definitions/"}
 export MNH_INT=${MNH_INT:-"4"}
 export LFI_INT=${LFI_INT:-8}
 export MNH_REAL=${MNH_REAL:-"8"}
@@ -347,7 +347,7 @@ export CC=gcc
                 export OPTLEVEL=${OPTLEVEL:-DEBUG}
                 export MVWORK=${MVWORK:-NO}
                 export VER_CDF=${VER_CDF:-CDFCTI}
-		export NEED_NCARG=${NEED_NCARG:-NO}
+		export NEED_NCARG=${NEED_NCARG:-YES}
 		export NEED_TOOLS=NO
               ;;
 'Linux nuwa'*)