```r
library(stockassessment)
load("run/model.RData")

# ---------------------------
# Settings
# ---------------------------
N.SIM <- 10001

# Reference points
Fmsy <- 0.23
Fpa  <- 0.60
Flim <- 0.85
Blim <- 9340
Bpa  <- 11627
Btrigger <- 11627

prev.advice <- 4644
prev.TAC <- 6353

# Storage
FC <- list()

# =========================================================
# FUNCTION TO RUN SCENARIOS (avoids repeating code 4 times)
# =========================================================
run_scenarios <- function(Ry, Fsq, label_prefix){

  yearbase <- 2025
  
  # ---------------------------
  # Initial run to get SSB
  # ---------------------------
  set.seed(12345)
  FC0 <- forecast(
    fit,
    fval=c(Fsq, Fsq, Fmsy, Fmsy),
    rec.years=Ry,
    processNoise=FALSE,
    addTSB=TRUE,
    year.base=yearbase
  )
  
  SSBNext <- attr(FC0,"shorttab")["ssb",3]
  
  # Advice rule
  F_advice <- Fmsy * min(SSBNext / Btrigger, 1)
  
  # ---------------------------
  # Scenarios
  # ---------------------------
  scen <- list(
    
    "SQ" = list(fval=c(Fsq, Fsq, Fsq, Fsq, Fsq, Fsq)),
    
    "F=0" = list(fval=c(Fsq, Fsq, 1e-6, 1e-6, 1e-6, 1e-6)),
    
    "Fmsy" = list(fval=c(Fsq, Fsq, Fmsy, Fmsy, Fmsy, Fmsy)),
    
    "Fpa" = list(fval=c(Fsq, Fsq, Fpa, Fpa, Fpa, Fpa)),
    
    "Flim" = list(fval=c(Fsq, Fsq, Flim, Flim, Flim, Flim)),
    
    "Advice rule" = list(fval=c(Fsq, Fsq, F_advice, Fmsy, Fmsy, Fmsy)),
    
    "SSB=Blim" = list(
      fval=c(Fsq, Fsq, NA, Fmsy),
      nextssb=c(NA, NA, Blim, NA)
    ),
    
    "SSB=Bpa" = list(
      fval=c(Fsq, Fsq, NA, Fmsy),
      nextssb=c(NA, NA, Bpa, NA)
    ),
    
    "Stable SSB" = list(
      fval=c(Fsq, Fsq, NA, Fmsy),
      nextssb=c(NA, NA, SSBNext, NA)
    ),
    
    "TAC" = list(
      fval=c(Fsq, Fsq, NA, Fmsy),
      catchval=c(NA, NA, prev.TAC, NA)
    )
  )
  
  # ---------------------------
  # Run all scenarios
  # ---------------------------
  out <- list()
  
  for(i in seq_along(scen)){
    set.seed(12345)
    
    args <- c(
      scen[[i]],
      list(
        fit = fit,
        rec.years = Ry,
        label = paste(label_prefix, names(scen)[i]),
        processNoise = FALSE,
        addTSB = TRUE,
        year.base = yearbase,
        savesim = TRUE,
        nosim = N.SIM
      )
    )
    
    out[[i]] <- try(do.call(forecast, args))
    
    if(class(out[[i]]) == "try-error"){
      out[[i]] <- NULL
    }
    
    print(paste("Done:", paste(label_prefix, names(scen)[i])))
  }
  
  names(out) <- paste(label_prefix, names(scen))
  
  return(out)
}

# =========================================================
# RUN ALL BLOCKS (your original structure preserved)
# =========================================================

# --- Block 1
FC1 <- run_scenarios(2020:2024, 0.798, "Rec 2020-2024")
FC <- c(FC, FC1)

# --- Block 2
FC2 <- run_scenarios(2020:2024, 0.341, "Rec_F 2020-2024")
FC <- c(FC, FC2)

# --- Block 3
FC3 <- run_scenarios(2015:2024, 0.328, "Rec&F 2015-2024")
FC <- c(FC, FC3)

# --- Block 4
FC4 <- run_scenarios(1997:2024, 0.798, "Rec 1997-2024")
FC <- c(FC, FC4)

# =========================================================
# CREATE SUMMARY TABLE (same as your first script)
# =========================================================

cs <- data.frame(
  Catch = NA,
  SSB = NA,
  F = NA,
  Risk = NA
)

for(i in seq_along(FC)){
  
  cs[i,1] <- attr(FC[[i]],"shorttab")["catch",3]
  cs[i,2] <- attr(FC[[i]],"shorttab")["ssb",4]
  cs[i,3] <- attr(FC[[i]],"shorttab")["fbar",3]
  
  # Risk: P(SSB < Blim)
  cs[i,4] <- 100 * sum(FC[[i]][[4]]$ssb < Blim) / length(FC[[i]][[4]]$ssb)
}

rownames(cs) <- names(FC)

# =========================================================
# SAVE OUTPUT
# =========================================================

save(FC, cs, file="run/forecast.RData")
```