diff --git a/DESCRIPTION b/DESCRIPTION
index 6c0292d248a0916ea53ca73c5c75b5c8d0a30d3b..c4c04a3826ed30a45254f1c6d67eb0f656be6876 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: SWATdoctR
 Type: Package
 Title: Finding the right diagnoses and treatments for SWAT+ models
-Version: 0.1.13
+Version: 0.1.14
 Author: c(person("Svajunas", "Plunge",
              email = "svajunas_plunge@sggw.edu.pl",
              role = c("aut")),
diff --git a/NAMESPACE b/NAMESPACE
index 9a651e517ee4dafb9b8cb11c80220ff09ef71a26..44d53bee96333ff38ae27da39d8f34a88b784a65 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,7 @@
 
 export(add_kill_op)
 export(check_hru_waterbalance)
+export(filter_hru_at_harvkill)
 export(get_hru_id_by_attribute)
 export(plot_basin_var)
 export(plot_climate_annual)
diff --git a/R/plot_mgt_harv.R b/R/plot_mgt_harv.R
index fdc9d7102b5f3e67c9fd75793e383e171ba1c863..ee22f57d8d0b12de29c7244fa1cc6be0d645ed67 100644
--- a/R/plot_mgt_harv.R
+++ b/R/plot_mgt_harv.R
@@ -147,3 +147,43 @@ var_at_harv <- function(mgt_out, years){
   }
   return(tbl_harv)
 }
+
+
+#' Filter HRUs by their variable values at harvest/kill
+#'
+#' The function can be complementary to \code{plot_variable_at_harvkill()} plots.
+#' If any issues are identified in the boxplots, the results can be filtered for
+#' further analysis with \code{filter_hru_at_harvkill()}.
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}
+#' as a single object or as a list of objects (representing \code{run_swat_verification()}
+#' function runs with different settings).
+#'   To plot the heat units at least the output option \code{outputs = 'mgt'}
+#'   must be set in \code{run_swat_verification()}.
+#' @param ... Boolean operations to be applied in the filter operation of the
+#'   management simulation outputs. Possible variables are \code{crop, phu, plant_bioms,
+#'   yield, water_stress, aero_stress, temp_stress, n_stress} and \code{p_stress}.
+#'
+#' @importFrom dplyr filter mutate select %>%
+#' @importFrom lubridate ymd
+#' @importFrom purrr set_names
+#'
+#' @export
+#' @examples
+#' \dontrun{
+#' # Filter all HRUs with unusually high PHU values at harvest kill
+#' filter_hru_at_harvkill(sim_verify, phu > 5)
+#'
+#' # Filter all HRUs with crop 'wbar' or 'wira' planted but no yield
+#' filter_hru_at_harvkill(sim_verify, crop %in% c('wbar', 'wira'), yield == 0)
+#' }
+filter_hru_at_harvkill <- function(sim_verify, ...) {
+  sim_verify$mgt_out %>%
+    mutate(date = ymd(paste(year, mon, day, sep = '-'))) %>%
+    filter(operation == 'HARVEST') %>%
+    select(hru, date, op_typ, phuplant, plant_bioms, op_var, var4, var5, var3, var1, var2) %>%
+    set_names(c('hru', 'date', 'crop', 'phu', 'plant_bioms', 'yield', 'water_stress', 'aero_stress',
+                'temp_stress', 'n_stress', 'p_stress')) %>%
+    filter(., ...)
+}
+
diff --git a/man/filter_hru_at_harvkill.Rd b/man/filter_hru_at_harvkill.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..69c626a4e7fda7c1be3a25f088b8cc2b9b13f152
--- /dev/null
+++ b/man/filter_hru_at_harvkill.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{filter_hru_at_harvkill}
+\alias{filter_hru_at_harvkill}
+\title{Filter HRUs by their variable values at harvest/kill}
+\usage{
+filter_hru_at_harvkill(sim_verify, ...)
+}
+\arguments{
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}
+as a single object or as a list of objects (representing \code{run_swat_verification()}
+function runs with different settings).
+  To plot the heat units at least the output option \code{outputs = 'mgt'}
+  must be set in \code{run_swat_verification()}.}
+
+\item{...}{Boolean operations to be applied in the filter operation of the
+management simulation outputs. Possible variables are \code{crop, phu, plant_bioms,
+yield, water_stress, aero_stress, temp_stress, n_stress} and \code{p_stress}.}
+}
+\description{
+The function can be complementary to \code{plot_variable_at_harvkill()} plots.
+If any issues are identified in the boxplots, the results can be filtered for
+further analysis with \code{filter_hru_at_harvkill()}.
+}
+\examples{
+\dontrun{
+# Filter all HRUs with unusually high PHU values at harvest kill
+filter_hru_at_harvkill(sim_verify, phu > 5)
+
+# Filter all HRUs with crop 'wbar' or 'wira' planted but no yield
+filter_hru_at_harvkill(sim_verify, crop \%in\% c('wbar', 'wira'), yield == 0)
+}
+}