diff --git a/DESCRIPTION b/DESCRIPTION
index ff706076409153e6b77291b649794aecdc2a6509..571142ae65a3a238b8faa2a390713723a65f6417 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,6 +19,7 @@ Imports:
   processx,
   purrr,
   readr,
+  stats,
   stringr,
   tibble,
   tidyr,
@@ -26,8 +27,7 @@ Imports:
   plotly,
   utils
 Suggests:
-  datasets,
-  stats
+  datasets
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
diff --git a/NAMESPACE b/NAMESPACE
index 5a3b7ff8dc152a7129dac37c7dc357690e294d40..920a3a793191e41d320ce911f4601a587c13d46d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,10 +5,13 @@ export(check_hru_waterbalance)
 export(get_hru_id_by_attribute)
 export(plot_basin_var)
 export(plot_climate_annual)
+export(plot_hru_pw)
 export(plot_hru_pw_day)
 export(plot_hru_var)
 export(plot_hru_var_aa)
 export(plot_monthly_snow)
+export(plot_ps)
+export(plot_qtile)
 export(plot_variable_at_harvkill)
 export(plot_water_partition)
 export(print_avannual_qtile)
@@ -22,9 +25,14 @@ import(patchwork)
 importFrom(data.table,fread)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,arrange)
+importFrom(dplyr,bind_rows)
 importFrom(dplyr,case_when)
+importFrom(dplyr,count)
 importFrom(dplyr,distinct)
+importFrom(dplyr,do)
+importFrom(dplyr,ends_with)
 importFrom(dplyr,filter)
+importFrom(dplyr,full_join)
 importFrom(dplyr,group_by)
 importFrom(dplyr,group_map)
 importFrom(dplyr,group_split)
@@ -35,6 +43,7 @@ importFrom(dplyr,mutate_all)
 importFrom(dplyr,mutate_at)
 importFrom(dplyr,n)
 importFrom(dplyr,rename)
+importFrom(dplyr,rename_with)
 importFrom(dplyr,row_number)
 importFrom(dplyr,select)
 importFrom(dplyr,slice)
@@ -44,12 +53,17 @@ importFrom(dplyr,summarise_all)
 importFrom(dplyr,tibble)
 importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
+importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_boxplot)
+importFrom(ggplot2,geom_density)
+importFrom(ggplot2,geom_histogram)
 importFrom(ggplot2,geom_hline)
+importFrom(ggplot2,geom_line)
 importFrom(ggplot2,geom_text)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,labs)
 importFrom(ggplot2,lims)
+importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_bw)
 importFrom(ggplot2,theme_void)
 importFrom(lubridate,floor_date)
@@ -78,7 +92,9 @@ importFrom(purrr,map_int)
 importFrom(purrr,map_lgl)
 importFrom(purrr,set_names)
 importFrom(readr,read_lines)
+importFrom(readr,write_delim)
 importFrom(readr,write_lines)
+importFrom(stats,density)
 importFrom(stringr,str_detect)
 importFrom(stringr,str_remove)
 importFrom(stringr,str_remove_all)
diff --git a/R/globals.R b/R/globals.R
index 9487d505feeed880fc388bdf318311852da0fa52..ed56251343bce44a6a4648bd32b49f4e0ce253cf 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -6,4 +6,4 @@ utils::globalVariables(c("%&&%", "%//%", ".", "crop", "day", "ecanopy", "eplant"
                          "tile", "time_out", "tmn", "tmpav", "tmx", "topo", "val_max", "val_mean",
                          "val_min", "value", "value_sum", "var", "var1", "var2", "var3", "var4", "var5",
                          "wndspd", "yr", "rhum", "rm_skp", "Date", "Values", "surq_gen", "latq", "perc",
-                         "description"))
+                         "description", "flo", "..density.."))
diff --git a/R/helper.R b/R/helper.R
index fde56bb2eb404a98e0442ba3dcab190fddf7424b..ad06264dc9400221fbff1da00290a3f7f063a175 100644
--- a/R/helper.R
+++ b/R/helper.R
@@ -75,3 +75,32 @@ hide_show <- function(graph){
                   args = list("visible", "legendonly"),
                   label = "hide all")))))
 }
+
+#' Remove tail endings for management, plant community, land use classes
+#'
+#' This function is required for SWATfarmR generated management files
+#' as management.sch, plant.ini, landuse.lum and hru-data.hru can get too
+#' long names for swat executable.
+#'
+#' @param f multiline character
+#' @param pattern pattern after which to remove tail
+#' @return corrected multiline character
+#' @keywords internal
+#'
+#' @examples
+#' \dontrun{
+#' library(readr)
+#' landuse <- read_lines(paste0(project_path,'/landuse.lum'), lazy = FALSE)
+#' landuse  <- remove_tail(landuse, "lum")
+#' landuse  <- remove_tail(landuse, "comm")
+#' landuse  <- remove_tail(landuse, "mgt")
+#' write_lines(landuse, paste0(project_path,'/landuse_new.lum'))
+#' }
+
+remove_tail <- function(f, pattern){
+  ind <- grep(paste0("_", pattern,"_"), f)
+  for(i in ind){
+    f[i] <- gsub(paste0("(_", pattern, ").*?\\s"), "\\1 ", f[i])
+  }
+  return(f)
+}
diff --git a/R/plot_hru_pw.R b/R/plot_hru_pw.R
index fbc84662a5f2b4b943e1dbd6b7186182edbbb8c6..e21529e1462a0b6f1665c7302128c4c019a10968 100644
--- a/R/plot_hru_pw.R
+++ b/R/plot_hru_pw.R
@@ -79,38 +79,6 @@ plot_hru_pw_day <- function(sim_verify, hru_id, var, title = "", years = 1900:21
           axis.title.y = element_blank())
 }
 
-#' Print the average annual qtile for HRUs
-#'
-#' print_avannual_qtile prints a table with the average annual qtile in mm
-#' for HRUs that used a tile flow parametrization in landuse.lum
-#'
-#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
-#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
-#'   be set in  \code{run_swat_verification()}
-#' @param exclude_lum Character vector to define land uses which are excluded
-#'   in the printed table.
-#'
-#' @importFrom dplyr arrange filter left_join rename select %>%
-#'
-#' @return Returns a table with hru ids average annual qtile and attributes.
-#'
-#' @export
-#'
-print_avannual_qtile <- function(sim_verify,
-                                 exclude_lum = c(
-                                   "urhd_lum", "urmd_lum", "urml_lum",
-                                   "urld_lum", "ucom_lum", "uidu_lum",
-                                   "utrn_lum", "uins_lum", "urbn_lum"
-                                 )) {
-
-  sim_verify$hru_wb_aa %>%
-    rename(id = unit) %>%
-    left_join(., sim_verify$lum_mgt, by = "id") %>%
-    filter(tile != 'null') %>%
-    filter(!lu_mgt %in% exclude_lum) %>%
-    select(id, qtile, lu_mgt, mgt, soil) %>%
-    arrange(qtile, id)
-}
 
 #' Aggregate and plot simulated variables saved in hru_pw_day
 #'
@@ -173,7 +141,11 @@ plot_hru_var <- function(sim_verify, hru_id, var, period = "day", fn_summarize =
 #' }
 
 plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
-  p <- paste(lum, mgt, soil)
+  p <- if(is.null(lum) & is.null(mgt) & is.null(soil)){
+    p <- "all"
+  } else {
+    p <- paste0(lum,"|",mgt,"|",soil)
+  }
   id <- get_hru_id_by_attribute(sim_verify, lum, mgt, soil)
   fig <- sim_verify$hru_wb_aa[sim_verify$hru_wb_aa$unit %in% id$id, -c(1:7)] %>%
     .[, colSums(.!= 0) > 0] %>%
@@ -182,7 +154,7 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
     group_by(p, var) %>%
     group_map(~plot_ly(., y=~Values, color = ~var, colors = "cyan4", type = 'box'), keep = TRUE) %>%
     subplot(nrows = 5) %>%
-    layout(title = paste("HRUs selected by", p))
+    layout(title = paste0("HRUs selected: ", p))
   return(fig)
 }
 
@@ -197,6 +169,8 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
 #' @param soil Optional character vector with soil type labels as defined in the soil data.
 #' @param exclude_lum Character vector to define land uses which are excluded
 #'   in the printed table.
+#' @param boxpoints Optional Boolean TRUE for displaying outliers, FALSE for hiding them.
+#'  \code{Default = TRUE}
 #' @return plotly figure object
 #' @importFrom dplyr %>% mutate group_by rename left_join summarise_all filter select
 #' @importFrom tidyr pivot_longer
@@ -211,7 +185,7 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
 plot_water_partition <- function(sim_verify, tile = NULL, lum = NULL, mgt = NULL, soil = NULL, exclude_lum = c(
   "urhd_lum", "urmd_lum", "urml_lum",
   "urld_lum", "ucom_lum", "uidu_lum",
-  "utrn_lum", "uins_lum", "urbn_lum")){
+  "utrn_lum", "uins_lum", "urbn_lum"), boxpoints = TRUE){
   df <- sim_verify$hru_wb_aa %>%
     rename(id = unit) %>%
     left_join(., sim_verify$lum_mgt, by = "id") %>%
@@ -261,7 +235,7 @@ plot_water_partition <- function(sim_verify, tile = NULL, lum = NULL, mgt = NULL
     add_pie(hole = 0.3)
   ##Preparing box plot
   box_pl <- plot_ly(df[c("var", "Values")], x=~Values,  color = ~var,  type = "box",  colors = pal,
-                    showlegend = F) %>%
+                    showlegend = F,  boxpoints = boxpoints) %>%
     layout(yaxis = list(autorange = "reversed"))
   ##Putting into one figure and annotations
   fig <- subplot(box_pl, pie_pl, nrows = 1, margin = 0.05) %>%
@@ -277,3 +251,62 @@ plot_water_partition <- function(sim_verify, tile = NULL, lum = NULL, mgt = NULL
   options(warn = -1)
   return(fig)
 }
+
+#' Aggregate and plot simulated variables saved in hru_pw_day
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'mgt'} must
+#'   be set in  \code{run_swat_verification()}.
+#' @param hru_id Numeric vector with HRU ids for which variables should be plotted.
+#' @param var Character vector that defines the variable names that are plotted.
+#' @param title Character for title to be put in the figure.
+#' @param period (optional) character describing, which time interval to display (default is "day",
+#' other examples are "week", "month", etc). \code{Default = "day"}
+#' @param fn_summarize (optional) function to recalculate to time interval (default is "mean", other examples
+#' are "median", "sum", etc). \code{Default = "mean"}
+#' @param interactive (optional) Boolean TRUE to provide output as plotly object, FALSE - as ggplot.  \code{Default = "FALSE"}
+#' @return ggplot or plotly object
+#' @importFrom lubridate floor_date
+#' @importFrom plotly plot_ly layout subplot
+#' @importFrom dplyr %>% rename summarise mutate group_by arrange do
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' id <- get_hru_id_by_attribute(sim_nostress, "wwht_lum")
+#' plot_hru_pw(sim_nostress, sample(id$id, 10), "lai")
+#' plot_hru_pw(sim_nostress, hru_id = c(2, 10), var = c("lai", "bioms"),
+#' title = "My great figure", period = "year", fn_summarize = "mean", interactive = TRUE)
+#' }
+
+plot_hru_pw <- function(sim_verify, hru_id, var,  title = "", period = "day", fn_summarize = "mean", interactive = FALSE){
+  options(dplyr.summarise.inform = FALSE)
+  df <- sim_verify$hru_pw_day[sim_verify$hru_pw_day$unit %in% hru_id, c("unit", "yr", "mon", "day", var)]
+  ##Aggregating data by time step
+  df$Date<- floor_date(ISOdate(df$yr, df$mon, df$day), period)
+  df <- df[c("Date", "unit", var)] %>%
+    mutate(unit = paste('hru:', unit)) %>%
+    pivot_longer(., cols = - c(Date, unit), names_to = 'var', values_to = 'Values') %>%
+    group_by(unit, Date, var) %>%
+    summarise(Values = get(fn_summarize)(Values))
+  if(interactive){
+    fig <- df %>%
+      group_by(var) %>%
+      do(p=plot_ly(., x = ~Date, y = ~Values, color = ~factor(unit), text = ~var, textposition = 'outside', colors = "Set2", type = 'scatter', mode = 'lines', connectgaps = FALSE) %>%
+           layout(showlegend = FALSE, xaxis = list(title = "Date"), yaxis = list(title = paste(var, "variable values")))) %>%
+      subplot(nrows = ifelse(length(var) %/% 1.75<1, 1, length(var) %/% 1.75), shareX = TRUE, shareY = FALSE, titleY= TRUE, margin = c(0.1,0,0.1,0)) %>%
+      layout(title = title)
+  } else {
+    fig <- ggplot(df) +
+      geom_line(aes(x = Date, y = Values, color = unit, lty = unit)) +
+      labs(x = 'Date', color = 'HRU', lty = 'HRU', title=title) +
+      facet_grid(rows = vars(all_of(var)), scales = 'free_y', switch = 'y') +
+      theme_bw() +
+      theme(legend.position = 'bottom',
+            strip.background = element_blank(),
+            strip.placement = 'outside',
+            strip.text = element_text(face = 'bold'),
+            axis.title.y = element_blank())
+  }
+  return(fig)
+}
diff --git a/R/plot_ps_tile.R b/R/plot_ps_tile.R
new file mode 100644
index 0000000000000000000000000000000000000000..7d6b2870e96bd9a1b098da61df54bb9daa355ea3
--- /dev/null
+++ b/R/plot_ps_tile.R
@@ -0,0 +1,140 @@
+#' Plot yearly simulated point source values (as yearly loads or average concentrations)
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param conc Boolean, TRUE to provide figure for point source average yearly pollutant concentrations,
+#' FALSE - for average yearly pollutant loads.
+#' @return ggplot object for point source yearly simulated data
+#' @importFrom dplyr %>% mutate group_by group_map select bind_rows case_when
+#' @importFrom ggplot2 ggplot geom_line facet_wrap labs theme_bw theme aes
+#' @importFrom tidyr pivot_longer
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' plot_ps(sim_nostress, TRUE)
+#' }
+
+plot_ps <- function(sim_verify, conc = FALSE){
+  if(!is.null(sim_verify$recall_yr)){
+    df <- sim_verify$recall_yr[, -c(1:3,6)] %>%
+      .[, colSums(.!= 0) > 0] %>%
+      mutate(name = gsub("hru00", "ps", name))
+    if(conc){
+      df <- df %>%
+        pivot_longer(-c(yr, name, flo), names_to = 'var', values_to = 'Values') %>%
+        mutate(Values = ifelse(var == "sed", 1000000*Values/flo, 1000*Values/flo)) %>%
+        select(yr, name, var, Values) %>%
+        bind_rows(df[c("yr", "name", "flo")] %>%
+                    mutate(var = "flo") %>%
+                    rename(Values = flo) %>%
+                    mutate(Values = Values/(24*60*60*365.25))) %>%
+        mutate(var = case_when(var == 'flo' ~ "flo m3/s",
+                               var %in% c("orgn", "no3", "nh3", "no2") ~ paste(var, "N mg/l"),
+                               var %in% c("sedp", "solp") ~ paste(var, "P mg/l"),
+                               TRUE ~  paste(var, "mg/l")))
+    } else {
+      df <- df %>%
+        pivot_longer(-c(yr, name), names_to = 'var', values_to = 'Values') %>%
+        mutate(var = case_when(var == 'flo' ~ "flo m3/d",
+                               var == 'sed' ~ "sed t/y",
+                               var %in% c("orgn", "no3", "nh3", "no2") ~ paste(var, "N kg/y"),
+                               var %in% c("sedp", "solp") ~ paste(var, "P kg/y"),
+                               TRUE ~  paste(var, "kg/y")))
+    }
+    fig <- ggplot(df, aes(x=yr, y=Values,  group=name, colour=name))+
+      geom_line(linewidth=1)+
+      facet_wrap(~var, scales = "free_y")+
+      labs(color='Point sources', x = 'Year') +
+      theme_bw()+
+      theme(strip.background = element_rect(fill = "grey50", colour = "grey80"),
+            strip.text = element_text(color = "white", face="bold"),
+            panel.border = element_rect(colour = "grey80"),
+            axis.text.x = element_text(angle = 25, hjust=1))
+
+    return(fig)
+  } else {
+    print("No point sources exists in this model setup!!!")
+  }
+}
+
+#' Print the average annual qtile for HRUs
+#'
+#' print_avannual_qtile prints a table with the average annual qtile in mm
+#' for HRUs that used a tile flow parametrization in landuse.lum
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param exclude_lum Character vector to define land uses which are excluded
+#'   in the printed table.
+#'
+#' @importFrom dplyr arrange filter left_join rename select %>%
+#'
+#' @return Returns a table with hru ids average annual qtile and attributes.
+#'
+#' @export
+#'
+print_avannual_qtile <- function(sim_verify,
+                                 exclude_lum = c(
+                                   "urhd_lum", "urmd_lum", "urml_lum",
+                                   "urld_lum", "ucom_lum", "uidu_lum",
+                                   "utrn_lum", "uins_lum", "urbn_lum"
+                                 )) {
+
+  sim_verify$hru_wb_aa %>%
+    rename(id = unit) %>%
+    left_join(., sim_verify$lum_mgt, by = "id") %>%
+    filter(tile != 'null') %>%
+    filter(!lu_mgt %in% exclude_lum) %>%
+    select(id, qtile, lu_mgt, mgt, soil) %>%
+    arrange(qtile, id)
+}
+
+#' Plot tile drain flow histogram and distribution curve
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param exclude_lum Character vector to define land uses which are excluded
+#'   in the printed table.
+#' @importFrom dplyr count
+#' @importFrom stats density
+#' @importFrom ggplot2 ggplot geom_histogram labs theme_bw theme geom_density aes
+#' @return Returns ggplot object with histogram and density curve for tile drain from distribution.
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' plot_qtile(sim_nostress)
+#' }
+
+plot_qtile <- function(sim_verify, exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum",
+                                                   "urld_lum", "ucom_lum", "uidu_lum",
+                                                   "utrn_lum", "uins_lum", "urbn_lum")){
+  df <- sim_verify$hru_wb_aa %>%
+    rename(id = unit) %>%
+    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.")
+
+  df <- df %>%
+    filter(!lu_mgt %in% exclude_lum) %>%
+    select(id, qtile, lu_mgt, mgt, soil) %>%
+    arrange(qtile, id)
+
+  fig <- ggplot(df, aes(x=qtile)) +
+    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.")+
+    theme_bw()+
+    theme(panel.border = element_blank(),
+          axis.line = element_line(color='black'),
+          axis.title.x=element_blank())
+
+  return(fig)
+}
diff --git a/R/print_mgt.R b/R/print_mgt.R
index 7bff2eb4652853716ffd9d20f15b6c8dc943ca3b..af00e2453544c4009ae853c499d5e6ed3b5b21b5 100644
--- a/R/print_mgt.R
+++ b/R/print_mgt.R
@@ -29,8 +29,7 @@ print_triggered_mgt <- function(sim_verify, hru_id, years = 1900:2100) {
     mutate(op_data3 = ifelse(operation != 'FERT', 0, op_data3)) %>%
     # mutate()
     # mutate(date = ymd(paste(year, mon, day, sep = '-'))) %>%
-    select(., year, mon, day, phuplant, operation, op_data1, op_data3) %>%
-    print(., n = Inf)
+    select(., year, mon, day, phuplant, operation, op_data1, op_data3)
 }
 
 #' Generate a report table that compares the scheduled and triggered managements
@@ -42,18 +41,21 @@ print_triggered_mgt <- function(sim_verify, hru_id, years = 1900:2100) {
 #' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
 #'   To print the management at least the output option \code{outputs = 'mgt'} must
 #'   be set in  \code{run_swat_verification()}
+#' @param write_report (optional) Boolean TRUE for writing output to 'schedule_report.txt' file,
+#' FALSE - not preparing this file. Default \code{write_report = FALSE}.
 #'
 #' @return Returns a tibble that summarises all management schedules for
 #'   which operations where scheduled, that were either not triggered of
 #'   for which operation properties differ..
 #'
-#' @importFrom dplyr filter group_by group_split lead left_join mutate rename select slice_sample summarise %>%
+#' @importFrom dplyr filter group_by group_split lead left_join mutate rename select slice_sample summarise %>% rename_with full_join ends_with
 #' @importFrom purrr map
-#' @importFrom stringr str_sub
+#' @importFrom stringr str_sub str_remove
+#' @importFrom readr write_delim write_lines
 #'
 #' @export
 #'
-report_mgt <- function(sim_verify) {
+report_mgt <- function(sim_verify, write_report = FALSE) {
   yr_start <- min(sim_verify$mgt_out$year)
   mgt_lbl <- unique(sim_verify$lum_mgt$mgt)
   mgt_lbl <- mgt_lbl[!is.na(mgt_lbl)]
@@ -90,12 +92,41 @@ report_mgt <- function(sim_verify) {
     mutate(op_typ = str_sub(op_typ, 1, 4) %>%  tolower(.),
            op_typ = ifelse(op_typ == 'plan', 'plnt', op_typ))
 
-  schdl_join <- left_join(schdl_mgt, mgt_i,
-                          by = c("schedule", "year",  "mon", "day", "op_typ")) %>%
+  ##Case for PHU scheduling
+  if(any(schdl_mgt$hu_sch>0)){
+    schdl_mgt_in <- schdl_mgt %>%
+      group_by(schedule, year) %>%
+      mutate(year = year - yr_start + 1) %>%
+      mutate(id = row_number())
+
+    df <- NULL
+    for(sch in unique(schdl_mgt_in$schedule)){ ##For each schedule in schdl_mgt
+      sch1 <- schdl_mgt_in[schdl_mgt_in$schedule == sch,] ##data separated
+      counter_max <- max(sch1$year) ##counter set
+      counter <- counter_min <- 1
+      for (n in seq(yr_start,max(sim_verify$mgt_out$year),1)){ ##For each modeling year
+        d_copy <- sch1[sch1$year == counter,]
+        d_copy$year <- n
+        if(!is.null(df)){df <- bind_rows(df, d_copy)}else{df <- d_copy}
+        if(counter>=counter_max){counter<-counter_min}else{counter<-counter+1} ##Counter for schedules with more than 1 year
+      }
+    }
+    schdl_join <- full_join(df, mgt_i %>%
+                              group_by(schedule, year) %>%
+                              mutate(id = row_number()),
+                            by = c("schedule", "year", "id", "op_typ", "op_data1" = "op_data1_trig"), keep = TRUE)
+  } else {
+    schdl_join <- full_join(schdl_mgt, mgt_i,
+                            by = c("schedule", "year",  "mon", "day", "op_typ", "op_data1" = "op_data1_trig"), keep = TRUE)
+  }
+
+  schdl_join <- schdl_join %>%
+    select(-ends_with(".y")) %>%
+    rename_with(~str_remove(., '.x')) %>%
     select(schedule, year, mon, day, op_typ, op_data1_trig, starts_with('op_data')) %>%
     mutate(op_issue = is.na(op_data1_trig) | op_data1_trig != op_data1,
            year = year - yr_start + 1) %>%
-    filter(op_issue)
+    filter(op_issue & year <= max(sim_verify$mgt_out$year) - yr_start)
 
   schdl_report <- schdl_join %>%
     select(schedule, op_issue) %>%
@@ -112,6 +143,20 @@ report_mgt <- function(sim_verify) {
   schdl_report <- schdl_report %>%
     mutate(schedule_report = ops_detail)
 
+  if(write_report){
+    print("Writing schedule_report.txt")
+    write(paste("File was written with SWATdoctR at", Sys.time( )), file = "schedule_report.txt")
+    for (i in seq(1, length(schdl_report$schedule_report))){
+      mgt <- schdl_report$schedule[[i]]
+      id <- get_hru_id_by_attribute(sim_verify, mgt = mgt)$id[1]
+      write_lines(" ", "schedule_report.txt", append = TRUE)
+      write_lines(paste("HRU number -", id, "- management name:", mgt), "schedule_report.txt", append = TRUE)
+      write_lines(" ", "schedule_report.txt", append = TRUE)
+      write_delim(schdl_report$schedule_report[[i]], "schedule_report.txt", delim = "\t", append = TRUE, col_names = TRUE)
+    }
+    print(paste("The file schedule_report.txt was written in", getwd( ), "directory."))
+  }
+
   return(schdl_report)
 }
 
diff --git a/R/run_swat_verify.R b/R/run_swat_verify.R
index a2bdfd3aa6993429c1d9965a32fa0f73632952e2..120042ab06cdc60d5d252e45ed80af6b6f2c4833 100644
--- a/R/run_swat_verify.R
+++ b/R/run_swat_verify.R
@@ -31,7 +31,7 @@
 #' @return Returns the simulation results for the defined output variables as a
 #'   list of tibbles.
 #'
-#' @importFrom dplyr %>%
+#' @importFrom dplyr %>% distinct
 #' @importFrom processx run
 #' @importFrom stringr str_split
 #' @importFrom tibble tibble
@@ -69,7 +69,6 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
     model_output <- err_msg
   } else if(nchar(msg$stderr) == 0) {
     model_output <- list()
-
     if ('plt' %in% outputs) {
       model_output$hru_pw_day <- read_tbl('hru_pw_day.txt', run_path, 3)
     }
@@ -77,6 +76,12 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
       model_output$basin_wb_day <- read_tbl('basin_wb_day.txt', run_path, 3)
       model_output$basin_pw_day <- read_tbl('basin_pw_day.txt', run_path, 3)
       model_output$hru_wb_aa <- read_tbl('hru_wb_aa.txt', run_path, 3)
+      tryCatch({
+        model_output$recall_yr <- read_tbl('recall_yr.txt', run_path, 3)
+      },
+      error = function(e) {
+        model_output$recall_yr <- NULL
+      })
     }
     if ('mgt' %in% outputs) {
       model_output$mgt_out <- read_mgt(run_path)
@@ -86,6 +91,7 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
       landuse_lum <- read_tbl('landuse.lum', run_path, 2)
       model_output$lum_mgt <- left_join(hru_data, landuse_lum,
                                         by = c("lu_mgt" = 'name')) %>%
+        distinct() %>%
         select(id, topo, hydro, soil, lu_mgt, plnt_com, mgt, tile)
     }
   }
diff --git a/R/update_input_files.R b/R/update_input_files.R
index 102605896e0b15df2c12704f2485335336abc8af..ba2733a8a6d7aad95f81a9deae9b11872cb1c67f 100644
--- a/R/update_input_files.R
+++ b/R/update_input_files.R
@@ -50,6 +50,7 @@ set_print_prt <- function(project_path, run_path, outputs, years_skip) {
     print_prt[11] <- "basin_wb                     y             n             n             n  "
     print_prt[14] <- "basin_pw                     y             n             n             n  "
     print_prt[33] <- "hru_wb                       n             n             n             y  "
+    print_prt[45] <- "recall                       n             n             y             y  "
   }
   if ('mgt' %in% outputs) {
     print_prt[9]  <- "n             y             n             n             "
diff --git a/man/plot_hru_pw.Rd b/man/plot_hru_pw.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..ec8a80a647be569821f25b6d738d87331401c123
--- /dev/null
+++ b/man/plot_hru_pw.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_hru_pw.R
+\name{plot_hru_pw}
+\alias{plot_hru_pw}
+\title{Aggregate and plot simulated variables saved in hru_pw_day}
+\usage{
+plot_hru_pw(
+  sim_verify,
+  hru_id,
+  var,
+  title = "",
+  period = "day",
+  fn_summarize = "mean",
+  interactive = FALSE
+)
+}
+\arguments{
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}.
+To plot the heat units at least the output option \code{outputs = 'mgt'} must
+be set in  \code{run_swat_verification()}.}
+
+\item{hru_id}{Numeric vector with HRU ids for which variables should be plotted.}
+
+\item{var}{Character vector that defines the variable names that are plotted.}
+
+\item{title}{Character for title to be put in the figure.}
+
+\item{period}{(optional) character describing, which time interval to display (default is "day",
+other examples are "week", "month", etc). \code{Default = "day"}}
+
+\item{fn_summarize}{(optional) function to recalculate to time interval (default is "mean", other examples
+are "median", "sum", etc). \code{Default = "mean"}}
+
+\item{interactive}{(optional) Boolean TRUE to provide output as plotly object, FALSE - as ggplot.  \code{Default = "FALSE"}}
+}
+\value{
+ggplot or plotly object
+}
+\description{
+Aggregate and plot simulated variables saved in hru_pw_day
+}
+\examples{
+\dontrun{
+id <- get_hru_id_by_attribute(sim_nostress, "wwht_lum")
+plot_hru_pw(sim_nostress, sample(id$id, 10), "lai")
+plot_hru_pw(sim_nostress, hru_id = c(2, 10), var = c("lai", "bioms"),
+title = "My great figure", period = "year", fn_summarize = "mean", interactive = TRUE)
+}
+}
diff --git a/man/plot_ps.Rd b/man/plot_ps.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..dba94e8de919358c510b335ea9fd6c3093df2a99
--- /dev/null
+++ b/man/plot_ps.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_ps_tile.R
+\name{plot_ps}
+\alias{plot_ps}
+\title{Plot yearly simulated point source values (as yearly loads or average concentrations)}
+\usage{
+plot_ps(sim_verify, conc = FALSE)
+}
+\arguments{
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}.
+To plot the heat units at least the output option \code{outputs = 'wb'} must
+be set in  \code{run_swat_verification()}}
+
+\item{conc}{Boolean, TRUE to provide figure for point source average yearly pollutant concentrations,
+FALSE - for average yearly pollutant loads.}
+}
+\value{
+ggplot object for point source yearly simulated data
+}
+\description{
+Plot yearly simulated point source values (as yearly loads or average concentrations)
+}
+\examples{
+\dontrun{
+plot_ps(sim_nostress, TRUE)
+}
+}
diff --git a/man/plot_qtile.Rd b/man/plot_qtile.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..046fd0def141ed8c40dcad5db8bde9ea75bc5b50
--- /dev/null
+++ b/man/plot_qtile.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_ps_tile.R
+\name{plot_qtile}
+\alias{plot_qtile}
+\title{Plot tile drain flow histogram and distribution curve}
+\usage{
+plot_qtile(
+  sim_verify,
+  exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum", "urld_lum", "ucom_lum", "uidu_lum",
+    "utrn_lum", "uins_lum", "urbn_lum")
+)
+}
+\arguments{
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}.
+To plot the heat units at least the output option \code{outputs = 'wb'} must
+be set in  \code{run_swat_verification()}}
+
+\item{exclude_lum}{Character vector to define land uses which are excluded
+in the printed table.}
+}
+\value{
+Returns ggplot object with histogram and density curve for tile drain from distribution.
+}
+\description{
+Plot tile drain flow histogram and distribution curve
+}
+\examples{
+\dontrun{
+plot_qtile(sim_nostress)
+}
+}
diff --git a/man/plot_water_partition.Rd b/man/plot_water_partition.Rd
index cba081aa17783a88120bcdbb2de2f440c828f5d5..4f9d63542cf60d9e1b76cdba54b024d00002ea46 100644
--- a/man/plot_water_partition.Rd
+++ b/man/plot_water_partition.Rd
@@ -11,7 +11,8 @@ plot_water_partition(
   mgt = NULL,
   soil = NULL,
   exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum", "urld_lum", "ucom_lum", "uidu_lum",
-    "utrn_lum", "uins_lum", "urbn_lum")
+    "utrn_lum", "uins_lum", "urbn_lum"),
+  boxpoints = TRUE
 )
 }
 \arguments{
@@ -29,6 +30,9 @@ be set in  \code{run_swat_verification()}}
 
 \item{exclude_lum}{Character vector to define land uses which are excluded
 in the printed table.}
+
+\item{boxpoints}{Optional Boolean TRUE for displaying outliers, FALSE for hiding them.
+\code{Default = TRUE}}
 }
 \value{
 plotly figure object
diff --git a/man/print_avannual_qtile.Rd b/man/print_avannual_qtile.Rd
index 4474d8ac91d6c7c8f5b15f26eb22fb0ffa506893..ed1a39b520e70e7921ca0f24f54017335a730b3f 100644
--- a/man/print_avannual_qtile.Rd
+++ b/man/print_avannual_qtile.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot_hru_pw.R
+% Please edit documentation in R/plot_ps_tile.R
 \name{print_avannual_qtile}
 \alias{print_avannual_qtile}
 \title{Print the average annual qtile for HRUs}
diff --git a/man/remove_tail.Rd b/man/remove_tail.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..faa1f4051195e2afc2f2f171a9337868b8461d72
--- /dev/null
+++ b/man/remove_tail.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper.R
+\name{remove_tail}
+\alias{remove_tail}
+\title{Remove tail endings for management, plant community, land use classes}
+\usage{
+remove_tail(f, pattern)
+}
+\arguments{
+\item{f}{multiline character}
+
+\item{pattern}{pattern after which to remove tail}
+}
+\value{
+corrected multiline character
+}
+\description{
+This function is required for SWATfarmR generated management files
+as management.sch, plant.ini, landuse.lum and hru-data.hru can get too
+long names for swat executable.
+}
+\examples{
+\dontrun{
+library(readr)
+landuse <- read_lines(paste0(project_path,'/landuse.lum'), lazy = FALSE)
+landuse  <- remove_tail(landuse, "lum")
+landuse  <- remove_tail(landuse, "comm")
+landuse  <- remove_tail(landuse, "mgt")
+write_lines(landuse, paste0(project_path,'/landuse_new.lum'))
+}
+}
+\keyword{internal}
diff --git a/man/report_mgt.Rd b/man/report_mgt.Rd
index 467be71e298f5878fe229608841c5c29a330d5e6..357c6d5fa2772ce178e475982ffd0178e57e73c5 100644
--- a/man/report_mgt.Rd
+++ b/man/report_mgt.Rd
@@ -4,12 +4,15 @@
 \alias{report_mgt}
 \title{Generate a report table that compares the scheduled and triggered managements}
 \usage{
-report_mgt(sim_verify)
+report_mgt(sim_verify, write_report = FALSE)
 }
 \arguments{
 \item{sim_verify}{Simulation output of the function \code{run_swat_verification()}.
 To print the management at least the output option \code{outputs = 'mgt'} must
 be set in  \code{run_swat_verification()}}
+
+\item{write_report}{(optional) Boolean TRUE for writing output to 'schedule_report.txt' file,
+FALSE - not preparing this file. Default \code{write_report = FALSE}.}
 }
 \value{
 Returns a tibble that summarises all management schedules for
diff --git a/template.Rmd b/template.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..3ef97746bd7749e4757bb959528836aea4dd7e96
--- /dev/null
+++ b/template.Rmd
@@ -0,0 +1,284 @@
+<style type="text/css">
+.main-container {
+  max-width: 100% !important;
+  margin: auto;
+}
+</style>
+
+---
+title: "Zglowiaczka catchment (PL)"
+author: "name surname"
+affiliation: "your affiliation"
+email: "your@email"
+date: "`r Sys.Date()`"
+output:
+  html_document:
+    df_print: paged
+---
+
+This is a document template for providing report for setup verification workflow and its output in html format. 
+
+## Running model
+
+Before starting working with `SWATdoctR` functions it is necessary to set paths to model to be verified and a folder, which could be used to save large files generated by the `run_swat_verification()` function.  
+
+```{r}
+library(SWATdoctR)
+##Set SWAT+ model path txtitout folder
+model_path <- "path_to_txtinout"
+##This path to save run_swat_verification function results
+result_path <- "path_to_run_swat_verification_results"
+model_name <- "model_a"
+
+model_path <- "../Test_models/Zglowiaczka"
+result_path <- "./test_data/"
+model_name <- "Zglowiaczka"
+rpath <- c()
+for(ns in c(0:2)){rpath <- c(rpath, paste0(result_path, model_name, "_ns", ns, ".rds"))}
+```
+
+Running `run_swat_verification()` function is important for `SWATdoctR` package as outputs of this function is used by other package functions. Code below runs function with different plant stress options activated and saves results in `.rds` files for the later examination. This is important as this step could be very time consuming (requires 3 model runs). Saving function results allows examining results without rerunning `run_swat_verification()` function. However, if the model was changed (by adjusting parameters, correcting/new input data, etc) function should be rerun.
+
+```{r, eval = FALSE}
+##Fixing management operations (changing hvkl to harv + kill). If there is not such issues nothing will be changed.
+add_kill_op(model_path)
+##Running loop to get and save run_swat_verification function output for different plant stress options
+for(ns in c(0:2)){
+  saveRDS(run_swat_verification(project_path = model_path, outputs = c('wb', 'mgt', 'plt'), 
+                                nostress = ns, keep_folder = TRUE), file = rpath[[ns+1]])}
+```
+
+Before starting model verification `run_swat_verification()`results should be loaded in the memory. The first two steps require function results without plant stresses activate (`nostress = 1`).
+
+```{r}
+sim_nostress <- readRDS(file = rpath[2])
+```
+
+```{r setup, include=FALSE}
+##Here adding some functions and libraries to be used in template
+library(DT)
+##Function to deliver table 
+create_dt <- function(table, pL=10){
+  DT::datatable(table,
+                extensions = 'Buttons',
+                fillContainer = FALSE, 
+                options = list(pageLength = pL, 
+                  initComplete = JS(
+                  "function(settings, json) {",
+                  "$(this.api().table().header()).css({'background-color': '#000', 'color': '#fff', 'font-size': '80%'});",
+                  "}"),
+                  dom = 'Blfrtip',
+                  buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
+                  lengthMenu = list(c(pL,25,50,-1),
+                                    c(pL,25,50,"All")),
+                  columnDefs = list(list(className = 'dt-center', targets="_all"))),
+                height = "100%")
+}
+```
+
+## Step 1: Climate & Water flows
+
+### Analysis of simulated annual climate variables
+
+The first step in verification workflow is to check climatic variables simulated by model. Multiple mistakes (as wrong units, mistake in inputs, etc) could cause a lot of problems in model. `plot_climate_annual()` function provides quick overview of climatic variables. 
+
+```{r fig.width = 20, fig.height = 7}
+plot_climate_annual(sim_nostress)
+```
+
+### Seasonal dynamics of the precipitation
+
+`plot_monthly_snow()` function generates overview for precipitation partition into snow fall and rain as well as provides summary of snowmelt monthly averages for simulation period.
+
+```{r fig.width = 12, fig.height = 10}
+plot_monthly_snow(sim_nostress)
+```
+
+### Water partition
+
+`plot_water_partition` function generates interactive figure plotting water pathway distribution after precipitation reached the ground.  
+
+```{r, out.width="80%"}
+plot_water_partition(sim_nostress)
+```
+
+### HRU variable distribution
+
+For detail examination of HRU variables `plot_hru_var_aa()` function could be applied with possibility to filter HRU by landuse (`lum` parameter in function), management (`mgt`) and soil (`soil`) labels. 
+
+```{r, out.width="80%"}
+plot_hru_var_aa(sim_nostress)
+```
+
+## Step 2: Management
+
+### Examining issues with management
+
+Management operation inputs in a SWAT+ model setup can be very complex and comprehensive (i.e. could include multiple years, plants, operations, even set to each field separately, etc). Scheduled operations point to several other input files which define the parameters of operations or inputs such as fertilizer or tillage types. Hence, the development of management schedules is highly error prone. Mistakes in management inputs usually do not stop a simulation or produce warnings in the model diagnostics, but lead to skipping certain operations in a simulation run.
+
+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
+print(paste("Issues were identified in", length(mgt_report$schedule), "schedules."))
+```
+
+### Case to examine in tables
+
+The function `report_mgt()` generates list of possible issues for management types. To examine them is detail on case-by-case basis following lines could be applied to receive case specific information and table about which management operations were triggered in model.
+
+```{r}
+##Select an example case to examine
+sel_nb <- 2
+if(length(mgt_report$schedule)>=sel_nb){
+  sel_mgt <- mgt_report[[sel_nb,"schedule"]]
+  ##Print selected case into interactive table
+  print(paste("Table of issues for selected management", sel_mgt))
+  create_dt(mgt_report$schedule_report[[sel_nb]])
+}
+```
+
+```{r}
+if(length(mgt_report$schedule)>=sel_nb){
+  ##Get HRU ids
+  id <- get_hru_id_by_attribute(sim_nostress, mgt = sel_mgt)
+  ##Print triggered management table
+  print(paste("HRU", id$id[1], "trigerred management table"))
+  create_dt(print_triggered_mgt(sim_nostress, id$id[1]))
+}
+```
+
+### Case to examine in figures
+
+Selected case can be examined in figures plotting selected HRU variables. Following variables are available for plotting.
+
+- *lai m^2^/m^2^ |average leaf area index during timestep*
+- *bioms kg/ha |average total plant biomass during timestep*
+- *yield kg/ha |harvested biomass yield (dry weight) during timestep*
+- *residue kg/ha |average surface residue cover during timestep*
+- *sol_tmp deg C |average temperature of soil layer 2 during timestep*
+- *strsw days |limiting water (drought) stress*
+- *strsa days |excess water (aeration) stress*
+- *strstmp days |temperature stress*
+- *strsn days |nitrogen stress*
+- *strsp days |phosphorus stress*
+- *nplnt kg N/ha |plant uptake of nitrogen*
+- *percn kg N/ha |nitrate NO3-N leached from bottom of soil profile*
+- *pplnt kg P/ha |plant uptake of phosphorus*
+- *tmx deg C |average maximum temperature during timestep*
+- *tmn deg C |average minimum temperature during timestep*
+- *tmpav deg C |average of average daily air temperature during timestep*
+- *solrad MJ/m^2 |average solar radiation during timestep*
+- *wndspd m/s |average windspeed during timestep*
+- *rhum none |average relative humidity during timestep*
+- *phubase0 deg c |base zero potential heat units*
+
+Function `plot_hru_pw()` can be used for this. It provides static and interactive options as well as aggregation by time interval (week, month, year, etc) and type of aggregation operation (mean, sum, max, etc).
+
+```{r fig.width = 10, fig.height = 6}
+plot_hru_pw(sim_nostress, id$id[1], "lai")
+plot_hru_pw(sim_nostress, id$id[1], c("lai", "bioms", "yield", "residue"), period = "month", fn_summarize = "max", interactive = TRUE)
+```
+
+## Step 3: Plant growth without stess factors
+
+The verification of plant growth is a two-tiered approach. In a first step plant growth is simulated
+and analyzed without simulating any limiting stress factors. Such analysis illustrates the potential biomass or yields, which plants can gain given the climatic and soil conditions of the simulated catchment. Moreover, it allows us to verify the duration of the scheduled growing period or if the selected crop parametrizations meet the climatic conditions.
+
+`plot_variable_at_harvkill()` function can be applied to extract relevant variables at the time of harvesting a crop (available options are 'phu', 'yield', 'bioms'). 
+
+### Check the stress factors
+
+It is all so important to check, if plant stress factors are actually switched off. If correct a figure should be plotting zeros.
+
+```{r fig.width = 10, fig.height = 6}
+plot_variable_at_harvkill(sim_nostress, variable = 'stress')
+```
+
+### Plot variables during harvest
+
+Plotting PHUs (plant heat units) and yield could help to verify plant growth. 
+
+```{r fig.width = 10, fig.height = 6}
+plot_variable_at_harvkill(sim_nostress, variable = 'phu')
+plot_variable_at_harvkill(sim_nostress, variable = 'yield')
+```
+
+### Step 4: Plant growth with stress active
+
+The next step includes activating potential sources for plant growth stresses, such as nutrient stress due to limited fertilizer inputs, or water stress due to limited water availability. An analysis of plant growth with active stresses can show issues in quantities of scheduled operations, such as the amounts of fertilizer inputs, or the problems with definition of irrigation schedules, or decision table rules. 
+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_stress_all <- readRDS(file = rpath[1])
+```
+
+### 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')
+plot_variable_at_harvkill(sim_stress_all, variable = 'stress')
+```
+
+### Check the difference for selected HRU
+
+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_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')
+```
+
+## Simulated point sources and tile drains
+
+The point source time series inputs define the water, sediment, and nutrient loads which are emitted by a point source into a spatial object. The common issues are wrong units of the defined fluxes or wrong time intervals for a certain accumulated flux. SWAT+ projects can potentially have a large number of unique point sources, which may be another reason for errors. Thus, it is good practice to verify the simulated influxes from point sources into the respective spatial objects.
+
+The verification of tile flow can focus on whether tile flow occurs or not. Also could be important to examine the rate of tile flow.
+
+Point source data could be examined using `plot_ps()` function. If `conc = FALSE`, yearly loads for different point sources will be provided in a figure. If `conc = TRUE`, concentration figure will be generated.
+
+
+```{r fig.width = 10, fig.height = 7}
+plot_ps(sim_stress_all)
+```
+
+```{r fig.width = 10, fig.height = 7}
+plot_ps(sim_stress_all, TRUE)
+```
+
+For examining tile drain flow two functions are available in the `SWATdoctR` package. `print_avannual_qtile()` function extracts tile flow for HRUs with installed tile drain systems.
+
+```{r}
+print_avannual_qtile(sim_stress_all)
+```
+
+`plot_qtile()` provides histogram and distribution curve in a figure for tile drain flows (in mm/year units).  
+
+```{r fig.width = 10, fig.height = 7}
+plot_qtile(sim_stress_all)
+```
+
+```{r, include = FALSE}
+rm(sim_stress_nutrient, sim_stress_all)
+```