diff --git a/plotting/layout.R b/plotting/layout.R
new file mode 100644
index 0000000000000000000000000000000000000000..a7ab458f0d985218d58ea3fd740d15be94ef43c1
--- /dev/null
+++ b/plotting/layout.R
@@ -0,0 +1,213 @@
+# Usefull library
+library(ggplot2)
+library(scales)
+library(qpdf)
+library(gridExtra)
+library(gridtext)
+library(dplyr)
+library(grid)
+library(ggh4x)
+library(RColorBrewer)
+
+
+# Sourcing R file
+source('plotting/panel.R', encoding='latin1')
+
+
+panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_opt='', filename_opt='', variable='', df_trend=NULL, p_threshold=0.1, unit2day=365.25, type='', period=NULL, missRect=FALSE, time_header=NULL, info_header=TRUE, header_ratio=2) {
+    
+    if (all(class(df_data) != 'list')) {
+        df_data = list(df_data)
+    }
+
+    nbp = length(df_data)
+
+    if (all(class(df_trend) != 'list')) {
+        df_trend = list(df_trend)
+        if (length(df_trend) == 1) {
+            df_trend = replicate(nbp, df_trend)
+        }}
+
+    if (all(class(p_threshold) != 'list')) {
+        p_threshold = list(p_threshold)
+        if (length(p_threshold) == 1) {
+            p_threshold = replicate(nbp, p_threshold)
+        }}
+  
+    if (all(class(unit2day) != 'list')) {
+        unit2day = list(unit2day)
+        if (length(unit2day) == 1) {
+            unit2day = replicate(nbp, unit2day)
+        }}
+
+    if (all(class(type) != 'list')) {
+        type = list(type)
+        if (length(type) == 1) {
+            type = replicate(nbp, type)
+        }}
+
+    if (all(class(missRect) != 'list')) {
+        missRect = list(missRect)
+        if (length(missRect) == 1) {
+            missRect = replicate(nbp, missRect)
+        }}
+
+    list_df2plot = vector(mode='list', length=nbp)
+    minTrend = c()
+    maxTrend = c()
+
+    for (i in 1:nbp) {
+        
+        df2plot = list(data=df_data[[i]], 
+                       trend=df_trend[[i]],
+                       p_threshold=p_threshold[[i]],
+                       unit2day=unit2day[[i]],
+                       type=type[[i]],
+                       missRect=missRect[[i]])
+        
+        okTrend = df_trend[[i]]$trend[df_trend[[i]]$p <= p_threshold[[i]]]
+
+        minTrend[i] = min(okTrend, na.rm=TRUE)
+        maxTrend[i] = max(okTrend, na.rm=TRUE)
+
+        list_df2plot[[i]] = df2plot
+    }
+
+
+    outfile = "Panels"
+    if (filename_opt != '') {
+        outfile = paste(outfile, '_', filename_opt, sep='')
+    }
+    outfile = paste(outfile, '.pdf', sep='')
+
+    # If there is not a dedicated figure directory it creats one
+    outdir = file.path(figdir, filedir_opt, sep='')
+    if (!(file.exists(outdir))) {
+        dir.create(outdir)
+    }
+
+    outdirTmp = file.path(outdir, 'tmp')
+    if (!(file.exists(outdirTmp))) {
+        dir.create(outdirTmp)
+    }
+
+
+    # Get all different stations code
+    Code = levels(factor(df_meta$code))
+    nCode = length(Code)
+    
+    for (code in Code) {
+        
+        # Print code of the station for the current plotting
+        print(paste("Plotting for sation :", code))
+        
+        nbh = as.numeric(info_header) + as.numeric(!is.null(time_header))
+        nbg = nbp + nbh
+
+        P = vector(mode='list', length=nbg)
+
+        if (info_header) {
+            Htext = text_panel(code, df_meta)
+            P[[1]] = Htext
+        }
+
+        if (!is.null(time_header)) {
+            time_header_code = time_header[time_header$code == code,]
+
+            Htime = time_panel(time_header_code, df_trend_code=NULL,
+                               period=period, missRect=TRUE,
+                               unit2day=365.25, type='Q')
+
+            P[[2]] = Htime
+        }
+
+
+        nbcol = ncol(as.matrix(layout_matrix))
+        for (i in 1:nbp) {
+            df_data = list_df2plot[[i]]$data
+            df_trend = list_df2plot[[i]]$trend
+            p_threshold = list_df2plot[[i]]$p_threshold
+            unit2day = list_df2plot[[i]]$unit2day
+            missRect = list_df2plot[[i]]$missRect
+            type = list_df2plot[[i]]$type
+            
+            df_data_code = df_data[df_data$code == code,] 
+            df_trend_code = df_trend[df_trend$code == code,]
+
+            if (df_trend_code$p <= p_threshold){
+                color_res = get_color(df_trend_code$trend, 
+                                      minTrend[i],
+                                      maxTrend[i], 
+                                      palette_name='perso',
+                                      reverse=FALSE)
+
+                color = color_res$color
+                palette = color_res$palette
+
+            } else {            
+                color = NULL
+                palette = NULL
+            }
+            
+            p = time_panel(df_data_code, df_trend_code, type=type,
+                           p_threshold=p_threshold, missRect=missRect,
+                           unit2day=unit2day, last=(i > nbp-nbcol),
+                           color=color)
+
+            P[[i+nbh]] = p
+
+        }
+        
+        layout_matrix = as.matrix(layout_matrix)
+        nel = nrow(layout_matrix)*ncol(layout_matrix)
+
+        idNA = which(is.na(layout_matrix), arr.ind=TRUE)
+
+        layout_matrix[idNA] = seq(max(layout_matrix, na.rm=TRUE) + 1,
+                                  max(layout_matrix, na.rm=TRUE) + 1 +
+                                  nel)
+
+        layout_matrix_H = layout_matrix + nbh
+
+
+        LM = c()
+        LMcol = ncol(layout_matrix_H)
+        LMrow = nrow(layout_matrix_H)
+        for (i in 1:(LMrow+nbh)) {
+
+            if (i <= nbh) {
+                LM = rbind(LM, rep(i, times=LMcol))
+            } else {
+                LM = rbind(LM, 
+                           matrix(rep(layout_matrix_H[i-nbh,],
+                                      times=header_ratio),
+                                  ncol=LMcol, byrow=TRUE))
+            }}
+
+        plot = grid.arrange(grobs=P, layout_matrix=LM)
+        
+        # plot = grid.arrange(rbind(cbind(ggplotGrob(P[[2]]), ggplotGrob(P[[2]])), cbind(ggplotGrob(P[[3]]), ggplotGrob(P[[3]]))), heights=c(1/3, 2/3))
+        
+
+        # Saving
+        ggsave(plot=plot, 
+               path=outdirTmp,
+               filename=paste(as.character(code), '.pdf', sep=''),
+               width=21, height=29.7, units='cm', dpi=100)
+
+    }
+
+    # mat = matrice_panel(list_df2plot, df_meta)
+
+    # # Saving matrix plot
+    # ggsave(plot=mat, 
+    #        path=outdirTmp,
+    #        filename=paste('matrix', '.pdf', sep=''),
+    #        width=21, height=29.7, units='cm', dpi=100)
+
+    # PDF combine
+    pdf_combine(input=file.path(outdirTmp, list.files(outdirTmp)),
+                output=file.path(outdir, outfile))
+    unlink(outdirTmp, recursive=TRUE)
+
+} 
diff --git a/plotting/panel.R b/plotting/panel.R
index 655f562a9bcd1f658faa8f9599fbc477aaef013a..835d9a7769f6852a04961cc79c852309c6537f72 100644
--- a/plotting/panel.R
+++ b/plotting/panel.R
@@ -184,10 +184,10 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
 
     # if (norm) {
     #     p = p +
-    #         ylab(bquote('débit ['*m^{3}*'.'*s^{-1}*']  x'~10^{.(as.character(power))}))
+    #         ylab(bquote('débit ['*m^{3}*'.'*s^{-1}*']  x'~10^{.(as.character(power))}))
     # } else {
     #     p = p +
-    #         ylab(expression(paste('débit [', m^{3}, '.', 
+    #         ylab(expression(paste('débit [', m^{3}, '.', 
     #                               s^{-1}, ']', sep='')))
     # }
 
@@ -218,7 +218,7 @@ text_panel = function(code, df_meta) {
     text = paste(
         "<span style='font-size:18pt'> station <b>", code, "</b></span><br>",
         "nom : ", df_meta_code$nom, "<br>", 
-        "région hydrographique : ", df_meta_code$region_hydro, "<br>",
+        "région hydrographique : ", df_meta_code$region_hydro, "<br>",
         "position : (", df_meta_code$L93X, "; ", df_meta_code$L93Y, ")", "<br>",
         "surface : ", df_meta_code$surface_km2, " km<sup>2</sup>",
         sep='')
diff --git a/script.R b/script.R
index b684f2732af24004fdb046ca7ec286c77da90746..56575dc1ccedfb7c5953ea64dee6dc4b402b1dc1 100644
--- a/script.R
+++ b/script.R
@@ -43,13 +43,13 @@ BHlistname =
 ### NIVALE ###
 # Path to the directory where NV data is stored
 NVfiledir = 
-    ""
-    # "France207"
+    # ""
+    "France207"
 
 # Name of the file that will be analysed from the NV directory
 NVfilename = 
-    ""
-    # "all"
+    # ""
+    "all"
 
 
 # Path to the list file of metadata about station that will be analysed
@@ -57,14 +57,13 @@ NVlistdir =
     ""
 
 NVlistname = 
-    ""
-    # "liste_bv_principaux_global.txt"
+    # ""
+    "liste_bv_principaux_global.txt"
 
 
 ### TREND ANALYSIS ###
 # Time period to analyse
 period = c("1980-01-01","2019-12-31")
-# period = c("1960-01-01","2020-01-01")
 
 
 
@@ -157,7 +156,7 @@ df_meta = df_join$meta
 res_QAtrend = get_QAtrend(df_data, period)
 
 # QMNA TREND #
-# res_QMNAtrend = get_QMNAtrend(df_data, period)
+res_QMNAtrend = get_QMNAtrend(df_data, period)
 
 # VCN10 TREND #
 res_VCN10trend = get_VCN10trend(df_data, df_meta, period)
@@ -176,29 +175,14 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta, period)
 #               figdir=figdir,
 #               filename_opt='time')
 
-# panels_layout(list(res_QAtrend$data, res_QMNAtrend$data,
-                   # res_VCN10trend$data), 
-              # layout_matrix=c(1, 2, 3),
-              # df_meta=df_meta, 
-              # df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend,
-                    # res_VCN10trend$trend), 
-              # type=list(bquote(Q[A]), bquote(Q[MNA]), bquote(V[CN10])),
-              # missRect=list(TRUE, TRUE, TRUE),
-              # period=period,
-              # info_header=TRUE,
-              # time_header=df_data,
-              # header_ratio=2,
-              # figdir=figdir,
-              # filename_opt='')
-
-
-panels_layout(list(res_QAtrend$data, res_VCN10trend$data), 
-              layout_matrix=c(1, 2),
+panels_layout(list(res_QAtrend$data, res_QMNAtrend$data,
+                   res_VCN10trend$data), 
+              layout_matrix=c(1, 2, 3),
               df_meta=df_meta, 
-              df_trend=list(res_QAtrend$trend,
+              df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend,
                     res_VCN10trend$trend), 
-              type=list(bquote(Q[A]), bquote(V[CN10])),
-              missRect=list(TRUE, TRUE),
+              type=list(bquote(Q[A]), bquote(Q[MNA]), bquote(V[CN10])),
+              missRect=list(TRUE, TRUE, TRUE),
               period=period,
               info_header=TRUE,
               time_header=df_data,
@@ -207,6 +191,21 @@ panels_layout(list(res_QAtrend$data, res_VCN10trend$data),
               filename_opt='')
 
 
+# panels_layout(list(res_QAtrend$data, res_VCN10trend$data), 
+#               layout_matrix=c(1, 2),
+#               df_meta=df_meta, 
+#               df_trend=list(res_QAtrend$trend,
+#                     res_VCN10trend$trend), 
+#               type=list(bquote(Q[A]), bquote(V[CN10])),
+#               missRect=list(TRUE, TRUE),
+#               period=period,
+#               info_header=TRUE,
+#               time_header=df_data,
+#               header_ratio=2,
+#               figdir=figdir,
+#               filename_opt='')
+
+
 ### /!\ Removed 185 row(s) containing missing values (geom_path) -> remove NA ###