diff --git a/.Rbuildignore b/.Rbuildignore
index a70e3dd6992644edce14079cf603a6d5d21e9730..c877860597bb70297205ca9c7a8cb1e739012adf 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,4 +1,5 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 test*
-template*
\ No newline at end of file
+template*
+env*
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 19c339f2b32399327830f516133905e62023dfbd..28e29458f67c6512edd0047ce6dfb75e386ccb7d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,4 +28,5 @@ vignettes/*.pdf
 rsconnect/
 .Rproj.user
 test*
-template*
\ No newline at end of file
+template*
+env*
\ No newline at end of file
diff --git a/R/plot_ps_tile.R b/R/plot_ps_tile.R
index 7d6b2870e96bd9a1b098da61df54bb9daa355ea3..78624c883ae71cfc28f10a2f0493b36c2b5fb1bb 100644
--- a/R/plot_ps_tile.R
+++ b/R/plot_ps_tile.R
@@ -118,8 +118,8 @@ plot_qtile <- function(sim_verify, exclude_lum = c("urhd_lum", "urmd_lum", "urml
     left_join(., sim_verify$lum_mgt, by = "id") %>%
     filter(tile != 'null')
 
-  hru_frac <- paste0(as.character(round(100*count(df)/count(sim_verify$hru_wb_aa),2)), " % HRUs are drained in this catchment.")
-
+  hru_frac <- paste0("The plot was created for HRUs with tile drainage active (",
+                      as.character(round(100*count(df)/count(sim_verify$hru_wb_aa),2)), "% of all HRUs).")
   df <- df %>%
     filter(!lu_mgt %in% exclude_lum) %>%
     select(id, qtile, lu_mgt, mgt, soil) %>%
@@ -129,8 +129,7 @@ plot_qtile <- function(sim_verify, exclude_lum = c("urhd_lum", "urmd_lum", "urml
     geom_histogram(aes(y=after_stat(density)), color="black", fill="blue", breaks = seq(min(df$qtile), max(df$qtile), 10))+
     geom_density(alpha=.3, fill="white", linewidth = 1, color = "grey25", linetype = "twodash")+
     labs(title = "Tile drain flow density mm/year",
-         subtitle = hru_frac,
-         caption = "Data in the figure represents only HRUs with active tile drainage.")+
+         subtitle = hru_frac)+
     theme_bw()+
     theme(panel.border = element_blank(),
           axis.line = element_line(color='black'),
diff --git a/template.Rmd b/template.Rmd
index 3ef97746bd7749e4757bb959528836aea4dd7e96..6aaa302210768dfe93f4a1da722ba91ce194fe73 100644
--- a/template.Rmd
+++ b/template.Rmd
@@ -119,7 +119,7 @@ Management operation inputs in a SWAT+ model setup can be very complex and compr
 Function `report_mgt()` can be applied to identify discrepancies between management operations in model input files and what operations are actually triggered in the model. If `write_report` parameter is set to TRUE, function also provides a report in *"schedule_report.txt"* text file. 
 
 ```{r}
-mgt_report <- report_mgt(sim_nostress, TRUE)
+mgt_report <- report_mgt(sim_nostress)
 mgt_report
 print(paste("Issues were identified in", length(mgt_report$schedule), "schedules."))
 ```
@@ -136,6 +136,8 @@ if(length(mgt_report$schedule)>=sel_nb){
   ##Print selected case into interactive table
   print(paste("Table of issues for selected management", sel_mgt))
   create_dt(mgt_report$schedule_report[[sel_nb]])
+} else {
+   id <- get_hru_id_by_attribute(sim_nostress)
 }
 ```
 
@@ -211,17 +213,19 @@ The next step includes activating potential sources for plant growth stresses, s
 Setting `nostress = 0` while running `run_swat_verification()` function will activate all stresses, however turning off the nutrient plant stress only can as well be a useful option for analyses (`nostress = 2`).  This is particularly useful for eliminating the fertilization impact on the plant growth and focusing only on the weather/climate and structural setting of the plant. Particularly, the aeration, temperature, and water stress, alongside yields are relevant outputs to be analyzed. A simulation with inactive nutrient stress will provide a good approximation of possible yields with an optimal fertilization and ideal plant nutrient supply. All other stresses will indicate the need of irrigation, drainage or plant-specific parameter adjustments for a plant to grow. 
 
 ```{r, include=FALSE}
-rm(sim_nostress)
-sim_stress_nutrient <- readRDS(file = rpath[3])
+sim_except_nutrient <- readRDS(file = rpath[3])
 sim_stress_all <- readRDS(file = rpath[1])
+sim_list <- list(no_stress = sim_nostress["mgt_out"],
+                 except_nutrient = sim_except_nutrient["mgt_out"],
+                 stress_all = sim_stress_all["mgt_out"])
+rm(sim_nostress)
 ```
 
 ### Examine stress factors 
 
 It is possible to plot each case side-by-side for examination with same functions applied. For instance `plot_variable_at_harvkill()` could be run to check how much stress factors affect in each case. 
 
-```{r fig.width = 10, fig.height = 6}
-plot_variable_at_harvkill(sim_stress_nutrient, variable = 'stress')
+```{r fig.width = 15, fig.height = 10}
 plot_variable_at_harvkill(sim_stress_all, variable = 'stress')
 ```
 
@@ -230,24 +234,16 @@ plot_variable_at_harvkill(sim_stress_all, variable = 'stress')
 Additionally, we can look how plant growth is different in the same HRUs. 
 
 ```{r fig.width = 10, fig.height = 6}
-plot_hru_pw(sim_stress_nutrient,  id$id[1], var = c('lai', 'bioms'))
+plot_hru_pw(sim_except_nutrient,  id$id[1], var = c('lai', 'bioms'))
 plot_hru_pw(sim_stress_all,  id$id[1], var = c('lai', 'bioms'))
 ```
 
-### Assess difference at the harvest for PHU
-
-```{r fig.width = 10, fig.height = 6}
-plot_variable_at_harvkill(sim_stress_nutrient, variable = 'phu')
-plot_variable_at_harvkill(sim_stress_all, variable = 'phu')
-```
-
 ### Assess difference for yields
 
 Could be useful to examine how harvest is affected by stresses factors. 
 
 ```{r fig.width = 10, fig.height = 6}
-plot_variable_at_harvkill(sim_stress_nutrient, variable = 'yield')
-plot_variable_at_harvkill(sim_stress_all, variable = 'yield')
+plot_variable_at_harvkill(sim_list, variable = 'yield')
 ```
 
 ## Simulated point sources and tile drains