diff --git a/DESCRIPTION b/DESCRIPTION
index b8922993..b22745aa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -24,6 +24,7 @@ Imports:
     dplyr,
     flextable,
     ggplot2,
+    glue,
     naniar,
     nmfspalette,
     stringr,
diff --git a/R/plot_biomass.R b/R/plot_biomass.R
index a3f052eb..5c216870 100644
--- a/R/plot_biomass.R
+++ b/R/plot_biomass.R
@@ -17,7 +17,6 @@
 #' @export
 #'
 plot_total_biomass <- function(dat,
-                                  model = "standard",
                                   show_warnings = FALSE,
                                   units = NULL,
                                   scaled = FALSE,
@@ -40,81 +39,79 @@ plot_total_biomass <- function(dat,
     bu <- "metric tons"
   }
 
-  if(model == "standard"){
-    output <- read.csv(dat)
-    totb <- output |>
-      dplyr::filter(label == "biomass",
-                    module_name == "DERIVED_QUANTITIES" | module_name == "t.series") |> # SS3 and BAM target module names
-      dplyr::mutate(estimate = as.numeric(estimate),
-                    year = as.numeric(year))
-    if (is.null(end_year)){
-      endyr <- max(totb$year)
+  output <- dat
+  totb <- output |>
+    dplyr::filter(label == "biomass",
+                  module_name == "DERIVED_QUANTITIES" | module_name == "t.series") |> # SS3 and BAM target module names
+    dplyr::mutate(estimate = as.numeric(estimate),
+                  year = as.numeric(year))
+  if (is.null(end_year)){
+    endyr <- max(totb$year)
+  }
+  # Select value for reference line and label
+  # update the target option later
+  if (any(grepl("target", output$label))) {
+    ref_line_val <- as.numeric(output[grep("(?=.*biomass)(?=.*target)", output$label, perl = TRUE),]$estimate)
+    ref_line_label <- "target"
+    if (scaled) {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[target])) # this might need to change
+    } else {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[target]))
     }
-    # Select value for reference line and label
-    # update the target option later
-    if (any(grepl("target", output$label))) {
-      ref_line_val <- as.numeric(output[grep("(?=.*biomass)(?=.*target)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "target"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[target])) # this might need to change
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[target]))
-      }
-    } else if (ref_line == "MSY" | ref_line == "msy") {
-      ref_line_val <- as.numeric(output[grep("(^biomass_msy)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "MSY"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[ref_line]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[MSY]))
-      }
-    } else if (ref_line == "unfished") {
-      ref_line_val <- as.numeric(output[grep("(^biomass_unfished)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "unfished"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[unfished]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[unfished]))
-      }
+  } else if (ref_line == "MSY" | ref_line == "msy") {
+    ref_line_val <- as.numeric(output[grep("(^biomass_msy)", output$label, perl = TRUE),]$estimate)
+    ref_line_label <- "MSY"
+    if (scaled) {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[ref_line]))
+    } else {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[MSY]))
     }
-    # Choose number of breaks for x-axis
-    x_n_breaks <- round(length(subset(totb, year<=endyr)$year)/10)
-    if (x_n_breaks <= 5) {
-      x_n_breaks <- round(length(subset(totb, year<=endyr)$year)/5)
+  } else if (ref_line == "unfished") {
+    ref_line_val <- as.numeric(output[grep("(^biomass_unfished)", output$label, perl = TRUE),]$estimate)
+    ref_line_label <- "unfished"
+    if (scaled) {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(B[unfished]))
+    } else {
+      ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(B[unfished]))
     }
-    if (relative) {
-      # plot relative TOTB
-      plt <- ggplot2::ggplot(data = subset(totb, year<=endyr)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/ref_line_val), linewidth = 1) +
-        # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/ref_line_val - stddev/ref_line_val), ymax = (value/ref_line_val + stddev/ref_line_val)), colour = "grey", alpha = 0.3) +
-        ggplot2::geom_hline(yintercept = ref_line_val/ref_line_val, linetype = 2) +
+  }
+  # Choose number of breaks for x-axis
+  x_n_breaks <- round(length(subset(totb, year<=endyr)$year)/10)
+  if (x_n_breaks <= 5) {
+    x_n_breaks <- round(length(subset(totb, year<=endyr)$year)/5)
+  }
+  if (relative) {
+    # plot relative TOTB
+    plt <- ggplot2::ggplot(data = subset(totb, year<=endyr)) +
+      ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/ref_line_val), linewidth = 1) +
+      # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/ref_line_val - stddev/ref_line_val), ymax = (value/ref_line_val + stddev/ref_line_val)), colour = "grey", alpha = 0.3) +
+      ggplot2::geom_hline(yintercept = ref_line_val/ref_line_val, linetype = 2) +
+      ggplot2::labs(x = "Year",
+                    y = paste("Biomass (", bu, ")", sep = "")) +
+      ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
+                                  guide = ggplot2::guide_axis(minor.ticks = TRUE))
+  } else {
+    if(scaled){
+      plt <- ggplot2::ggplot(data = totb) +
+        ggplot2::geom_line(ggplot2::aes(x = year, y = estimate), linewidth = 1) +
+        ggplot2::geom_hline(yintercept = ref_line_val, linetype = 2) +
         ggplot2::labs(x = "Year",
                       y = paste("Biomass (", bu, ")", sep = "")) +
         ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
                                     guide = ggplot2::guide_axis(minor.ticks = TRUE))
+      plt <- plt + ann_add
     } else {
-      if(scaled){
-        plt <- ggplot2::ggplot(data = totb) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = estimate), linewidth = 1) +
-          ggplot2::geom_hline(yintercept = ref_line_val, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Biomass (", bu, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-        plt <- plt + ann_add
-      } else {
-        plt <- ggplot2::ggplot(data = totb) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/1000), linewidth = 1) +
-          # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/1000 - stddev/1000), ymax = (value/1000 + stddev/1000)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val/1000, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Biomass (", bu, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-        plt <- plt + ann_add
-      }
+      plt <- ggplot2::ggplot(data = totb) +
+        ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/1000), linewidth = 1) +
+        # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/1000 - stddev/1000), ymax = (value/1000 + stddev/1000)), colour = "grey", alpha = 0.3) +
+        ggplot2::geom_hline(yintercept = ref_line_val/1000, linetype = 2) +
+        ggplot2::labs(x = "Year",
+                      y = paste("Biomass (", bu, ")", sep = "")) +
+        ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
+                                    guide = ggplot2::guide_axis(minor.ticks = TRUE))
+      plt <- plt + ann_add
     }
-    plt_fin <- add_theme(plt)
   }
+  plt_fin <- add_theme(plt)
   return(plt_fin)
 }
diff --git a/R/plot_landings.R b/R/plot_landings.R
index c3b867b1..225f6ccf 100644
--- a/R/plot_landings.R
+++ b/R/plot_landings.R
@@ -9,14 +9,9 @@
 #' @export
 #'
 plot_landings <- function(dat,
-                          model = "standard",
                           units = NULL){
-  # check to make sure file works with fxn
-  if (grepl(".sso|.rdat", dat)) {
-    stop("File type not compatible with function. Please use standard csv format of output files. An example can be found at https://github.com/nmfs-ost/satf")
-  }
   # read standard data file and extract target quantity
-  land <- utils::read.csv(dat) |>
+  land <- dat |>
     dplyr::filter(module_name == "t.series" | module_name == "CATCH", # t.series is associated with a conversion from BAM output and CATCH with SS3 converted output
                   grepl("landings", label) | label == "obs") |>
     dplyr::mutate(estimate = as.numeric(estimate),
diff --git a/R/plot_recruitment.R b/R/plot_recruitment.R
index 8e7fb9d6..26daba26 100644
--- a/R/plot_recruitment.R
+++ b/R/plot_recruitment.R
@@ -1,7 +1,6 @@
 #' Plot Recruitment
 #'
-#' @param dat .csv file and path where the data is located of standard stock assessment output
-#' @param model the model in which the output came from (standard, SS3, BAM, ect)
+#' @param dat A data frame returned from `asar::convert_output()`.
 #' @param params Print/export the parameters of the stock recruitment function?
 #' @param params_only Only export the stock recruitment function or both the parameters and the plot(s)?
 #' @param units If units are not available in the output file, in metric tons,
@@ -19,9 +18,7 @@
 #' reference line, stock recruitment curve, and other related figures.
 #' @export
 #'
-
 plot_recruitment <- function(dat,
-                             model = "standard",
                              params = FALSE,
                              params_only = FALSE,
                              units = c(sb = "metric tons", recruitment = "metric tons"),
@@ -46,70 +43,10 @@ plot_recruitment <- function(dat,
     sbu <- "metric tons"
   }
 
-  if (model == "standard"){
-    output <- utils::read.csv(dat)
-    if (scaled) {
-      rec <- output |>
-        dplyr::filter(label == "recruitment",
-                      module_name == "TIME_SERIES" | module_name == "t.series",
-                      !is.na(year),
-                      is.na(fleet) | length(unique(fleet)) <= 1,
-                      is.na(sex) | length(unique(sex)) <= 1,
-                      is.na(area) | length(unique(area)) <= 1,
-                      is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
-                      year != "S/Rcurve" | year != "Init" | year != "Virg"
-        ) |> # SS3 and BAM target module names
-        dplyr::mutate(estimate = as.numeric(estimate),
-                      year = as.numeric(year)) |>
-        dplyr::rename(recruitment = estimate) |>
-        dplyr::select(-c(module_name, label))
-      sb <- output |>
-        dplyr::filter(label == "spawning_biomass",
-                      module_name == "TIME_SERIES" | module_name == "t.series",
-                      !is.na(year),
-                      is.na(fleet) | length(unique(fleet)) <= 1,
-                      is.na(sex) | length(unique(sex)) <= 1,
-                      is.na(area) | length(unique(area)) <= 1,
-                      is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
-                      year != "S/Rcurve" | year != "Init" | year != "Virg"
-        ) |> # SS3 and BAM target module names
-        dplyr::mutate(estimate = as.numeric(estimate),
-                      year = as.numeric(year)) |>
-        dplyr::rename(spawning_biomass = estimate) |>
-        dplyr::select(-c(module_name, label))
-    } else {
-      rec <- output |>
-        dplyr::filter(label == "recruitment",
-                      module_name == "TIME_SERIES" | module_name == "t.series",
-                      !is.na(year),
-                      is.na(fleet) | length(unique(fleet)) <= 1,
-                      is.na(sex) | length(unique(sex)) <= 1,
-                      is.na(area) | length(unique(area)) <= 1,
-                      is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
-                      year != "S/Rcurve" | year != "Init" | year != "Virg"
-        ) |> # SS3 and BAM target module names
-        dplyr::mutate(estimate = as.numeric(estimate),
-                      year = as.numeric(year)) |>
-        dplyr::rename(recruitment = estimate) |>
-        dplyr::select(-c(module_name, label))
-      sb <- output |>
-        dplyr::filter(label == "spawning_biomass",
-                      module_name == "TIME_SERIES" | module_name == "t.series",
-                      !is.na(year),
-                      is.na(fleet) | length(unique(fleet)) <= 1,
-                      is.na(sex) | length(unique(sex)) <= 1,
-                      is.na(area) | length(unique(area)) <= 1,
-                      is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
-                      year != "S/Rcurve" | year != "Init" | year != "Virg"
-        ) |> # SS3 and BAM target module names
-        dplyr::mutate(estimate = as.numeric(estimate),
-                      year = as.numeric(year)) |>
-        dplyr::rename(spawning_biomass = estimate) |>
-        dplyr::select(-c(module_name, label))
-    }
-
-    rec_devs <- output |>
-      dplyr::filter(label == "recruitment_deviations" | label == "log_recruitment_deviations",
+  output <- dat
+  if (scaled) {
+    rec <- output |>
+      dplyr::filter(label == "recruitment",
                     module_name == "TIME_SERIES" | module_name == "t.series",
                     !is.na(year),
                     is.na(fleet) | length(unique(fleet)) <= 1,
@@ -120,203 +57,111 @@ plot_recruitment <- function(dat,
       ) |> # SS3 and BAM target module names
       dplyr::mutate(estimate = as.numeric(estimate),
                     year = as.numeric(year)) |>
-      dplyr::rename(recruitment_deviations = estimate) |>
+      dplyr::rename(recruitment = estimate) |>
       dplyr::select(-c(module_name, label))
-    if(return == "recruitment_deviations" & nrow(rec_devs)==0){
-      stop("No recruitment deviations found in data.")
-    }
-
-    if (is.null(end_year)){
-      endyr <- max(rec$year)
-    } else {
-      endyr <- end_year
-    }
-    stryr <- min(rec$year)
-
-    # merge DF
-    sr <- dplyr::full_join(sb, rec)
-
-    # Choose number of breaks for x-axis
-    x_n_breaks <- round(length(subset(sr, year<=endyr)$year)/10)
-    if (x_n_breaks <= 5) {
-      x_n_breaks <- round(length(subset(sr, year<=endyr)$year)/5)
-    }
-
-    if (return == "stock_recruitment"){
-      plt <- ggplot2::ggplot(data = sr) +
-        ggplot2::geom_line(ggplot2::aes(x = spawning_biomass, y = recruitment/1000), linewidth = 1) +
-        ggplot2::labs(x = paste("Spawning Biomass (", sbu, ")", sep = ""),
-                      y = paste("Recruitment (", ru, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none")
-      # sr_plt <- add_theme(sr_plt)
-    } else if (return == "recruitment") {
-      plt <- ggplot2::ggplot(data = sr) +
-        ggplot2::geom_point(ggplot2::aes(x = year, y = recruitment)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = recruitment), linewidth = 1) +
-        ggplot2::labs(x = "Year",
-                      y = paste("Recruitment (", ru, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none") +
-        ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                    guide = ggplot2::guide_axis(minor.ticks = TRUE))
-    } else if (return == "recruitment_deviations") {
-      plt <- ggplot2::ggplot(data = rec_devs) +
-        # ggplot2::geom_point(ggplot2::aes(x = year, y = log_rec_dev), shape = 1, size = 2.5) +
-        ggplot2::geom_pointrange(ggplot2::aes(x = year, y = estimate, ymax = Value, ymin = 0),  fatten = 1, size = 2, shape = 1) +
-        ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
-        ggplot2::labs(x = "Year",
-                      y = "Recruitment Deviations")
-    }
-    plt_fin <- add_theme(plt)
-  } else if(model == "SS3") {
-    # Read rep file or rename
-    if(grepl(".sso", dat)){
-      get_ncol <- function(file, skip = 0) {
-        nummax <- max(utils::count.fields(file,
-                                          skip = skip, quote = "",
-                                          comment.char = ""
-        )) + 1
-        return(nummax)
-      }
-
-      output <- utils::read.table(
-        file = dat, col.names = 1:get_ncol(dat), fill = TRUE, quote = "",
-        nrows = -1, comment.char = "", blank.lines.skip = FALSE
-      )
-    } else {
-      output <- dat
-    }
-    # extract recruitment
-    sr_info <- SS3_extract_df(output, "SPAWN_RECRUIT")
-    sr <- sr_info[-c(1:14),]
-    colnames(sr) <- tolower(sr_info[13,])
-    sr <- sr  |>
-      dplyr::mutate(year = yr,
-                    spawnbio = as.numeric(spawnbio),
-                    exp_recr = as.numeric(exp_recr),
-                    with_regime = as.numeric(with_regime),
-                    bias_adjusted = as.numeric(bias_adjusted),
-                    pred_recr = as.numeric(pred_recr),
-                    # dev = as.numeric(dev),
-                    # biasadjuster = as.numeric(biasadjuster)
-                    )
-
-    # Store virgin and initial recruitment for reference in plotting
-    time_series <- SS3_extract_df(output, "TIME_SERIES")
-    colnames(time_series) <- tolower(time_series[2,])
-    time_series <- time_series[-c(1:2),]
-    time_series[, 5:ncol(time_series)] <- lapply(5:ncol(time_series), function(x) suppressWarnings(as.numeric(time_series[[x]])))
-
-    # Adapted from r4ss
-    B0 <- sum(time_series[["spawnbio"]][time_series[["era"]] == "VIRG"], na.rm = TRUE)
-    B1 <- sum(time_series[["spawnbio"]][time_series[["era"]] == "INIT"], na.rm = TRUE)
-    R0 <- sum(time_series[["recruit_0"]][time_series[["era"]] == "VIRG"], na.rm = TRUE)
-    R1 <- sum(time_series[["recruit_0"]][time_series[["era"]] == "INIT"], na.rm = TRUE)
-
-    # Units - change this statement to fit new convention
-    sb_units = spawning_biomass_units
-    rec_units = recruitment_units
-    message("Default units for both SB and R are metric tons.")
-
-    # remove inital and virg recruitment
-    sr <- sr |>
-      dplyr::filter(year!="Init" & year!="Virg" & era != "Forecast") |>
-      dplyr::mutate(year = as.numeric(year))
-
-    if (return == "stock_recruitment") {
-      # need to rescale to 1000s rather than xe.
-      plt <- ggplot2::ggplot(data = sr) +
-        ggplot2::geom_line(ggplot2::aes(x = spawnbio/1000, y = exp_recr/1000), linewidth = 1) + # exp. R
-        # add exp R after bias adjustment (dotted line)
-        ggplot2::geom_point(ggplot2::aes(x = spawnbio/1000, y = pred_recr/1000, color = year)) + # change colors
-        # ggplot2::geom_text() +
-        ggplot2::labs(x = paste("Spawning Biomass (", sb_units, ")", sep = ""),
-             y = paste("Recruitment (", rec_units, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none")
-    } else if (return == "recruitment") {
-      plt <- ggplot2::ggplot(data = sr) +
-        ggplot2::geom_point(ggplot2::aes(x = year, y = pred_recr)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = pred_recr), linewidth = 1) +
-        ggplot2::labs(x = "Year",
-                      y = paste("Recruitment (", rec_units, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none")
-    } else if (return == "recruitment_deviations") {
-      # recruitment deviations
-      params <- SS3_extract_df(output, "PARAMETERS")
-      colnames(params) <- params[2,]
-      params2 <- params[-c(1:2),] |>
-        dplyr::filter(grepl('RecrDev', Label)) |>
-        dplyr::select(Label, Value) |>
-        tidyr::separate_wider_delim(Label, delim = "_", names = c("Era", "Recr", "Year"))
-      params_proj <-  params[-c(1:2),] |>
-        dplyr::filter(grepl('ForeRecr', Label)) |>
-        dplyr::select(Label, Value) |>
-        tidyr::separate_wider_delim(Label, delim = "_", names = c("Recr", "Year")) |>
-        dplyr::mutate(Era = "Fore")
-
-      params_fin <- rbind(params2, params_proj)
-
-      sr2 <- dplyr::left_join(sr, params_fin, by = c("yr" = "Year")) |>
-        dplyr::mutate(Value = as.numeric(Value))
-
-      plt <- ggplot2::ggplot(data = sr2) +
-        # ggplot2::geom_point(ggplot2::aes(x = year, y = log_rec_dev), shape = 1, size = 2.5) +
-        ggplot2::geom_pointrange(ggplot2::aes(x = year, y = exp_recr, ymax = Value, ymin = 0),  fatten = 1, size = 2, shape = 1) +
-        ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
-        ggplot2::labs(x = "Year",
-                      y = "logR Deviations")
-    }
-    plt_fin <- add_theme(plt)
-  } # close SS3 if statement
+    sb <- output |>
+      dplyr::filter(label == "spawning_biomass",
+                    module_name == "TIME_SERIES" | module_name == "t.series",
+                    !is.na(year),
+                    is.na(fleet) | length(unique(fleet)) <= 1,
+                    is.na(sex) | length(unique(sex)) <= 1,
+                    is.na(area) | length(unique(area)) <= 1,
+                    is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
+                    year != "S/Rcurve" | year != "Init" | year != "Virg"
+      ) |> # SS3 and BAM target module names
+      dplyr::mutate(estimate = as.numeric(estimate),
+                    year = as.numeric(year)) |>
+      dplyr::rename(spawning_biomass = estimate) |>
+      dplyr::select(-c(module_name, label))
+  } else {
+    rec <- output |>
+      dplyr::filter(label == "recruitment",
+                    module_name == "TIME_SERIES" | module_name == "t.series",
+                    !is.na(year),
+                    is.na(fleet) | length(unique(fleet)) <= 1,
+                    is.na(sex) | length(unique(sex)) <= 1,
+                    is.na(area) | length(unique(area)) <= 1,
+                    is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
+                    year != "S/Rcurve" | year != "Init" | year != "Virg"
+      ) |> # SS3 and BAM target module names
+      dplyr::mutate(estimate = as.numeric(estimate),
+                    year = as.numeric(year)) |>
+      dplyr::rename(recruitment = estimate) |>
+      dplyr::select(-c(module_name, label))
+    sb <- output |>
+      dplyr::filter(label == "spawning_biomass",
+                    module_name == "TIME_SERIES" | module_name == "t.series",
+                    !is.na(year),
+                    is.na(fleet) | length(unique(fleet)) <= 1,
+                    is.na(sex) | length(unique(sex)) <= 1,
+                    is.na(area) | length(unique(area)) <= 1,
+                    is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
+                    year != "S/Rcurve" | year != "Init" | year != "Virg"
+      ) |> # SS3 and BAM target module names
+      dplyr::mutate(estimate = as.numeric(estimate),
+                    year = as.numeric(year)) |>
+      dplyr::rename(spawning_biomass = estimate) |>
+      dplyr::select(-c(module_name, label))
+  }
 
-  if(model == "BAM"){
-    # read in output
-    output <- dget(dat)
+  rec_devs <- output |>
+    dplyr::filter(label == "recruitment_deviations" | label == "log_recruitment_deviations",
+                  module_name == "TIME_SERIES" | module_name == "t.series",
+                  !is.na(year),
+                  is.na(fleet) | length(unique(fleet)) <= 1,
+                  is.na(sex) | length(unique(sex)) <= 1,
+                  is.na(area) | length(unique(area)) <= 1,
+                  is.na(growth_pattern) | length(unique(growth_pattern)) <= 1,
+                  year != "S/Rcurve" | year != "Init" | year != "Virg"
+    ) |> # SS3 and BAM target module names
+    dplyr::mutate(estimate = as.numeric(estimate),
+                  year = as.numeric(year)) |>
+    dplyr::rename(recruitment_deviations = estimate) |>
+    dplyr::select(-c(module_name, label))
+  if(return == "recruitment_deviations" & nrow(rec_devs)==0){
+    stop("No recruitment deviations found in data.")
+  }
 
-    R_virg <- output$parms$R.virgin.bc
-    R0 <- output$parms$R0
-    SSBmsy <- output$parms$SSBmsy
-    SSB0 <- output$parms$SSB0
+  if (is.null(end_year)){
+    endyr <- max(rec$year)
+  } else {
+    endyr <- end_year
+  }
+  stryr <- min(rec$year)
 
-    # SSB time series
-    sr <- data.frame(year = output$t.series$year,
-                     spawn_bio = output$t.series$SSB,
-                     pred_recr = output$t.series$recruits,
-                     # exp_recr = output$t.series,
-                     log_rec_dev = output$t.series$logR.dev)
+  # merge DF
+  sr <- dplyr::full_join(sb, rec)
 
-    # Check if units were declared
-    sb_units = spawning_biomass_units
-    rec_units = recruitment_units
-    message("Default units for both SB and R are metric tons.")
+  # Choose number of breaks for x-axis
+  x_n_breaks <- round(length(subset(sr, year<=endyr)$year)/10)
+  if (x_n_breaks <= 5) {
+    x_n_breaks <- round(length(subset(sr, year<=endyr)$year)/5)
+  }
 
-    if (return == "stock_recruitment") {
-      # plot stock recruitment in 1000
-      plt <- ggplot2::ggplot(data = sr) +
-        # ggplot2::geom_line(ggplot2::aes(x = spawn_bio, y = exp_recr), linewidth = 1) + # exp. R
-        # add exp R after bias adjustment (dotted line)
-        ggplot2::geom_point(ggplot2::aes(x = spawn_bio, y = pred_recr, color = year)) + # change colors
-        # ggplot2::geom_text() +
-        ggplot2::labs(x = paste("Spawning Biomass (", sb_units, ")", sep = ""),
-                      y = paste("Recruitment (", rec_units, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none")
-    } else if (return == "recruitment") {
-      plt <- ggplot2::ggplot(data = sr) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = pred_recr), linewidth = 1) + # exp. R
-        ggplot2::geom_point(ggplot2::aes(x = year, y = pred_recr)) + # change colors
-        # ggplot2::geom_text() +
-        ggplot2::labs(x = "Year",
-                      y = paste("Recruitment (", rec_units, ")", sep = "")) +
-        ggplot2::theme(legend.position = "none")
-    } else if (return == "recruitment_deviations") {
-      plt <- ggplot2::ggplot(data = sr) +
-        # ggplot2::geom_point(ggplot2::aes(x = year, y = log_rec_dev), shape = 1, size = 2.5) +
-        ggplot2::geom_pointrange(ggplot2::aes(x = year, y = log_rec_dev, ymax = log_rec_dev, ymin = 0),  fatten = 1, size = 2, shape = 1) +
-        ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
-        ggplot2::labs(x = "Year",
-                      y = "logR Deviations")
-    }
-    plt_fin <- add_theme(plt)
-  } # close BAM if statement
+  if (return == "stock_recruitment"){
+    plt <- ggplot2::ggplot(data = sr) +
+      ggplot2::geom_line(ggplot2::aes(x = spawning_biomass, y = recruitment/1000), linewidth = 1) +
+      ggplot2::labs(x = paste("Spawning Biomass (", sbu, ")", sep = ""),
+                    y = paste("Recruitment (", ru, ")", sep = "")) +
+      ggplot2::theme(legend.position = "none")
+    # sr_plt <- add_theme(sr_plt)
+  } else if (return == "recruitment") {
+    plt <- ggplot2::ggplot(data = sr) +
+      ggplot2::geom_point(ggplot2::aes(x = year, y = recruitment)) +
+      ggplot2::geom_line(ggplot2::aes(x = year, y = recruitment), linewidth = 1) +
+      ggplot2::labs(x = "Year",
+                    y = paste("Recruitment (", ru, ")", sep = "")) +
+      ggplot2::theme(legend.position = "none") +
+      ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
+                                  guide = ggplot2::guide_axis(minor.ticks = TRUE))
+  } else if (return == "recruitment_deviations") {
+    plt <- ggplot2::ggplot(data = rec_devs) +
+      # ggplot2::geom_point(ggplot2::aes(x = year, y = log_rec_dev), shape = 1, size = 2.5) +
+      ggplot2::geom_pointrange(ggplot2::aes(x = year, y = estimate, ymax = Value, ymin = 0),  fatten = 1, size = 2, shape = 1) +
+      ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
+      ggplot2::labs(x = "Year",
+                    y = "Recruitment Deviations")
+  }
+  plt_fin <- add_theme(plt)
   return(plt_fin)
 }
-
diff --git a/R/plot_spawning_biomass.R b/R/plot_spawning_biomass.R
index 8cd0ef95..caa324ed 100644
--- a/R/plot_spawning_biomass.R
+++ b/R/plot_spawning_biomass.R
@@ -1,363 +1,130 @@
-#' Plot Spawning Biomass
+#' Plot spawning biomass (SB)
+#'
+#' Plot spawning biomass with a reference line as a dashed line. The figure can
+#' also be made relative to this reference line rather than in absolute units.
 #'
 #' @inheritParams plot_recruitment
-#' @param show_warnings Option to suppress warnings
-#' @param biomass_units units for biomass
-#' @param ref_line choose with reference point to plot a reference line and use
-#' in relative sb calculations
-#' @param end_year input the end year of the stock assessment data (not including
-#' projections). This parameter will be deprecated once the output converter is fully developed.
-#' @param relative Plot relative spawning biomass. Ref line indicates which reference point to use
+#' @param ref_line A string specifying the type of reference you want to
+#'   compare spawning biomass to. The default is `"target"`, which looks for
+#'   `"spawning_biomass_target"` in the `"label"` column of `dat`. The actual
+#'   searching in `dat` is case agnostic and will work with either upper- or
+#'   lower-case letters but you must use one of the options specified in the
+#'   default list to ensure that the label on the figure looks correct
+#'   regardless of how it is specified in `dat`.
+#' @param relative A logical value specifying if the resulting figures should
+#'   be relative spawning biomass. The default is `FALSE`. `ref_line` indicates
+#'   which reference point to use.
 #'
-#' @return Plot spawning biomass from a stock assessment model as found in a NOAA
-#' stock assessment report. Units of spawning biomass can either be manually added
-#' or will be extracted from the provided file if possible. In later releases, model will not
+#' @return
+#' Plot spawning biomass from the results of an assessment model translated to
+#' the standard output. The {ggplot2} object is returned for further
+#' modifications if needed.
 #' @export
 #'
-plot_spawning_biomass <- function(dat,
-                         model = "standard",
-                         show_warnings = FALSE,
-                         units = NULL,
-                         # biomass_units = NULL,
-                         spawning_biomass_units = NULL,
-                         scaled = FALSE,
-                         scale_amount = 1000,
-                         ref_line = c("target", "MSY", "msy", "unfished"),
-                         end_year = NULL,
-                         relative = FALSE
-                         ){
+plot_spawning_biomass <- function(
+  dat,
+  unit_label = "metric ton",
+  scale_amount = 1,
+  ref_line = c("target", "unfished"),
+  end_year = NULL,
+  relative = FALSE
+) {
+  ref_line <- match.arg(ref_line)
+  # TODO: Fix the unit label if scaling. Maybe this is up to the user to do if
+  #       they want something scaled then they have to supply a better unit name
+  #       or we create a helper function to do this.
+  spawning_biomass_label <- ifelse(
+    relative,
+    yes = "Relative spawning biomass",
+    no = glue::glue("Spawning biomass ({unit_label})")
+  )
 
-  if(length(ref_line)>1){
-    ref_line = "target"
-  } else {
-    ref_line <- match.arg(ref_line, several.ok = FALSE)
+  output <- dat
+  # Determine the end year
+  all_years <- output[["year"]][grepl("^[0-9\\.]+$", output[["year"]])]
+  if (is.null(end_year)) {
+    end_year <- as.numeric(max(all_years, na.rm = TRUE))
   }
+  stopifnot(any(end_year >= all_years))
 
-  # check units
-  # spawning biomass
-  if(!is.null(spawning_biomass_units)){
-    sbu <- spawning_biomass_units
-  } else {
-    sbu <- "metric tons"
+  # Select value for reference line and label
+  ref_line_val <- as.numeric(output[
+    grep(
+      pattern = glue::glue("^spawning_biomass.*{tolower(ref_line)}"),
+      x = output[["label"]]
+    ),
+    "estimate"
+  ])
+  if (length(ref_line_val) == 0) {
+    stop(glue::glue(
+      "The resulting reference value of `spawning_biomass_{ref_line}` was
+      not found in `dat[[\"label\"]]`."
+    ))
   }
+  sb <- output |>
+    dplyr::filter(
+      label == "spawning_biomass",
+      module_name %in% c("DERIVED_QUANTITIES", "t.series"),
+      year <= end_year
+    ) |>
+    dplyr::mutate(
+      estimate = as.numeric(estimate),
+      year = as.numeric(year),
+      estimate_y = estimate / ifelse(relative, ref_line_val, scale_amount),
+      # TODO: Determine what unit uncertainty is in, following is for normal sd
+      # TODO: This fails for Bayesian estimation
+      estimate_lower = (estimate - 1.96 * uncertainty) /
+        ifelse(relative, ref_line_val, scale_amount),
+      estimate_upper = (estimate + 1.96 * uncertainty) /
+        ifelse(relative, ref_line_val, scale_amount)
+    )
 
-  if(model == "standard"){
-    output <- utils::read.csv(dat)
-    sb <- output |>
-      dplyr::filter(label == "spawning_biomass",
-                    module_name == "DERIVED_QUANTITIES" | module_name == "t.series") |> # SS3 and BAM target module names
-      dplyr::mutate(estimate = as.numeric(estimate),
-                    year = as.numeric(year))
-    if (is.null(end_year)){
-      endyr <- max(sb$year)
-    } else {
-      endyr <- end_year
-    }
-    # Select value for reference line and label
-    if (any(grepl("target", output$label))) {
-      ref_line_val <- as.numeric(output[grep("(?=.*spawning_biomass)(?=.*target)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "target"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(SB[target])) # this might need to change
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(SB[target]))
-      }
-    } else if (ref_line == "MSY" | ref_line == "msy") {
-      ref_line_val <- as.numeric(output[grep("(?=.*spawning_biomass)(?=.*msy)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "MSY"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(SB[MSY]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(SB[MSY]))
-      }
-    } else if (ref_line == "unfished") {
-      ref_line_val <- as.numeric(output[grep("(?=.*spawning_biomass)(?=.*unfished)", output$label, perl = TRUE),]$estimate)
-      ref_line_label <- "unfished"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val, label = bquote(SB[unfished]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr + 0.05, y=ref_line_val/1000, label = bquote(SB[unfished]))
-      }
-    }
-    # Choose number of breaks for x-axis
-    x_n_breaks <- round(length(subset(sb, year<=endyr)$year)/10)
-    if (x_n_breaks <= 5) {
-      x_n_breaks <- round(length(subset(sb, year<=endyr)$year)/5)
-    }
-    if (relative) {
-    # plot relative SB
-      plt <- ggplot2::ggplot(data = subset(sb, year<=endyr)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/ref_line_val), linewidth = 1) +
-        # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/ref_line_val - stddev/ref_line_val), ymax = (value/ref_line_val + stddev/ref_line_val)), colour = "grey", alpha = 0.3) +
-        ggplot2::geom_hline(yintercept = ref_line_val/ref_line_val, linetype = 2) +
-        ggplot2::labs(x = "Year",
-                      y = paste("Spawning Biomass (", sbu, ")", sep = "")) +
-        ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                    guide = ggplot2::guide_axis(minor.ticks = TRUE))
-    } else {
-      if(scaled){
-        plt <- ggplot2::ggplot(data = sb) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = estimate), linewidth = 1) +
-          ggplot2::geom_hline(yintercept = ref_line_val, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sbu, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-        plt <- plt + ann_add
-      } else {
-        plt <- ggplot2::ggplot(data = sb) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = estimate/1000), linewidth = 1) +
-          # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/1000 - stddev/1000), ymax = (value/1000 + stddev/1000)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val/1000, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sbu, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-        plt <- plt + ann_add
-      }
-    }
-    plt_fin <- add_theme(plt)
-  } else if (model == "SS3") {
-    # load SS3 data file
-    # check if dat parameter is a report file or a mutates df
-    if(grepl(".sso", dat)){
-      get_ncol <- function(file, skip = 0) {
-        nummax <- max(utils::count.fields(file,
-                                          skip = skip, quote = "",
-                                          comment.char = ""
-        )) + 1
-        return(nummax)
-      }
-
-      output <- utils::read.table(
-        file = dat, col.names = 1:get_ncol(dat), fill = TRUE, quote = "",
-        colClasses = "character", nrows = -1, comment.char = "",
-        blank.lines.skip = FALSE
-      )
-    } else {
-      output <- dat
-    }
-
-    bio_info <- SS3_extract_df(output, "DERIVED_QUANTITIES")[-c(1:4),]
-    colnames(bio_info) <- bio_info[1,]
-    bio_info <- bio_info[-1,] |>
-      dplyr::mutate(year = stringr::str_extract(Label, "[0-9]+$"),
-                    Label = stringr::str_remove(Label, "_[0-9]+$"))
-
-    sb <- bio_info |>
-      dplyr::filter(Label == "SSB") |>
-      dplyr::mutate(Value = as.numeric(Value),
-                    year = as.numeric(year),
-                    StdDev = as.numeric(StdDev))
-    if(is.null(end_year)){
-      endyr <- max(sb$year)
-    }
-
-    # Check if units were declared
-    if(!is.null(units)){
-      if(length(units)>1){
-        sb_units <- units[1]
-        rec_units <- units[2]
-        message("Please check the units on your axes are correct. If they are flipped, change the order of names in the units argument.")
-      } else {
-        if(grepl("eggs", units)){
-          sb_units <- "1e10 eggs"
-          rec_units <- "metric tons"
-        } else if(grepl("number", units)){
-          rec_units <- "number of fish"
-          sb_units <- "metric tons"
-        } else {
-          warning("Unit type is not defined for this function. Please leave an issue at https://github.com/nmfs-ost/satf/issues")
-        }
-      }
-    } else {
-      sb_units <- "metric tons"
-      rec_units <- "metric tons"
-      message("Default units for both SB and R are metric tons.")
-    }
-
-    # Select value for reference line and label
-    if (ref_line == "target") {
-      ref_line_val <- as.numeric(bio_info[['Value']][bio_info[['Label']]=="SSB_Btgt"])
-      ref_line_label <- "target"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[Btarget]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[Btarget]))
-      }
-    } else if (ref_line == "MSY" | ref_line == "msy") {
-      ref_line_val <- as.numeric(bio_info[['Value']][bio_info[['Label']]=="SSB_MSY"])
-      ref_line_label <- "MSY"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[MSY]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[MSY]))
-      }
-    } else if (ref_line == "unfished") {
-      ref_line_val <- as.numeric(bio_info[['Value']][bio_info[['Label']]=="SSB_unfished"])
-      ref_line_label <- "unfished"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[unfished]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[unfished]))
-      }
-    } else if (ref_line == "spr") {
-      ref_line_val <- as.numeric(bio_info[['Value']][bio_info[['Label']]=="SSB_SPR"])
-      ref_line_label <- "SPR"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[SPR]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[SPR]))
-      }
-    } else {
-      ref_line_name <- paste("SSB_", ref_line, sep = "")
-      ref_line_val <- as.numeric(bio_info[['Value']][bio_info[['Label']]==ref_line_name])
-      ref_line_label <- ref_line
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[ref]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[ref]))
-      }
-    }
-    # Choose number of breaks for x-axis
-    x_n_breaks <- round(length(subset(sb, year<=endyr)$year)/10)
-    if (x_n_breaks <= 5) {
-      x_n_breaks <- round(length(subset(sb, year<=endyr)$year)/5)
-    }
-    # Plot of rel. SB
-    if (relative) {
-      plt <- ggplot2::ggplot(data = subset(sb, year <= endyr)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = Value/ref_line_val), linewidth = 1) +
-        ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (Value/ref_line_val - StdDev/ref_line_val), ymax = (Value/ref_line_val + StdDev/ref_line_val)), colour = "grey", alpha = 0.3) +
-        ggplot2::geom_hline(yintercept = ref_line_val/ref_line_val, linetype = 2) +
-        ggplot2::labs(x = "Year",
-                      y = paste("Rel. Spawning Biomass", sep = "")) +
-        ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                    guide = ggplot2::guide_axis(minor.ticks = TRUE))
-    } else {
-      if (scaled) {
-        plt <- ggplot2::ggplot(data = subset(sb, year < endyr + 1)) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = Value), linewidth = 1) +
-          ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (Value - StdDev), ymax = (Value + StdDev)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sb_units, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-      } else {
-        plt <- ggplot2::ggplot(data = subset(sb, year < endyr + 1)) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = Value/1000), linewidth = 1) +
-          ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (Value/1000 - StdDev/1000), ymax = (Value/1000 + StdDev/1000)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val/1000, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sb_units, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-      }
-      plt <- plt + ann_add
-    }
-    plt_fin <- add_theme(plt)
-  } # close SS3 if statement
-
-  if (model == "BAM"){
-    # extract spawning biomass from output
-    output <- dget(dat)
-    sb <- data.frame(sapply(output$t.series, c)) |>
-      dplyr::select(year, SSB)
-
-    # Projection years
-    proj <- data.frame(sapply(output$proj.t.series, c)) |>
-      dplyr::select(year, SSB.proj) |>
-      dplyr::rename(SSB = SSB.proj)
-
-    sb2 <- rbind(sb, proj) |>
-      dplyr::rename(value = SSB)
-
-    # max_yr <- max(unique(sb2$year))
-    if(is.null(end_year)){
-      endyr <- output$parms$endyr
-    }
-    stryr <- output$parms$styr
-
-    # Check if units were declared
-    if(!is.null(units)){
-      sb_units <- units[1]
-      rec_units <- units[2]
-      message("Please check the units on your axes are correct. If they are flipped, change the order of names in the units argument.")
-    } else {
-      sb_units <- output$info$units.ssb
-      rec_units <- output$info$units.rec
-    }
+  # Choose number of breaks for x-axis
+  x_n_breaks <- round(length(sb[["year"]]) / 10)
+  if (x_n_breaks <= 5) {
+    x_n_breaks <- round(length(sb[["year"]]) / 5)
+  }
 
-    # Pull reference pts
-    # unfished <- output$parms$SSB0
-    # msy <- output$parms$SSBmsy
-    # tgt <- output$parms$SSB.F30 # SSBF30 but calling tgt to use same plotting for all
+  plt <- ggplot2::ggplot(data = sb) +
+    ggplot2::geom_line(
+      ggplot2::aes(
+        x = year,
+        y = estimate_y
+      ),
+      linewidth = 1
+    ) +
+    ggplot2::geom_hline(
+      yintercept = ref_line_val / ifelse(relative, ref_line_val, scale_amount),
+      linetype = 2
+    ) +
+    # Only add confidence intervals for the non NA estimates
+    # which allows for no warnings if uncertainty = NA
+    ggplot2::geom_ribbon(
+      data = sb |> dplyr::filter(!is.na(estimate_lower)),
+      ggplot2::aes(
+        x = year,
+        ymin = estimate_lower,
+        ymax = estimate_upper
+      ),
+      colour = "grey",
+      alpha = 0.3,
+    ) +
+    ggplot2::labs(
+      x = "Year",
+      y = spawning_biomass_label
+    ) +
+    ggplot2::scale_x_continuous(
+      n.breaks = x_n_breaks,
+      guide = ggplot2::guide_axis(minor.ticks = TRUE)
+    ) +
+    ggplot2::annotate(
+      geom = "text",
+      x = end_year + 0.05,
+      y = ref_line_val / ifelse(relative, ref_line_val, scale_amount),
+      label = list(bquote(SB[.(ref_line)])),
+      parse = TRUE
+    )
 
-    # Select value for reference line and label
-    if (ref_line == "target") {
-      ref_line_val <- output$parms$SSB.Fproxy
-      ref_line_label <- "target"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[F30]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[F30])) # this might need to change
-      }
-    } else if (ref_line == "MSY" | ref_line == "msy") {
-      ref_line_val <- output$parms$SSBmsy
-      ref_line_label <- "MSY"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[MSY]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[MSY]))
-      }
-    } else if (ref_line == "unfished") {
-      ref_line_val <- output$parms$SSB0
-      ref_line_label <- "unfished"
-      if (scaled) {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val, label = bquote(SB[unfished]))
-      } else {
-        ann_add <- ggplot2::annotate("text", x = endyr+0.5, y=ref_line_val/1000, label = bquote(SB[unfished]))
-      }
-    }
-    # Choose number of breaks for x-axis
-    x_n_breaks <- round(length(subset(sb2, year<=endyr)$year)/10)
-    if (x_n_breaks <= 5) {
-      x_n_breaks <- round(length(subset(sb2, year<=endyr)$year)/5)
-    }
-    if (relative) {
-      # plot relative SB
-      plt <- ggplot2::ggplot(data = subset(sb2, year<=endyr)) +
-        ggplot2::geom_line(ggplot2::aes(x = year, y = value/ref_line_val), linewidth = 1) +
-        # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/ref_line_val - stddev/ref_line_val), ymax = (value/ref_line_val + stddev/ref_line_val)), colour = "grey", alpha = 0.3) +
-        ggplot2::geom_hline(yintercept = ref_line_val/ref_line_val, linetype = 2) +
-        ggplot2::labs(x = "Year",
-                      y = paste("Spawning Biomass (", sb_units, ")", sep = "")) +
-        ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                    guide = ggplot2::guide_axis(minor.ticks = TRUE))
-    } else {
-      if (scaled) {
-        plt <- ggplot2::ggplot(data = subset(sb2, year<=endyr)) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = value), linewidth = 1) +
-          # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/1000 - stddev/1000), ymax = (value/1000 + stddev/1000)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sb_units, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = xn_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-      } else {
-        plt <- ggplot2::ggplot(data = subset(sb2, year<=endyr)) +
-          ggplot2::geom_line(ggplot2::aes(x = year, y = value/1000), linewidth = 1) +
-          # ggplot2::geom_ribbon(ggplot2::aes(x = year, ymin = (value/1000 - stddev/1000), ymax = (value/1000 + stddev/1000)), colour = "grey", alpha = 0.3) +
-          ggplot2::geom_hline(yintercept = ref_line_val/1000, linetype = 2) +
-          ggplot2::labs(x = "Year",
-                        y = paste("Spawning Biomass (", sb_units, ")", sep = "")) +
-          ggplot2::scale_x_continuous(n.breaks = x_n_breaks,
-                                      guide = ggplot2::guide_axis(minor.ticks = TRUE))
-      }
-      plt <- plt + ann_add
-    }
-    plt_fin <- suppressWarnings(add_theme(plt))
-  }
+  plt_fin <- suppressWarnings(add_theme(plt))
   return(plt_fin)
 }
diff --git a/R/table_indices.R b/R/table_indices.R
index 018de1e1..14f85f62 100644
--- a/R/table_indices.R
+++ b/R/table_indices.R
@@ -6,153 +6,28 @@
 #' @export
 #'
 
-table_indices <- function(dat,
-                          model = "standard"){
-  if (model == "standard"){
-    output <- utils::read.csv(dat)
+table_indices <- function(dat) {
+  output <- dat
+  output <- output |>
+    dplyr::filter(module_name == "INDEX_2" | module_name == "t.series")
+  if (any(unique(output$module_name=="INDEX_2"))) {
     output <- output |>
-      dplyr::filter(module_name == "INDEX_2" | module_name == "t.series")
-    if (any(unique(output$module_name=="INDEX_2"))) {
-      output <- output |>
-        dplyr::filter(grepl("obs", label))
-    } else if (any(unique(output$module_name=="t.series"))) {
-      output <- output |>
-        dplyr::filter(grepl("cpue", label))
-    }
-    fleet_names <- unique(output$fleet)
-    factors <- c("year", "fleet", "fleet_name", "age", "sex", "area", "seas", "season", "time", "era", "subseas", "subseason", "platoon", "platoo","growth_pattern", "gp")
-    # re-structure df for table
-    indices <- output |>
-      dplyr::rename(!!unique(output$label) := estimate,
-                    !!unique(output$uncertainty_label) := uncertainty) |>
-      tidyr::pivot_wider(
-        id_cols = -intersect(colnames(output), factors),
-        names_from = fleet,
-        values_from = c(unique(output$label), unique(output$uncertainty_label))
-      ) # stated internal error for tidyr and asks to report - try again monday
-
-  } else if (model == "SS3"){
-    # Read report file
-    if(grepl(".sso|.txt", dat)){
-      get_ncol <- function(file, skip = 0) {
-        nummax <- max(utils::count.fields(file,
-                                          skip = skip, quote = "",
-                                          comment.char = ""
-        )) + 1
-        return(nummax)
-      }
-
-      output <- utils::read.table(
-        file = dat, col.names = 1:get_ncol(dat), fill = TRUE, quote = "",
-        colClasses = "character", nrows = -1, comment.char = "",
-        blank.lines.skip = FALSE
-      )
-    }
-    # Extract Fleet Names
-    fleet_info <- SS3_extract_df(output, "Fleet")[,1]
-    colnames(fleet_info) <- fleet_info[1,1]
-    fleet_info <- fleet_info[-1,]
-    fleet_names <- unlist(as.vector(fleet_info))
-    fleet_names <- paste("Fleet_", fleet_names, sep ="")
-    # fleet_names <- as.character(fleet_info[fleet_info$X1=="fleet_names:",][,-1])
-    # fleet_names <- gsub("_", " ", fleet_names)
-
-    # Index 1
-    index1 <- SS3_extract_df(output, "INDEX_1")
-    colnames(index1) <- index1[2,]
-    index1 <- index1[-c(1:2),]
-
-    # Index 2
-    index2 <- SS3_extract_df(output, "INDEX_2")
-    colnames(index2) <- index2[2,]
-    index2 <- index2[-c(1:2),] |>
-      dplyr::select(Fleet, Yr, Seas, Area, Subseas, Month, Obs, SE, Exp) |>
-      dplyr::mutate(Obs = round(as.numeric(Obs), digits = 3),
-                    SE = round(as.numeric(SE), digits = 3),
-                    Exp = round(as.numeric(Exp), digits = 3))
-    index2_obs <- index2 |>
-      tidyr::pivot_wider(id_cols = Yr,
-                         names_from = Fleet,
-                         values_from = Obs) |>
-      dplyr::rename_at(dplyr::vars(-1), ~ paste(., "obs", sep = "_"))
-    index2_se <- index2 |>
-      tidyr::pivot_wider(id_cols = Yr,names_from = Fleet, values_from = SE) |>
-      dplyr::rename_at(dplyr::vars(-1), ~ paste(., "se", sep = "_"))
-
-    # ind_fleets <- unique(index2$Fleet) |>
-    #   strsplit(split = "_")
-    ind_fleets <- paste("Fleet_", unique(index2$Fleet), sep = "")
-
-    # Recombine fleet names - header name ready
-    # select_one <- function(x){
-    #   a <- lapply(x, tail, -1)
-    #   sapply(a, paste, collapse = " ")
-    # }
-    #
-    # ind_fleets <- select_one(ind_fleets)
-
-    indices <- merge(index2_obs, index2_se) |>
-      dplyr::rename(year = Yr)
-    indices <- indices[, sort(names(indices))] |>
-      dplyr::select("year", dplyr::everything())
-    indices[is.na(indices)] <- "-"
-
-    tab <- indices |>
-      flextable::flextable() |>
-      flextable::set_header_labels(values = c(
-        "Year", rep(c("Obs.", "SE"),
-                    length(ind_fleets))
-      )) |>
-      flextable::add_header_row(top = TRUE, values = c(
-        "", rep(intersect(ind_fleets,fleet_names), each = 2))
-      )
-    tab <- add_theme(tab)
-
-
-  } # close SS3 if statement
-
-  if(model == "BAM"){
-    output <- dget(dat)
-    indices <- output$t.series |>
-      dplyr::select(year, dplyr::contains("U.") & contains(".ob") | contains("cv.U"))
-
-    # Create function to reorder column names so ordered by fleet
-    fleet_names <- function(x){
-      extract_names <- colnames(x)[-1] |>
-        stringr::str_replace("U.", "") |>
-        stringr::str_replace("cv.", "") |>
-        stringr::str_replace(".ob", "") |>
-        unique()
-      return(extract_names)
-    }
-    col_order <- function(x){
-      extract_names <- fleet_names(x)
-      out_df <- x |> dplyr::select(year)
-      for(i in 1:length(extract_names)){
-        cols <- grep(extract_names[i], names(x), value = TRUE)
-        cols_extract <- x[, cols]
-        out_df <- cbind(out_df, cols_extract)
-      }
-      return(out_df)
-    }
-
-    indices <- col_order(indices)
-    indices[is.na(indices)] <- "-"
-    ind_fleets <- fleet_names(indices)
-
-    tab <- indices |>
-      dplyr::mutate(year = as.factor(year)) |>
-      flextable::flextable() |>
-      flextable::set_header_labels(values = c(
-        "Year", rep(c("Obs.", "CV"),
-                    length(ind_fleets))
-      )) |>
-      flextable::add_header_row(top = TRUE, values = c(
-        "", rep(ind_fleets, each = 2))
-      )
-    tab <- add_theme(tab)
-
-  } # close BAM if statement
+      dplyr::filter(grepl("obs", label))
+  } else if (any(unique(output$module_name=="t.series"))) {
+    output <- output |>
+      dplyr::filter(grepl("cpue", label))
+  }
+  fleet_names <- unique(output$fleet)
+  factors <- c("year", "fleet", "fleet_name", "age", "sex", "area", "seas", "season", "time", "era", "subseas", "subseason", "platoon", "platoo","growth_pattern", "gp")
+  # re-structure df for table
+  indices <- output |>
+    dplyr::rename(!!unique(output$label) := estimate,
+                  !!unique(output$uncertainty_label) := uncertainty) |>
+    tidyr::pivot_wider(
+      id_cols = -intersect(colnames(output), factors),
+      names_from = fleet,
+      values_from = c(unique(output$label), unique(output$uncertainty_label))
+    ) # stated internal error for tidyr and asks to report - try again monday
 
   return(tab)
 }
diff --git a/man/plot_landings.Rd b/man/plot_landings.Rd
index deebd377..4b5f8703 100644
--- a/man/plot_landings.Rd
+++ b/man/plot_landings.Rd
@@ -4,12 +4,10 @@
 \alias{plot_landings}
 \title{Plot observed landings by fleet}
 \usage{
-plot_landings(dat, model = "standard", units = NULL)
+plot_landings(dat, units = NULL)
 }
 \arguments{
-\item{dat}{.csv file and path where the data is located of standard stock assessment output}
-
-\item{model}{the model in which the output came from (standard, SS3, BAM, ect)}
+\item{dat}{A data frame returned from `asar::convert_output()`.}
 
 \item{units}{indicate the name of the units of landings as to label the axis}
 }
diff --git a/man/plot_recruitment.Rd b/man/plot_recruitment.Rd
index 6588e88e..c8f6da90 100644
--- a/man/plot_recruitment.Rd
+++ b/man/plot_recruitment.Rd
@@ -6,7 +6,6 @@
 \usage{
 plot_recruitment(
   dat,
-  model = "standard",
   params = FALSE,
   params_only = FALSE,
   units = c(sb = "metric tons", recruitment = "metric tons"),
@@ -20,9 +19,7 @@ plot_recruitment(
 )
 }
 \arguments{
-\item{dat}{.csv file and path where the data is located of standard stock assessment output}
-
-\item{model}{the model in which the output came from (standard, SS3, BAM, ect)}
+\item{dat}{A data frame returned from `asar::convert_output()`.}
 
 \item{params}{Print/export the parameters of the stock recruitment function?}
 
diff --git a/man/plot_spawning_biomass.Rd b/man/plot_spawning_biomass.Rd
index 2599a959..eb014f17 100644
--- a/man/plot_spawning_biomass.Rd
+++ b/man/plot_spawning_biomass.Rd
@@ -2,53 +2,42 @@
 % Please edit documentation in R/plot_spawning_biomass.R
 \name{plot_spawning_biomass}
 \alias{plot_spawning_biomass}
-\title{Plot Spawning Biomass}
+\title{Plot spawning biomass (SB)}
 \usage{
 plot_spawning_biomass(
   dat,
-  model = "standard",
-  show_warnings = FALSE,
-  units = NULL,
-  spawning_biomass_units = NULL,
-  scaled = FALSE,
-  scale_amount = 1000,
-  ref_line = c("target", "MSY", "msy", "unfished"),
+  unit_label = "metric ton",
+  scale_amount = 1,
+  ref_line = c("target", "unfished"),
   end_year = NULL,
   relative = FALSE
 )
 }
 \arguments{
-\item{dat}{.csv file and path where the data is located of standard stock assessment output}
-
-\item{model}{the model in which the output came from (standard, SS3, BAM, ect)}
-
-\item{show_warnings}{Option to suppress warnings}
-
-\item{units}{If units are not available in the output file, in metric tons,
-or are different for SB and R, then report them here starting with SB units
-and following with R units.}
-
-\item{spawning_biomass_units}{units of spawning biomass if different from biomass}
-
-\item{scaled}{T/F; indicate whether the output values for biomass and recruitment are scaled}
+\item{dat}{A data frame returned from `asar::convert_output()`.}
 
 \item{scale_amount}{indicate the exact amount of scale (i.e. 1000)}
 
-\item{ref_line}{choose with reference point to plot a reference line and use
-in relative sb calculations}
-
-\item{end_year}{input the end year of the stock assessment data (not including
-projections). This parameter will be deprecated once the output converter is fully developed.}
+\item{ref_line}{A string specifying the type of reference you want to
+compare spawning biomass to. The default is `"target"`, which looks for
+`"spawning_biomass_target"` in the `"label"` column of `dat`. The actual
+searching in `dat` is case agnostic and will work with either upper- or
+lower-case letters but you must use one of the options specified in the
+default list to ensure that the label on the figure looks correct
+regardless of how it is specified in `dat`.}
 
-\item{relative}{Plot relative spawning biomass. Ref line indicates which reference point to use}
+\item{end_year}{last year of assessment}
 
-\item{biomass_units}{units for biomass}
+\item{relative}{A logical value specifying if the resulting figures should
+be relative spawning biomass. The default is `FALSE`. `ref_line` indicates
+which reference point to use.}
 }
 \value{
-Plot spawning biomass from a stock assessment model as found in a NOAA
-stock assessment report. Units of spawning biomass can either be manually added
-or will be extracted from the provided file if possible. In later releases, model will not
+Plot spawning biomass from the results of an assessment model translated to
+the standard output. The {ggplot2} object is returned for further
+modifications if needed.
 }
 \description{
-Plot Spawning Biomass
+Plot spawning biomass with a reference line as a dashed line. The figure can
+also be made relative to this reference line rather than in absolute units.
 }
diff --git a/man/plot_total_biomass.Rd b/man/plot_total_biomass.Rd
index ce9302a0..9cb7f5a9 100644
--- a/man/plot_total_biomass.Rd
+++ b/man/plot_total_biomass.Rd
@@ -6,7 +6,6 @@
 \usage{
 plot_total_biomass(
   dat,
-  model = "standard",
   show_warnings = FALSE,
   units = NULL,
   scaled = FALSE,
@@ -17,9 +16,7 @@ plot_total_biomass(
 )
 }
 \arguments{
-\item{dat}{.csv file and path where the data is located of standard stock assessment output}
-
-\item{model}{the model in which the output came from (standard, SS3, BAM, ect)}
+\item{dat}{A data frame returned from `asar::convert_output()`.}
 
 \item{show_warnings}{Option to suppress warnings}
 
diff --git a/man/table_indices.Rd b/man/table_indices.Rd
index 37a5f7df..392b69af 100644
--- a/man/table_indices.Rd
+++ b/man/table_indices.Rd
@@ -4,12 +4,10 @@
 \alias{table_indices}
 \title{Create Indices of Abundance Table}
 \usage{
-table_indices(dat, model = "standard")
+table_indices(dat)
 }
 \arguments{
-\item{dat}{.csv file and path where the data is located of standard stock assessment output}
-
-\item{model}{the model in which the output came from (standard, SS3, BAM, ect)}
+\item{dat}{A data frame returned from `asar::convert_output()`.}
 }
 \value{
 Create table of observed annual indices of abundance plus error