diff --git a/DESCRIPTION b/DESCRIPTION
index 126e281..4eeb115 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -18,8 +18,7 @@ RoxygenNote: 7.3.3
 Suggests:
   knitr,
   rmarkdown,
-  testthat,
-  vdiffr
+  testthat
 Imports:
   circlize,
   curl,
@@ -30,7 +29,6 @@ Imports:
   grid,
   gridExtra,
   lazyeval,
-  magrittr,
   proj4,
   purrr,
   RColorBrewer,
diff --git a/NAMESPACE b/NAMESPACE
index 829157a..74adaee 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -75,5 +75,4 @@ export(scan_prm)
 export(str_split_twice)
 export(theme_atlantis)
 export(write_diet)
-importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
diff --git a/R/calculate-consumed-biomass.R b/R/calculate-consumed-biomass.R
index 7432d7d..eab5a64 100644
--- a/R/calculate-consumed-biomass.R
+++ b/R/calculate-consumed-biomass.R
@@ -96,7 +96,7 @@ calculate_consumed_biomass <- function(eat, grazing, dm, vol, bio_conv) {
       atoutput = atoutput * bio_conv
     ) |>
     # Step2: Combine with diet contribution. We need a full join to make sure no data is lost!
-    dplyr::full_join(dm, by = c("species" = "pred", "time", "agecl")) %>%
+    dplyr::full_join(dm, by = c("species" = "pred", "time", "agecl")) |>
     # Restrict timesteps to netcdf data! Last timestep is weird in Dietcheck.txt.
     dplyr::filter(time %in% ts_eat) |>
     dplyr::rename(pred = species)
diff --git a/R/calculate-spatial-overlap.R b/R/calculate-spatial-overlap.R
index cc6bd9c..15a1cbe 100644
--- a/R/calculate-spatial-overlap.R
+++ b/R/calculate-spatial-overlap.R
@@ -146,7 +146,7 @@ schoener <- function(predgrp, ageclass, biomass, avail) {
     df_avail,
     df_pred,
     by = c("pred" = "species", "pred_stanza" = "species_stanza")
-  ) %>%
+  ) |>
     dplyr::inner_join(
       biomass_clean,
       by = c(
diff --git a/R/check-growth.R b/R/check-growth.R
index ad7a8b4..dcde491 100644
--- a/R/check-growth.R
+++ b/R/check-growth.R
@@ -28,15 +28,6 @@ check_growth <- function(data, yearly = FALSE) {
     df
   }
 
-  # Divide output with initial value!
-  # ref <- data[data$time == min(data$time), ]
-  # ref$time <- NULL
-  # names(ref)[names(ref) == "atoutput"] <- "atoutput_ref"
-  # result <- data %>%
-  #   dplyr::left_join(ref) %>%
-  #   dplyr::mutate(atoutput = atoutput / atoutput_ref)
-  # result$atoutput[result$atoutput_ref == 0] <- 0
-  # outcomment in case lm procedure is used! This is a bit messy.
   result <- data
 
   # Split dataframe into species and age specific subdataframes!
diff --git a/R/combine-ages.R b/R/combine-ages.R
index 4b98fb9..5908867 100644
--- a/R/combine-ages.R
+++ b/R/combine-ages.R
@@ -43,10 +43,10 @@ combine_ages <- function(data, grp_col, agemat, value_col = "atoutput") {
   data_stanza$stanza <- ifelse(data_stanza$agecl < data_stanza$age_mat, 1, 2)
   data_stanza$stanza[is.na(data_stanza$stanza)] <- 1 # Not sure if this is correct!
 
-  result <- data_stanza %>%
+  result <- data_stanza |>
     agg_data(
       col = value_col,
-      groups = names(.)[!names(.) %in% c(value_col, "agecl")],
+      groups = setdiff(colnames(data_stanza), c(value_col, "agecl")),
       out = value_col,
       fun = sum
     )
diff --git a/R/combine-groups.R b/R/combine-groups.R
index 5f34421..c16ef91 100644
--- a/R/combine-groups.R
+++ b/R/combine-groups.R
@@ -29,8 +29,8 @@ combine_groups <- function(
   # Arrange by group_col and group by groups to select the 1:combinthethresh
   # group in group_col for each grouping combination.
   if (length(groups) > 0) {
-    comb_grps <- comb_grps %>%
-      as.data.frame() %>%
+    comb_grps <- comb_grps |>
+      as.data.frame() |>
       group_data(groups = groups)
   }
   imp_species <- comb_grps |>
@@ -47,7 +47,7 @@ combine_groups <- function(
     # Only combine groups if necessary!
     low_contrib[, group_col] <- "Rest"
 
-    new_data <- dplyr::inner_join(data, imp_species) %>%
+    new_data <- dplyr::inner_join(data, imp_species) |>
       rbind(low_contrib)
 
     # This should never happen...
diff --git a/R/load-init-age.R b/R/load-init-age.R
index 5316bd9..79cab9c 100644
--- a/R/load-init-age.R
+++ b/R/load-init-age.R
@@ -237,7 +237,7 @@ load_init_weight <- function(init, fgs, bboxes) {
     fgs = fgs,
     select_variable = "ResN",
     bboxes = bboxes
-  ) %>%
+  ) |>
     dplyr::filter(!is.na(atoutput)) |>
     dplyr::select(atoutput, species, agecl) |>
     dplyr::rename(rn = atoutput) |>
@@ -247,7 +247,7 @@ load_init_weight <- function(init, fgs, bboxes) {
     fgs = fgs,
     select_variable = "StructN",
     bboxes = bboxes
-  ) %>%
+  ) |>
     dplyr::filter(!is.na(atoutput)) |>
     dplyr::select(atoutput, species, agecl) |>
     dplyr::rename(sn = atoutput) |>
diff --git a/R/load-mort.R b/R/load-mort.R
index 0a6e4ac..de10453 100644
--- a/R/load-mort.R
+++ b/R/load-mort.R
@@ -32,8 +32,8 @@ load_mort <- function(mortFile, prm_run, fgs, convert_names = F) {
   mort <- load_txt(file = mortFile, id_col = c("Time"))
 
   # separate species code.mortalityType into two columns
-  mort <- mort %>%
-    tidyr::separate(.data$code, into = c("code", "source"), sep = "\\.") %>%
+  mort <- mort |>
+    tidyr::separate(.data$code, into = c("code", "source"), sep = "\\.") |>
     tibble::as_tibble()
 
   # First time step only has 0s as entry!
@@ -43,12 +43,11 @@ load_mort <- function(mortFile, prm_run, fgs, convert_names = F) {
   # Convert species codes to longnames!
   if (convert_names) {
     data_fgs <- load_fgs(fgs = fgs)
-    mort <- mort %>%
+    mort <- mort |>
       dplyr::left_join(
-        .,
         data_fgs[, c("Code", "LongName")],
         by = c("code" = "Code")
-      ) %>%
+      ) |>
       dplyr::rename(species = .data$LongName)
   }
 
diff --git a/R/load-spec-mort.R b/R/load-spec-mort.R
index 0ebbd75..20acd45 100644
--- a/R/load-spec-mort.R
+++ b/R/load-spec-mort.R
@@ -39,18 +39,17 @@ load_spec_mort <- function(
     df_txt = df,
     into = c("code", "agecl", "empty_col", "mort"),
     removeZeros = removeZeros
-  ) %>%
+  ) |>
     tibble::as_tibble()
 
   # Convert species codes to longnames!
   if (convert_names) {
     data_fgs <- load_fgs(fgs = fgs)
-    mort <- mort %>%
+    mort <- mort |>
       dplyr::left_join(
-        .,
         data_fgs[, c("Code", "LongName")],
         by = c("code" = "Code")
-      ) %>%
+      ) |>
       dplyr::rename(species = .data$LongName)
   }
 
diff --git a/R/plot-spatial-box.R b/R/plot-spatial-box.R
index 27ee30b..e0d95c7 100644
--- a/R/plot-spatial-box.R
+++ b/R/plot-spatial-box.R
@@ -138,7 +138,7 @@ plot_spatial_box <- function(
   }
 
   # Step2: Apply predator and stanza specific plot function
-  dfs_spatial <- select_time(perc_bio, timesteps = timesteps) %>%
+  dfs_spatial <- select_time(perc_bio, timesteps = timesteps) |>
     split_dfs(cols = c("species", "species_stanza"))
   plots_spatial <- lapply(
     dfs_spatial,
diff --git a/R/sc-init.R b/R/sc-init.R
index 4cdf3e6..d5f670d 100644
--- a/R/sc-init.R
+++ b/R/sc-init.R
@@ -171,14 +171,14 @@ sc_init <- function(
 
   # Comment in to extract numbers from output!
   # nums <- load_nc(dir = dir, nc = nc, bps = bps, select_variable = "Nums",
-  #                 fgs = fgs, select_groups = groups_age, bboxes = bboxes) %>%
+  #                 fgs = fgs, select_groups = groups_age, bboxes = bboxes) |>
   #   dplyr::filter(time == 0)
   nums <- load_init_age(
     init = init,
     fgs = fgs,
     select_variable = "Nums",
     bboxes = bboxes
-  ) %>%
+  ) |>
     dplyr::inner_join(surface, by = c("polygon", "layer")) # not needed in case numbers are already only in surface in init file
 
   # Add stanzas for all age-based groups!
@@ -211,7 +211,7 @@ sc_init <- function(
 
   # Comment in to get data from output file
   # preydens_invert <- load_nc(dir = dir, nc = nc, bps = bps, fgs = fgs, select_groups = groups_rest,
-  #                            select_variable = "N", bboxes = bboxes) %>%
+  #                            select_variable = "N", bboxes = bboxes) |>
   #   dplyr::filter(time == 0)
   preydens_invert <- load_init_nonage(
     init = init,
@@ -239,7 +239,7 @@ sc_init <- function(
     dplyr::rename(prey = species, preydens = atoutput)
 
   # Extract availability matrix and combine with assimilation types
-  ass_type <- fgs_data %>%
+  ass_type <- fgs_data |>
     dplyr::select(Code, dplyr::any_of(c("GroupType", "InvertType")))
 
   names(ass_type) <- c("prey", "grp")
diff --git a/R/pipe.R b/data-raw/pipe.R
similarity index 100%
rename from R/pipe.R
rename to data-raw/pipe.R
diff --git a/man/pipe.Rd b/man/pipe.Rd
deleted file mode 100644
index bab9ab0..0000000
--- a/man/pipe.Rd
+++ /dev/null
@@ -1,11 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/pipe.R
-\name{\%>\%}
-\alias{\%>\%}
-\title{Pipeoperator}
-\description{
-Atlantistools makes heavy use of dplyr data transformations.
-Therefore it is advisable to import the pipeoperator \code{\%>\%}
-from magrittr.
-}
-\keyword{internal}
diff --git a/tests/testthat/test-combine-groups.R b/tests/testthat/test-combine-groups.R
index 9e75732..40a1048 100644
--- a/tests/testthat/test-combine-groups.R
+++ b/tests/testthat/test-combine-groups.R
@@ -1,8 +1,8 @@
 context("combine_groups tests")
 
 wuwu <- combine_groups(ref_dm, group_col = "pred", combine_thresh = 2)
-wawa <- wuwu %>%
-  dplyr::group_by(prey, agecl) %>%
+wawa <- wuwu |>
+  dplyr::group_by(prey, agecl) |>
   dplyr::summarise(count = dplyr::n_distinct(pred))
 
 # no grouping variable!
diff --git a/tests/testthat/test-load-init-age.R b/tests/testthat/test-load-init-age.R
index 87acffc..109a7f4 100644
--- a/tests/testthat/test-load-init-age.R
+++ b/tests/testthat/test-load-init-age.R
@@ -62,9 +62,9 @@ dl <- load_init_stanza(init = init, fgs = fgs, bboxes = bboxes)
 
 dm <- load_init_weight(init = init, fgs = fgs, bboxes = bboxes)
 
-data <- de %>%
-  dplyr::filter(!is.na(atoutput)) %>%
-  dplyr::group_by(species, agecl) %>%
+data <- de |>
+  dplyr::filter(!is.na(atoutput)) |>
+  dplyr::group_by(species, agecl) |>
   dplyr::summarise(out = unique(atoutput))
 
 test_that("test output numbers", {
diff --git a/vignettes/model-calibration.Rmd b/vignettes/model-calibration.Rmd
index d183d80..17a1e6e 100644
--- a/vignettes/model-calibration.Rmd
+++ b/vignettes/model-calibration.Rmd
@@ -21,7 +21,6 @@ to the list of dataframes at the end of the preprocess vignette and change `mode
 library("atlantistools")
 library("ggplot2")
 library("gridExtra")
-library("magrittr")
 
 fig_height2 <- 11
 gen_labels <- list(x = "Time [years]", y = "Biomass [t]")
@@ -112,7 +111,7 @@ update_labels(plot, gen_labels)
 
 ## Biomass benchmark 2
 ```{r}
-plot <- plot_line(result$biomass) %>% update_labels(labels = gen_labels)
+plot <- plot_line(result$biomass) |>  update_labels(labels = gen_labels)
 plot_add_range(plot, ex_bio)
 ```
 
@@ -127,12 +126,13 @@ custom_grid(plot, grid_x = "polygon", grid_y = "variable")
 
 # Physics
 ```{r, results = 'asis', fig.width = 30, fig.height = 12}
-physics <- result$physics %>%
-  flip_layers() %>%
-  split(., .$variable)
+physics <- result$physics |>
+  flip_layers() |>
+  dplyr::group_split(variable)
 
-plots <- lapply(physics, plot_line, wrap = NULL) %>% 
-  lapply(., custom_grid, grid_x = "polygon", grid_y = "layer")
+plots <- physics |> 
+  lapply(plot_line, wrap = NULL) |> 
+  lapply(custom_grid, grid_x = "polygon", grid_y = "layer")
 
 for (i in seq_along(plots)) {
   cat(paste0("## ", names(plots)[i]), sep = "\n")
@@ -145,7 +145,7 @@ for (i in seq_along(plots)) {
 
 ## Fluxes 1
 ```{r, fig.width = 30, fig.height = 12}
-plot <- flip_layers(result$flux) %>% 
+plot <- flip_layers(result$flux) |> 
   plot_line(wrap = NULL, col = "variable")
 custom_grid(plot, grid_x = "polygon", grid_y = "layer")
 ```
@@ -153,16 +153,16 @@ custom_grid(plot, grid_x = "polygon", grid_y = "layer")
 
 ## Fluxes 2
 ```{r, fig.width = 30, fig.height = 12}
-plot <- flip_layers(result$sink) %>% 
+plot <- flip_layers(result$sink) |>  
   plot_line(wrap = NULL, col = "variable")
 custom_grid(plot, grid_x = "polygon", grid_y = "layer")
 ```
 
 ## Relative change of water column height compared to nominal_dz
 ```{r}
-check_dz <- result$dz %>% 
-  dplyr::left_join(result$nominal_dz, by = c("polygon", "layer")) %>% 
-  dplyr::mutate(check_dz = atoutput.x / atoutput.y) %>% 
+check_dz <- result$dz |> 
+  dplyr::left_join(result$nominal_dz, by = c("polygon", "layer")) |> 
+  dplyr::mutate(check_dz = atoutput.x / atoutput.y) |> 
   dplyr::filter(!is.na(check_dz)) # remove sediment layer
 
 plot <- plot_line(check_dz, x = "time", y = "check_dz", wrap = "polygon", col = "layer")
diff --git a/vignettes/model-comparison.Rmd b/vignettes/model-comparison.Rmd
index 2f021a5..192dd5d 100644
--- a/vignettes/model-comparison.Rmd
+++ b/vignettes/model-comparison.Rmd
@@ -19,7 +19,7 @@ to the list of dataframes at the end of the preprocess vignette and change `mode
 library("atlantistools")
 library("ggplot2")
 library("gridExtra")
-library("magrittr")
+
 gen_labels <- list(x = "Time [years]", y = "Biomass [t]")
 
 # You should be able to build the vignette either by clicking on "Knit PDF" in RStudio or with
@@ -61,7 +61,7 @@ update_labels(plot, gen_labels)
 
 # Biomass timeseries
 ```{r}
-plot_line(result$biomass, col = "run", ncol = 4) %>% update_labels(gen_labels)
+plot_line(result$biomass, col = "run", ncol = 4)  |>  update_labels(gen_labels)
 ```
 
 
diff --git a/vignettes/model-preprocess.Rmd b/vignettes/model-preprocess.Rmd
index 56e6ad6..ae9711a 100644
--- a/vignettes/model-preprocess.Rmd
+++ b/vignettes/model-preprocess.Rmd
@@ -24,7 +24,7 @@ library("atlantistools")
 library("ggplot2")
 library("gridExtra")
 library("dplyr")
-library("magrittr")
+
 
 # You should be able to build the vignette either by clicking on "Knit" in RStudio or with
 # rmarkdown::render("model-preprocess.Rmd")
@@ -100,8 +100,8 @@ vol_dz <- load_nc_physics(nc = nc_gen, select_physics = c("volume", "dz"),
 dz <- dplyr::filter(vol_dz, variable == "dz")
 vol <- dplyr::filter(vol_dz, variable == "volume")
 
-nominal_dz <- load_init(init = init, vars = "nominal_dz") %>% 
-  as.data.frame() %>% 
+nominal_dz <- load_init(init = init, vars = "nominal_dz") |>  
+  as.data.frame() |> 
   dplyr::filter(!is.na(layer))
 
 # Read in Dietcheck
@@ -126,11 +126,11 @@ bio_sp <- calculate_biomass_spatial(nums = dfs_gen[[1]], sn = dfs_gen[[2]], rn =
 bio_sp_stanza <- combine_ages(bio_sp, grp_col = "species", agemat = df_agemat)
 
 # Aggregate biomass
-biomass <- bio_sp %>%
+biomass <- bio_sp |> 
   agg_data(groups = c("species", "time"), fun = sum)
 
-biomass_age <- bio_sp %>%
-  filter(agecl > 2) %>%
+biomass_age <- bio_sp |> 
+  filter(agecl > 2) |> 
   agg_data(groups = c("species", "agecl", "time"), fun = sum)
 
 # Aggregate Numbers! This is done seperately since numbers need to be summed!
@@ -147,7 +147,7 @@ growth_age  <- agg_data(data = dfs_prod[[3]], groups = c("species", "time", "age
 
 # Calculate consumed biomass
 bio_cons <- calculate_consumed_biomass(eat = dfs_prod[[1]], grazing = dfs_prod[[2]], dm = df_dm,
-                                       vol = vol, bio_conv = bio_conv) %>%
+                                       vol = vol, bio_conv = bio_conv) |> 
   agg_data(groups = c("pred", "agecl", "time", "prey"), fun = sum)
 
 # Calculate spatial overlap
@@ -158,26 +158,36 @@ rec_weight <- prm_to_df(prm_biol = prm_biol, fgs = fgs,
                         group = get_age_acronyms(fgs = fgs), 
                         parameter = c("KWRR", "KWSR", "AgeClassSize"))
 
-pd <- load_init_weight(init = init, fgs = fgs, bboxes = bboxes) %>%
-  left_join(rec_weight) %>%
-  split(.$species)
-
-# Calculate weight difference from one ageclass to the next!
-for (i in seq_along(pd)) {
-  pd[[i]]$wdiff <- c((pd[[i]]$rn[1] + pd[[i]]$sn[1]) - (pd[[i]]$kwrr[1] + pd[[i]]$kwsr[1]),
-                     diff(pd[[i]]$rn + pd[[i]]$sn))
-}
-pd <- do.call(rbind, pd)
-pd$growth_req <- pd$wdiff / (365 * pd$ageclasssize)
+pd <- load_init_weight(init = init, fgs = fgs, bboxes = bboxes) |> 
+  dplyr::left_join(rec_weight) |> 
+  
+  # Ensure structural order before running any sequential math
+  dplyr::arrange(species, agecl) |> 
+  dplyr::mutate(
+    # Total combined weight for the current row
+    total_w = rn + sn,
+    # Base weight threshold for the first age class
+    base_w = kwrr + kwsr,
+    # Calculate weight difference across age classes
+    wdiff = dplyr::if_else(
+      agecl == min(agecl),
+      total_w - base_w,                # Ageclass 1: Compare to baseline constants
+      total_w - dplyr::lag(total_w)    # Ageclass 2+: Subtract previous ageclass total weight
+    ),
+    # Calculate growth requirements using your dynamic wdiff
+    growth_req = wdiff / (365 * ageclasssize),
+    
+    .by = species
+  )
 if (any(pd$growth_req < 0)) {
   warning("Required growth negative for some groups. Please check your initial conditions files.")
 }
 
-gr_req <- pd %>%
+gr_req <- pd |> 
   select(species, agecl, growth_req)
 
-gr_rel_init <- growth_age %>%
-  left_join(gr_req) %>%
+gr_rel_init <- growth_age |> 
+  left_join(gr_req) |> 
   mutate(gr_rel = (atoutput - growth_req) / growth_req)
 
 # Aggregate volume vertically.
diff --git a/vignettes/package-demo.Rmd b/vignettes/package-demo.Rmd
index 21e5f3d..bf70f65 100644
--- a/vignettes/package-demo.Rmd
+++ b/vignettes/package-demo.Rmd
@@ -198,15 +198,15 @@ liking afterwards.
 
 ```{r}
 # Aggregate spatial biomass!
-biomass <- df_bio_spatial %>%
+biomass <- df_bio_spatial |> 
   agg_data(groups = c("species", "time"), fun = sum)
 
 plot_line(biomass, ncol = 3)
 plot_line(biomass, col = "species", ncol = 3)
 
 # Aggregate spatial biomass for fully age structured groups!
-biomass_age <- df_bio_spatial %>%
-  filter(agecl > 2) %>%
+biomass_age <- df_bio_spatial |> 
+  filter(agecl > 2) |> 
   agg_data(groups = c("species", "agecl", "time"), fun = sum)
 
 plot_line(biomass_age, col = "agecl")