diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
new file mode 100644
index 0000000..a97d04c
--- /dev/null
+++ b/.github/workflows/R-CMD-check.yaml
@@ -0,0 +1,88 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    paths-ignore:
+      - 'README.md'
+  pull_request:
+
+name: build-ss-R-CMD-check
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: macOS-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v3
+
+      - name: Get the latest SS3 executable for linux and move to expected location
+        if: matrix.config.os == 'ubuntu-latest'
+        run: |
+          curl https://api.github.com/repos/nmfs-stock-synthesis/stock-synthesis/releases/latest | grep "browser_download_url" | grep -Eo 'https://[^\"]*' | grep "ss_linux" | xargs wget
+          mv ss_linux ss
+          sudo chmod a+x ss
+          cp ss inst/extdata/ss
+          rm ss
+
+      - name: Get the latest SS3 executable for windows and move to expected location
+        if: matrix.config.os == 'windows-latest' 
+        run:  |
+          Invoke-WebRequest -Uri https://github.com/nmfs-stock-synthesis/stock-synthesis/releases/download/v3.30.21/ss_win.exe -Outfile "D:\a\ss3diags\ss3diags\ss_win.exe"
+          mv ss_win.exe ss.exe
+          cp ss.exe inst/extdata/ss.exe
+
+      - name: Get the latest SS3 executable for macOS and move to expected location
+        if: matrix.config.os == 'macOS-latest' 
+        run:  |
+          curl https://api.github.com/repos/nmfs-stock-synthesis/stock-synthesis/releases/latest | grep "browser_download_url" | grep -Eo 'https://[^\"]*' | grep "ss_osx" | xargs wget
+          mv ss_osx ss
+          sudo chmod a+x ss
+          cp ss inst/extdata/ss
+          rm ss
+      
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: |
+            any::rcmdcheck
+            any::remotes
+
+      - uses: r-lib/actions/check-r-package@v2
+        with: 
+          error-on: '"error"'
+
+      - name: remove ss executable
+        run: |
+          if ["$RUNNER_OS" == "Linux" ]; then
+            rm inst/extdata/ss.exe
+          else 
+            rm inst/extdata/ss
+          fi
+        shell: bash
+
+
+    
+
+
+
diff --git a/.github/workflows/call-r-cmd-check.yml b/.github/workflows/call-r-cmd-check.yml
deleted file mode 100644
index 1f4d5d7..0000000
--- a/.github/workflows/call-r-cmd-check.yml
+++ /dev/null
@@ -1,13 +0,0 @@
-# Run r cmd check
-name: call-r-cmd-check
-# on specifies the build triggers. See more info at https://docs.github.com/en/actions/learn-github-actions/events-that-trigger-workflows
-on:
-# The default build trigger is to run the action on every push and pull request, for any branch
-  push:
-  pull_request:
-  # To run the default repository branch weekly on sunday, uncomment the following 2 lines
-  #schedule:
-    #- cron: '0 0 * * 0'
-jobs:
-  call-workflow:
-    uses: nmfs-fish-tools/ghactions4r/.github/workflows/r-cmd-check.yml@main
diff --git a/DESCRIPTION b/DESCRIPTION
index 4357147..47b3973 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,57 +2,41 @@ Package: ss3diags
 Title: Stock Synthesis Model Diagnostics for Intergated Stock Assessments
 Version: 2.1.1
 Date: 2023-06-07
-Authors@R: 
-  c(person(given = "Henning",
-    family = "Winker",
-    role = c("aut"),
-    email = "henning.winker@ec.europa.eu"),
-    person(given = "Felipe",
-    family = "Carvalho",
-    role = c("aut"),
-    email = "felipe.carvalho@noaa.gov"),
-    person(given = "Massimiliano",
-    family = "Cardinale",
-    role = c("aut"),
-    email = "massimiliano.cardinale@slu.se"),
-    person(given = "Laurence",
-    family = "Kell",
-    role = c("aut"),
-    email = "laurie@seaplusplus.co.uk"),
-    person(given = "Megumi",
-    family = "Oshima",
-    role = c("aut","cre"),
-    email = "megumi.oshima@noaa.gov"),
-    person(given = "Eric",
-    family = "Fletcher",
-    role = c("aut"),
-    email = "eric.fletcher@noaa")
-  ) 
-Description: Compilation of diagnostic functions for Stock Synthesis models, which
-  developed for the manuscript "Using Model Diagnostics in Integrated Stock Assessments".
+Authors@R: c(
+    person("Henning", "Winker", , "henning.winker@ec.europa.eu", role = "aut"),
+    person("Felipe", "Carvalho", , "felipe.carvalho@noaa.gov", role = "aut"),
+    person("Massimiliano", "Cardinale", , "massimiliano.cardinale@slu.se", role = "aut"),
+    person("Laurence", "Kell", , "laurie@seaplusplus.co.uk", role = "aut"),
+    person("Megumi", "Oshima", , "megumi.oshima@noaa.gov", role = c("aut", "cre")),
+    person("Eric", "Fletcher", , "eric.fletcher@noaa", role = "aut")
+  )
+Description: Compilation of diagnostic functions for Stock Synthesis
+    models, which developed for the manuscript "Using Model Diagnostics in
+    Integrated Stock Assessments".
+License: EUPL
+URL: https://github.com/PIFSCstockassessments/ss3diags
+BugReports: https://github.com/PIFSCstockassessments/ss3diags/issues
 Depends:
-  R (>= 3.5)
+    R (>= 3.5)
 Imports:
-    reshape2,
-    mvtnorm,
-    randtests,
-    r4ss(>= 1.46.1),
+    dplyr,
     gplots,
+    lifecycle,
     magrittr,
-    dplyr,
-    rlang,
-    lifecycle
-Remotes: r4ss/r4ss#821
-License: EUPL
+    mvtnorm,
+    r4ss(>= 1.46.1),
+    randtests,
+    reshape2,
+    rlang
+Suggests:
+    knitr,
+    markdown,
+    rmarkdown,
+    testthat (>= 3.0.0)
+Remotes: 
+    r4ss/r4ss#821
 Encoding: UTF-8
 LazyData: True
 LazyDataCompression: gzip
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
-Suggests:
-  testthat (>= 3.0.0),
-  knitr,
-  rmarkdown,
-  markdown
-URL: https://github.com/PIFSCstockassessments/ss3diags
-BugReports: https://github.com/PIFSCstockassessments/ss3diags/issues
diff --git a/README.Rmd b/README.Rmd
index 3b76de0..91b1f7a 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -18,6 +18,7 @@ knitr::opts_chunk$set(
 
 <!-- badges: start -->
 [![R-CMD-check](https://github.com/r4ss/r4ss/workflows/R-CMD-check/badge.svg)](https://github.com/PIFSCstockassessments/ss3diags/actions)
+[![R-CMD-check](https://github.com/PIFSCstockassessments/ss3diags/workflows/R-CMD-check/badge.svg)](https://github.com/PIFSCstockassessments/ss3diags/actions)
 <!-- badges: end -->
 
 ### Build Status
diff --git a/inst/extdata/control.ss b/inst/extdata/control.ss
new file mode 100644
index 0000000..c56e819
--- /dev/null
+++ b/inst/extdata/control.ss
@@ -0,0 +1,314 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#C growth parameters are estimated
+#C spawner-recruitment bias adjustment Not tuned For optimality
+#_data_and_control_files: data.ss // control.ss
+0  # 0 means do not read wtatage.ss; 1 means read and use wtatage.ss and also read and use growth parameters
+1  #_N_Growth_Patterns (Growth Patterns, Morphs, Bio Patterns, GP are terms used interchangeably in SS3)
+1 #_N_platoons_Within_GrowthPattern 
+#_Cond 1 #_Platoon_within/between_stdev_ratio (no read if N_platoons=1)
+#_Cond  1 #vector_platoon_dist_(-1_in_first_val_gives_normal_approx)
+#
+4 # recr_dist_method for parameters:  2=main effects for GP, Area, Settle timing; 3=each Settle entity; 4=none (only when N_GP*Nsettle*pop==1)
+1 # not yet implemented; Future usage: Spawner-Recruitment: 1=global; 2=by area
+1 #  number of recruitment settlement assignments 
+0 # unused option
+#GPattern month  area  age (for each settlement assignment)
+ 1 1 1 0
+#
+#_Cond 0 # N_movement_definitions goes here if Nareas > 1
+#_Cond 1.0 # first age that moves (real age at begin of season, not integer) also cond on do_migration>0
+#_Cond 1 1 1 2 4 10 # example move definition for seas=1, morph=1, source=1 dest=2, age1=4, age2=10
+#
+0 #_Nblock_Patterns
+#_Cond 0 #_blocks_per_pattern 
+# begin and end years of blocks
+#
+# controls for all timevary parameters 
+1 #_time-vary parm bound check (1=warn relative to base parm bounds; 3=no bound check); Also see env (3) and dev (5) options to constrain with base bounds
+#
+# AUTOGEN
+ 0 0 0 0 0 # autogen: 1st element for biology, 2nd for SR, 3rd for Q, 4th reserved, 5th for selex
+# where: 0 = autogen time-varying parms of this category; 1 = read each time-varying parm line; 2 = read then autogen if parm min==-12345
+#
+#_Available timevary codes
+#_Block types: 0: P_block=P_base*exp(TVP); 1: P_block=P_base+TVP; 2: P_block=TVP; 3: P_block=P_block(-1) + TVP
+#_Block_trends: -1: trend bounded by base parm min-max and parms in transformed units (beware); -2: endtrend and infl_year direct values; -3: end and infl as fraction of base range
+#_EnvLinks:  1: P(y)=P_base*exp(TVP*env(y));  2: P(y)=P_base+TVP*env(y);  3: P(y)=f(TVP,env_Zscore) w/ logit to stay in min-max;  4: P(y)=2.0/(1.0+exp(-TVP1*env(y) - TVP2))
+#_DevLinks:  1: P(y)*=exp(dev(y)*dev_se;  2: P(y)+=dev(y)*dev_se;  3: random walk;  4: zero-reverting random walk with rho;  5: like 4 with logit transform to stay in base min-max
+#_DevLinks(more):  21-25 keep last dev for rest of years
+#
+#_Prior_codes:  0=none; 6=normal; 1=symmetric beta; 2=CASAL's beta; 3=lognormal; 4=lognormal with biascorr; 5=gamma
+#
+# setup for M, growth, wt-len, maturity, fecundity, (hermaphro), recr_distr, cohort_grow, (movement), (age error), (catch_mult), sex ratio 
+#_NATMORT
+0 #_natM_type:_0=1Parm; 1=N_breakpoints;_2=Lorenzen;_3=agespecific;_4=agespec_withseasinterpolate;_5=BETA:_Maunder_link_to_maturity;_6=Lorenzen_range
+  #_no additional input for selected M option; read 1P per morph
+#
+1 # GrowthModel: 1=vonBert with L1&L2; 2=Richards with L1&L2; 3=age_specific_K_incr; 4=age_specific_K_decr; 5=age_specific_K_each; 6=NA; 7=NA; 8=growth cessation
+0 #_Age(post-settlement)_for_L1;linear growth below this
+25 #_Growth_Age_for_L2 (999 to use as Linf)
+-999 #_exponential decay for growth above maxage (value should approx initial Z; -999 replicates 3.24; -998 to not allow growth above maxage)
+0  #_placeholder for future growth feature
+#
+0 #_SD_add_to_LAA (set to 0.1 for SS2 V1.x compatibility)
+0 #_CV_Growth_Pattern:  0 CV=f(LAA); 1 CV=F(A); 2 SD=F(LAA); 3 SD=F(A); 4 logSD=F(A)
+#
+1 #_maturity_option:  1=length logistic; 2=age logistic; 3=read age-maturity matrix by growth_pattern; 4=read age-fecundity; 5=disabled; 6=read length-maturity
+1 #_First_Mature_Age
+1 #_fecundity_at_length option:(1)eggs=Wt*(a+b*Wt);(2)eggs=a*L^b;(3)eggs=a*Wt^b; (4)eggs=a+b*L; (5)eggs=a+b*W
+0 #_hermaphroditism option:  0=none; 1=female-to-male age-specific fxn; -1=male-to-female age-specific fxn
+1 #_parameter_offset_approach for M, G, CV_G:  1- direct, no offset**; 2- male=fem_parm*exp(male_parm); 3: male=female*exp(parm) then old=young*exp(parm)
+#_** in option 1, any male parameter with value = 0.0 and phase <0 is set equal to female parameter
+#
+#_growth_parms
+#_ LO HI INIT PRIOR PR_SD PR_type PHASE env_var&link dev_link dev_minyr dev_maxyr dev_PH Block Block_Fxn
+# Sex: 1  BioPattern: 1  NatMort
+ 0.05 0.15 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # NatM_uniform_Fem_GP_1
+# Sex: 1  BioPattern: 1  Growth
+ -10 45 22.769 36 10 0 2 0 0 0 0 0 0 0 # L_at_Amin_Fem_GP_1
+ 40 90 71.8072 70 10 0 4 0 0 0 0 0 0 0 # L_at_Amax_Fem_GP_1
+ 0.05 0.25 0.142165 0.15 0.8 0 4 0 0 0 0 0 0 0 # VonBert_K_Fem_GP_1
+ 0.05 0.25 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # CV_young_Fem_GP_1
+ 0.05 0.25 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # CV_old_Fem_GP_1
+# Sex: 1  BioPattern: 1  WtLen
+ -3 3 2.44e-06 2.44e-06 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_1_Fem_GP_1
+ -3 4 3.34694 3.34694 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_2_Fem_GP_1
+# Sex: 1  BioPattern: 1  Maturity&Fecundity
+ 50 60 55 55 0.8 0 -3 0 0 0 0 0 0 0 # Mat50%_Fem_GP_1
+ -3 3 -0.25 -0.25 0.8 0 -3 0 0 0 0 0 0 0 # Mat_slope_Fem_GP_1
+ -3 3 1 1 0.8 0 -3 0 0 0 0 0 0 0 # Eggs/kg_inter_Fem_GP_1
+ -3 3 0 0 0.8 0 -3 0 0 0 0 0 0 0 # Eggs/kg_slope_wt_Fem_GP_1
+# Sex: 2  BioPattern: 1  NatMort
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # NatM_uniform_Mal_GP_1
+# Sex: 2  BioPattern: 1  Growth
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # L_at_Amin_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # L_at_Amax_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # VonBert_K_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # CV_young_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # CV_old_Mal_GP_1
+# Sex: 2  BioPattern: 1  WtLen
+ -3 3 2.44e-06 2.44e-06 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_1_Mal_GP_1
+ -3 4 3.34694 3.34694 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_2_Mal_GP_1
+# Hermaphroditism
+#  Recruitment Distribution 
+#  Cohort growth dev base
+ 0.1 10 1 1 1 0 -1 0 0 0 0 0 0 0 # CohortGrowDev
+#  Movement
+#  Age Error from parameters
+#  catch multiplier
+#  fraction female, by GP
+ 1e-06 0.999999 0.5 0.5 0.5 0 -99 0 0 0 0 0 0 0 # FracFemale_GP_1
+#  M2 parameter for each predator fleet
+#
+#_no timevary MG parameters
+#
+#_seasonal_effects_on_biology_parms
+ 0 0 0 0 0 0 0 0 0 0 #_femwtlen1,femwtlen2,mat1,mat2,fec1,fec2,Malewtlen1,malewtlen2,L1,K
+#_ LO HI INIT PRIOR PR_SD PR_type PHASE
+#_Cond -2 2 0 0 -1 99 -2 #_placeholder when no seasonal MG parameters
+#
+3 #_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm
+0  # 0/1 to use steepness in initial equ recruitment calculation
+0  #  future feature:  0/1 to make realized sigmaR a function of SR curvature
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn #  parm_name
+             6            12        8.9274          10.3            10             0          1          0          0          0          0          0          0          0 # SR_LN(R0)
+           0.2             1           0.8           0.8             1             0         -4          0          0          0          0          0          0          0 # SR_BH_steep
+             0             2           0.6           0.8           0.8             0         -4          0          0          0          0          0          0          0 # SR_sigmaR
+            -5             5             0             0             1             0         -4          0          0          0          0          0          0          0 # SR_regime
+             0             0             0             0             0             0        -99          0          0          0          0          0          0          0 # SR_autocorr
+#_no timevary SR parameters
+1 #do_recdev:  0=none; 1=devvector (R=F(SSB)+dev); 2=deviations (R=F(SSB)+dev); 3=deviations (R=R0*dev; dev2=R-f(SSB)); 4=like 3 with sum(dev2) adding penalty
+2011 # first year of main recr_devs; early devs can preceed this era
+2022 # last year of main recr_devs; forecast devs start in following year
+2 #_recdev phase 
+1 # (0/1) to read 13 advanced options
+ 0 #_recdev_early_start (0=none; neg value makes relative to recdev_start)
+ -4 #_recdev_early_phase
+ 0 #_forecast_recruitment phase (incl. late recr) (0 value resets to maxphase+1)
+ 1 #_lambda for Fcast_recr_like occurring before endyr+1
+ 1900 #_last_yr_nobias_adj_in_MPD; begin of ramp
+ 1900 #_first_yr_fullbias_adj_in_MPD; begin of plateau
+ 2001 #_last_yr_fullbias_adj_in_MPD
+ 2002 #_end_yr_for_ramp_in_MPD (can be in forecast to shape ramp, but SS3 sets bias_adj to 0.0 for fcast yrs)
+ 1 #_max_bias_adj_in_MPD (typical ~0.8; -3 sets all years to 0.0; -2 sets all non-forecast yrs w/ estimated recdevs to 1.0; -1 sets biasadj=1.0 for all yrs w/ recdevs)
+ 0 #_period of cycles in recruitment (N parms read below)
+ -5 #min rec_dev
+ 5 #max rec_dev
+ 0 #_read_recdevs
+#_end of advanced SR options
+#
+#_placeholder for full parameter lines for recruitment cycles
+# read specified recr devs
+#_Yr Input_value
+#
+# all recruitment deviations
+#  2011R 2012R 2013R 2014R 2015R 2016R 2017R 2018R 2019R 2020R 2021R 2022R 2023F 2024F 2025F 2026F 2027F 2028F 2029F 2030F 2031F 2032F
+#  -0.0600525 -0.0128487 0.134787 0.114101 -0.083254 0.145568 -0.12974 -0.183081 0.0244276 -0.0897981 -0.0781232 0.218013 0 0 0 0 0 0 0 0 0 0
+#
+#Fishing Mortality info 
+0.3 # F ballpark value in units of annual_F
+-2022 # F ballpark year (neg value to disable)
+3 # F_Method:  1=Pope midseason rate; 2=F as parameter; 3=F as hybrid; 4=fleet-specific parm/hybrid (#4 is superset of #2 and #3 and is recommended)
+2.95 # max F (methods 2-4) or harvest fraction (method 1)
+4  # N iterations for tuning in hybrid mode; recommend 3 (faster) to 5 (more precise if many fleets)
+#
+#_initial_F_parms; for each fleet x season that has init_catch; nest season in fleet; count = 1
+#_for unconstrained init_F, use an arbitrary initial catch and set lambda=0 for its logL
+#_ LO HI INIT PRIOR PR_SD  PR_type  PHASE
+ 0.05 0.15 0.1 0.1 0.8 0 -3 # InitF_seas_1_flt_1FISHERY
+#
+# F rates by fleet x season
+# Yr:  2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032
+# seas:  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+# FISHERY 0.220229 0.239068 0.263759 0.287798 0.316758 0.350726 0.404085 0.405179 0.390484 0.37618 0.33635 0.226617 0.0405108 0.0466812 0.051486 0.0553062 0.0583986 0.0609127 0.0629524 0.0646111 0.0659699 0.067099
+#
+#_Q_setup for fleets with cpue or survey data
+#_1:  fleet number
+#_2:  link type: (1=simple q, 1 parm; 2=mirror simple q, 1 mirrored parm; 3=q and power, 2 parm; 4=mirror with offset, 2 parm)
+#_3:  extra input for link, i.e. mirror fleet# or dev index number
+#_4:  0/1 to select extra sd parameter
+#_5:  0/1 for biasadj or not
+#_6:  0/1 to float
+#_   fleet      link link_info  extra_se   biasadj     float  #  fleetname
+         2         1         0         1         0         0  #  SURVEY1
+         3         1         0         0         0         0  #  SURVEY2
+-9999 0 0 0 0 0
+#
+#_Q_parms(if_any);Qunits_are_ln(q)
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn  #  parm_name
+            -7             5      0.476759             0             1             0          1          0          0          0          0          0          0          0  #  LnQ_base_SURVEY1(2)
+             0           0.5             0          0.05             1             0         -4          0          0          0          0          0          0          0  #  Q_extraSD_SURVEY1(2)
+            -7             5      -6.30407             0             1             0          1          0          0          0          0          0          0          0  #  LnQ_base_SURVEY2(3)
+#_no timevary Q parameters
+#
+#_size_selex_patterns
+#Pattern:_0;  parm=0; selex=1.0 for all sizes
+#Pattern:_1;  parm=2; logistic; with 95% width specification
+#Pattern:_5;  parm=2; mirror another size selex; PARMS pick the min-max bin to mirror
+#Pattern:_11; parm=2; selex=1.0  for specified min-max population length bin range
+#Pattern:_15; parm=0; mirror another age or length selex
+#Pattern:_6;  parm=2+special; non-parm len selex
+#Pattern:_43; parm=2+special+2;  like 6, with 2 additional param for scaling (average over bin range)
+#Pattern:_8;  parm=8; double_logistic with smooth transitions and constant above Linf option
+#Pattern:_9;  parm=6; simple 4-parm double logistic with starting length; parm 5 is first length; parm 6=1 does desc as offset
+#Pattern:_21; parm=2+special; non-parm len selex, read as pairs of size, then selex
+#Pattern:_22; parm=4; double_normal as in CASAL
+#Pattern:_23; parm=6; double_normal where final value is directly equal to sp(6) so can be >1.0
+#Pattern:_24; parm=6; double_normal with sel(minL) and sel(maxL), using joiners
+#Pattern:_2;  parm=6; double_normal with sel(minL) and sel(maxL), using joiners, back compatibile version of 24 with 3.30.18 and older
+#Pattern:_25; parm=3; exponential-logistic in length
+#Pattern:_27; parm=special+3; cubic spline in length; parm1==1 resets knots; parm1==2 resets all 
+#Pattern:_42; parm=special+3+2; cubic spline; like 27, with 2 additional param for scaling (average over bin range)
+#_discard_options:_0=none;_1=define_retention;_2=retention&mortality;_3=all_discarded_dead;_4=define_dome-shaped_retention
+#_Pattern Discard Male Special
+ 1 0 0 0 # 1 FISHERY
+ 1 0 0 0 # 2 SURVEY1
+ 0 0 0 0 # 3 SURVEY2
+#
+#_age_selex_patterns
+#Pattern:_0; parm=0; selex=1.0 for ages 0 to maxage
+#Pattern:_10; parm=0; selex=1.0 for ages 1 to maxage
+#Pattern:_11; parm=2; selex=1.0  for specified min-max age
+#Pattern:_12; parm=2; age logistic
+#Pattern:_13; parm=8; age double logistic. Recommend using pattern 18 instead.
+#Pattern:_14; parm=nages+1; age empirical
+#Pattern:_15; parm=0; mirror another age or length selex
+#Pattern:_16; parm=2; Coleraine - Gaussian
+#Pattern:_17; parm=nages+1; empirical as random walk  N parameters to read can be overridden by setting special to non-zero
+#Pattern:_41; parm=2+nages+1; // like 17, with 2 additional param for scaling (average over bin range)
+#Pattern:_18; parm=8; double logistic - smooth transition
+#Pattern:_19; parm=6; simple 4-parm double logistic with starting age
+#Pattern:_20; parm=6; double_normal,using joiners
+#Pattern:_26; parm=3; exponential-logistic in age
+#Pattern:_27; parm=3+special; cubic spline in age; parm1==1 resets knots; parm1==2 resets all 
+#Pattern:_42; parm=2+special+3; // cubic spline; with 2 additional param for scaling (average over bin range)
+#Age patterns entered with value >100 create Min_selage from first digit and pattern from remainder
+#_Pattern Discard Male Special
+ 0 0 0 0 # 1 FISHERY
+ 0 0 0 0 # 2 SURVEY1
+ 11 0 0 0 # 3 SURVEY2
+#
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn  #  parm_name
+# 1   FISHERY LenSelex
+            19            80       56.5138            50            99             0          2          0          0          0          0          0          0          0  #  Size_inflection_FISHERY(1)
+          0.01            60       20.1793            15            99             0          3          0          0          0          0          0          0          0  #  Size_95%width_FISHERY(1)
+# 2   SURVEY1 LenSelex
+            19            70       36.0233            30            99             0          2          0          0          0          0          0          0          0  #  Size_inflection_SURVEY1(2)
+          0.01            60       5.33502            10            99             0          3          0          0          0          0          0          0          0  #  Size_95%width_SURVEY1(2)
+# 3   SURVEY2 LenSelex
+# 1   FISHERY AgeSelex
+# 2   SURVEY1 AgeSelex
+# 3   SURVEY2 AgeSelex
+             0            20             0             5            99             0        -99          0          0          0          0          0          0          0  #  minage@sel=1_SURVEY2(3)
+             0            20             0             6            99             0        -99          0          0          0          0          0          0          0  #  maxage@sel=1_SURVEY2(3)
+#_No_Dirichlet parameters
+#_no timevary selex parameters
+#
+0   #  use 2D_AR1 selectivity(0/1)
+#_no 2D_AR1 selex offset used
+#
+# Tag loss and Tag reporting parameters go next
+0  # TG_custom:  0=no read and autogen if tag data exist; 1=read
+#_Cond -6 6 1 1 2 0.01 -4 0 0 0 0 0 0 0  #_placeholder if no parameters
+#
+# no timevary parameters
+#
+#
+# Input variance adjustments factors: 
+ #_1=add_to_survey_CV
+ #_2=add_to_discard_stddev
+ #_3=add_to_bodywt_CV
+ #_4=mult_by_lencomp_N
+ #_5=mult_by_agecomp_N
+ #_6=mult_by_size-at-age_N
+ #_7=mult_by_generalized_sizecomp
+#_Factor  Fleet  Value
+ -9999   1    0  # terminator
+#
+4 #_maxlambdaphase
+1 #_sd_offset; must be 1 if any growthCV, sigmaR, or survey extraSD is an estimated parameter
+# read 3 changes to default Lambdas (default value is 1.0)
+# Like_comp codes:  1=surv; 2=disc; 3=mnwt; 4=length; 5=age; 6=SizeFreq; 7=sizeage; 8=catch; 9=init_equ_catch; 
+# 10=recrdev; 11=parm_prior; 12=parm_dev; 13=CrashPen; 14=Morphcomp; 15=Tag-comp; 16=Tag-negbin; 17=F_ballpark; 18=initEQregime
+#like_comp fleet  phase  value  sizefreq_method
+ 1 2 2 1 1
+ 4 2 2 1 1
+ 4 2 3 1 1
+-9999  1  1  1  1  #  terminator
+#
+# lambdas (for info only; columns are phases)
+#  0 0 0 0 #_CPUE/survey:_1
+#  1 1 1 1 #_CPUE/survey:_2
+#  1 1 1 1 #_CPUE/survey:_3
+#  1 1 1 1 #_lencomp:_1
+#  1 1 1 1 #_lencomp:_2
+#  0 0 0 0 #_lencomp:_3
+#  1 1 1 1 #_agecomp:_1
+#  1 1 1 1 #_agecomp:_2
+#  0 0 0 0 #_agecomp:_3
+#  1 1 1 1 #_size-age:_1
+#  1 1 1 1 #_size-age:_2
+#  0 0 0 0 #_size-age:_3
+#  1 1 1 1 #_init_equ_catch1
+#  1 1 1 1 #_init_equ_catch2
+#  1 1 1 1 #_init_equ_catch3
+#  1 1 1 1 #_recruitments
+#  1 1 1 1 #_parameter-priors
+#  1 1 1 1 #_parameter-dev-vectors
+#  1 1 1 1 #_crashPenLambda
+#  0 0 0 0 # F_ballpark_lambda
+1 # (0/1/2) read specs for more stddev reporting: 0 = skip, 1 = read specs for reporting stdev for selectivity, size, and numbers, 2 = add options for M,Dyn. Bzero, SmryBio
+ 1 1 -1 5 # Selectivity: (1) 0 to skip or fleet, (2) 1=len/2=age/3=combined, (3) year, (4) N selex bins; NOTE: combined reports in age bins
+ 1 5 # Growth: (1) 0 to skip or growth pattern, (2) growth ages; NOTE: does each sex
+ 1 -1 5 # Numbers-at-age: (1) 0 or area(-1 for all), (2) year, (3) N ages;  NOTE: sums across morphs
+ 5 15 25 35 43 # vector with selex std bins (-1 in first bin to self-generate)
+ 1 2 14 20 20 # vector with growth std ages picks (-1 in first bin to self-generate)
+ 1 2 14 20 20 # vector with NatAge std ages (-1 in first bin to self-generate)
+999
+
diff --git a/inst/extdata/data.ss b/inst/extdata/data.ss
new file mode 100644
index 0000000..07c785a
--- /dev/null
+++ b/inst/extdata/data.ss
@@ -0,0 +1,213 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#_Start_time: Tue Feb 28 13:44:05 2023
+#_echo_input_data
+#C data file for simple example
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+2011 #_StartYr
+2022 #_EndYr
+1 #_Nseas
+ 12 #_months/season
+2 #_Nsubseasons (even number, minimum is 2)
+1 #_spawn_month
+2 #_Ngenders: 1, 2, -1  (use -1 for 1 sex setup with SSB multiplied by female_frac parameter)
+20 #_Nages=accumulator age, first age is always age 0
+1 #_Nareas
+3 #_Nfleets (including surveys)
+#_fleet_type: 1=catch fleet; 2=bycatch only fleet; 3=survey; 4=predator(M2) 
+#_sample_timing: -1 for fishing fleet to use season-long catch-at-age for observations, or 1 to use observation month;  (always 1 for surveys)
+#_fleet_area:  area the fleet/survey operates in 
+#_units of catch:  1=bio; 2=num (ignored for surveys; their units read later)
+#_catch_mult: 0=no; 1=yes
+#_rows are fleets
+#_fleet_type fishery_timing area catch_units need_catch_mult fleetname
+ 1 -1 1 1 0 FISHERY  # 1
+ 3 1 1 2 0 SURVEY1  # 2
+ 3 1 1 2 0 SURVEY2  # 3
+#Bycatch_fleet_input_goes_next
+#a:  fleet index
+#b:  1=include dead bycatch in total dead catch for F0.1 and MSY optimizations and forecast ABC; 2=omit from total catch for these purposes (but still include the mortality)
+#c:  1=Fmult scales with other fleets; 2=bycatch F constant at input value; 3=bycatch F from range of years
+#d:  F or first year of range
+#e:  last year of range
+#f:  not used
+# a   b   c   d   e   f 
+#_Catch data: yr, seas, fleet, catch, catch_se
+#_catch_se:  standard error of log(catch)
+#_NOTE:  catch data is ignored for survey fleets
+-999 1 1 4153.75 0.2
+2011 1 1 10015.9 0.01
+2012 1 1 10026.4 0.01
+2013 1 1 10124.8 0.01
+2014 1 1 10036.8 0.01
+2015 1 1 9969.51 0.01
+2016 1 1 9893.53 0.01
+2017 1 1 10107 0.01
+2018 1 1 9003.8 0.01
+2019 1 1 7882.21 0.01
+2020 1 1 7053.95 0.01
+2021 1 1 6000.76 0.01
+2022 1 1 4024.62 0.01
+-9999 0 0 0 0
+#
+#_CPUE_and_surveyabundance_and_index_observations
+#_Units: 0=numbers; 1=biomass; 2=F; 30=spawnbio; 31=exp(recdev); 36=recdev; 32=spawnbio*recdev; 33=recruitment; 34=depletion(&see Qsetup); 35=parm_dev(&see Qsetup)
+#_Errtype:  -1=normal; 0=lognormal; >0=T
+#_SD_Report: 0=no sdreport; 1=enable sdreport
+#_note that link functions are specified in Q_setup section of control file
+#_Fleet Units Errtype SD_Report
+1 1 0 0 # FISHERY
+2 1 0 1 # SURVEY1
+3 0 0 0 # SURVEY2
+#_yr month fleet obs stderr
+2013 7 2 144745 0.3 #_ SURVEY1
+2016 7 2 63760.3 0.3 #_ SURVEY1
+2019 7 2 59242.9 0.3 #_ SURVEY1
+2022 7 2 49649.7 0.3 #_ SURVEY1
+2011 7 3 11.5668 0.7 #_ SURVEY2
+2012 7 3 13.9955 0.7 #_ SURVEY2
+2013 7 3 12.5783 0.7 #_ SURVEY2
+2014 7 3 16.7479 0.7 #_ SURVEY2
+2015 7 3 7.7595 0.7 #_ SURVEY2
+2016 7 3 9.36206 0.7 #_ SURVEY2
+2017 7 3 16.9079 0.7 #_ SURVEY2
+2018 7 3 6.90196 0.7 #_ SURVEY2
+2019 7 3 14.6227 0.7 #_ SURVEY2
+2020 7 3 7.4737 0.7 #_ SURVEY2
+2021 7 3 7.60085 0.7 #_ SURVEY2
+2022 7 3 15.8286 0.7 #_ SURVEY2
+-9999 1 1 1 1 # terminator for survey observations 
+#
+0 #_N_fleets_with_discard
+#_discard_units (1=same_as_catchunits(bio/num); 2=fraction; 3=numbers)
+#_discard_errtype:  >0 for DF of T-dist(read CV below); 0 for normal with CV; -1 for normal with se; -2 for lognormal; -3 for trunc normal with CV
+# note: only enter units and errtype for fleets with discard 
+# note: discard data is the total for an entire season, so input of month here must be to a month in that season
+#_Fleet units errtype
+# -9999 0 0 0.0 0.0 # terminator for discard data 
+#
+0 #_use meanbodysize_data (0/1)
+#_COND_0 #_DF_for_meanbodysize_T-distribution_like
+# note:  type=1 for mean length; type=2 for mean body weight 
+#_yr month fleet part type obs stderr
+#  -9999 0 0 0 0 0 0 # terminator for mean body size data 
+#
+# set up population length bin structure (note - irrelevant if not using size data and using empirical wtatage
+2 # length bin method: 1=use databins; 2=generate from binwidth,min,max below; 3=read vector
+2 # binwidth for population size comp 
+10 # minimum size in the population (lower edge of first bin and size at age 0.00) 
+94 # maximum size in the population (lower edge of last bin) 
+1 # use length composition data (0/1/2) where 2 invokes new comp_control format
+#_mintailcomp: upper and lower distribution for females and males separately are accumulated until exceeding this level.
+#_addtocomp:  after accumulation of tails; this value added to all bins
+#_combM+F: males and females treated as combined gender below this bin number 
+#_compressbins: accumulate upper tail by this number of bins; acts simultaneous with mintailcomp; set=0 for no forced accumulation
+#_Comp_Error:  0=multinomial, 1=dirichlet using Theta*n, 2=dirichlet using beta, 3=MV_Tweedie
+#_ParmSelect:  consecutive index for dirichlet or MV_Tweedie
+#_minsamplesize: minimum sample size; set to 1 to match 3.24, minimum value is 0.001
+#
+#_Using old format for composition controls
+#_mintailcomp addtocomp combM+F CompressBins CompError ParmSelect minsamplesize
+-1 0.0001 0 0 0 0 0.1 #_fleet:1_FISHERY
+-1 0.0001 0 0 0 0 0.1 #_fleet:2_SURVEY1
+-1 0.0001 0 0 0 0 0.1 #_fleet:3_SURVEY2
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+25 #_N_LengthBins; then enter lower edge of each length bin
+ 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 68 72 76 80 90
+#_yr month fleet sex part Nsamp datavector(female-male)
+ 2011 7 1 3 0 50 0 0 0 0 0 0 0 1 0 1 1 2 0 0 0 0 3 5 1 3 5 1 2 0 1 0 0 1 0 0 0 0 0 0 2 2 1 2 0 3 0 0 0 3 5 2 3 0 0 0
+ 2012 7 1 3 0 50 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 3 1 4 1 5 1 3 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 2 2 3 4 4 4 2 0 1 0 0
+ 2013 7 1 3 0 50 0 0 0 0 0 0 0 0 0 0 2 0 3 0 0 2 0 4 3 1 2 2 0 1 0 0 0 0 0 1 1 1 1 0 2 0 0 1 2 1 2 0 2 6 6 1 1 1 1 0
+ 2014 7 1 3 0 50 0 0 0 0 1 0 0 0 0 1 5 0 1 1 2 3 2 0 2 3 3 1 1 0 1 0 0 1 0 0 0 1 1 0 2 0 0 1 2 1 1 2 1 1 5 1 1 1 1 0
+ 2015 7 1 3 0 50 0 0 0 0 0 1 0 0 0 0 1 0 0 2 0 0 1 0 2 4 4 2 2 0 0 0 0 1 1 0 0 0 0 1 2 0 3 1 1 2 1 1 3 4 2 4 1 2 1 0
+ 2016 7 1 3 0 50 1 0 0 0 0 0 0 0 0 2 0 2 1 2 5 3 1 0 4 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 2 0 0 2 1 2 3 3 1 2 2 0 0
+ 2017 7 1 3 0 50 0 0 0 0 1 0 0 0 0 0 1 0 1 2 3 2 1 3 2 6 2 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 1 2 0 1 1 2 3 3 3 2 2 0 2 0
+ 2018 7 1 3 0 50 0 0 0 0 0 0 1 0 0 1 0 3 2 1 7 3 2 2 0 4 4 1 0 1 0 0 0 0 0 1 0 0 1 0 1 2 0 1 1 2 1 0 0 1 3 2 2 0 0 0
+ 2019 7 1 3 0 50 1 0 0 0 0 0 0 0 1 2 0 0 1 4 5 1 6 1 1 3 4 0 1 0 0 0 0 0 0 0 1 1 0 0 3 0 0 0 1 1 1 0 1 2 3 4 1 0 0 0
+ 2020 7 1 3 0 50 0 0 0 0 1 2 0 0 1 2 1 1 1 2 4 1 2 1 0 6 0 1 0 0 0 0 0 0 1 0 1 0 1 1 1 0 1 0 3 4 0 2 4 1 0 1 1 1 1 0
+ 2021 7 1 3 0 50 0 0 0 1 1 0 0 1 1 2 2 0 1 1 1 2 2 2 1 2 3 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 3 4 2 2 3 1 2 3 3 0 0 0 0
+ 2022 7 1 3 0 50 0 0 0 0 0 0 1 0 0 2 2 0 1 2 2 1 0 2 1 4 0 1 1 0 0 1 0 0 0 0 0 0 1 1 0 2 1 1 3 3 2 3 1 2 2 2 4 1 0 0
+ 2013 7 2 3 0 50 0 0 0 0 0 2 0 3 1 4 2 2 0 4 1 2 2 0 1 3 2 0 0 2 0 0 0 0 1 1 0 0 1 1 1 2 2 0 1 3 0 2 0 0 2 1 1 0 0 0
+ 2016 7 2 3 0 50 0 0 0 0 2 0 1 0 0 2 0 2 2 2 1 0 1 4 1 2 3 1 0 0 0 0 0 0 0 1 1 2 1 4 2 1 0 4 0 2 1 0 2 1 0 3 0 1 0 0
+ 2019 7 2 3 0 50 0 0 0 0 1 1 0 4 1 1 4 0 0 3 1 2 2 1 0 2 0 0 0 0 0 1 0 1 0 0 2 1 2 3 1 3 4 3 1 0 2 1 0 1 0 0 1 0 0 0
+ 2022 7 2 3 0 50 0 0 0 1 3 0 2 0 0 0 0 1 6 2 2 0 4 0 1 4 1 0 1 0 0 0 0 0 0 1 1 0 4 0 1 1 1 3 1 2 2 1 0 1 1 2 0 0 0 0
+-9999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+#
+15 #_N_age_bins
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+2 #_N_ageerror_definitions
+ 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5
+ 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
+ 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5
+ 0.5 0.65 0.67 0.7 0.73 0.76 0.8 0.84 0.88 0.92 0.97 1.03 1.09 1.16 1.23 1.32 1.41 1.51 1.62 1.75 1.89
+#_mintailcomp: upper and lower distribution for females and males separately are accumulated until exceeding this level.
+#_addtocomp:  after accumulation of tails; this value added to all bins
+#_combM+F: males and females treated as combined gender below this bin number 
+#_compressbins: accumulate upper tail by this number of bins; acts simultaneous with mintailcomp; set=0 for no forced accumulation
+#_Comp_Error:  0=multinomial, 1=dirichlet using Theta*n, 2=dirichlet using beta, 3=MV_Tweedie
+#_ParmSelect:  consecutive index for dirichlet or MV_Tweedie
+#_minsamplesize: minimum sample size; set to 1 to match 3.24, minimum value is 0.001
+#
+#_mintailcomp addtocomp combM+F CompressBins CompError ParmSelect minsamplesize
+-1 0.0001 0 0 0 0 0.1 #_fleet:1_FISHERY
+-1 0.0001 0 0 0 0 0.1 #_fleet:2_SURVEY1
+-1 0.0001 0 0 0 0 0.1 #_fleet:3_SURVEY2
+1 #_Lbin_method_for_Age_Data: 1=poplenbins; 2=datalenbins; 3=lengths
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+#_yr month fleet sex part ageerr Lbin_lo Lbin_hi Nsamp datavector(female-male)
+ 2011 7 1 3 0 2 1 -1 25 0 1 0 0 0 2 2 0 0 1 0 0 0 0 4 0 0 2 0 2 1 1 1 1 0 2 0 0 1 4
+ 2012 7 1 3 0 2 1 -1 25 0 0 2 0 1 1 1 2 1 1 0 0 1 0 2 0 1 1 1 0 1 0 2 1 0 2 1 1 1 1
+ 2013 7 1 3 0 2 1 -1 25 1 0 2 1 2 1 1 1 0 1 1 0 0 0 3 0 1 0 0 1 1 3 0 1 0 0 1 0 2 1
+ 2014 7 1 3 0 2 1 -1 25 1 1 1 2 1 1 1 0 0 1 2 0 1 0 3 0 0 1 2 2 0 0 1 0 1 0 1 0 0 2
+ 2015 7 1 3 0 2 1 -1 25 0 0 0 1 0 2 1 0 2 0 2 1 0 1 3 0 0 0 3 1 0 1 0 3 0 0 0 1 0 3
+ 2016 7 1 3 0 2 1 -1 25 1 1 1 1 2 2 3 2 1 1 0 0 2 1 1 0 0 1 0 1 0 0 0 0 0 1 1 0 0 2
+ 2017 7 1 3 0 2 1 -1 25 0 1 2 2 0 2 1 2 0 0 0 0 1 1 1 0 0 1 0 1 1 2 2 0 2 1 0 0 1 1
+ 2018 7 1 3 0 2 1 -1 25 0 1 1 3 3 0 2 0 1 1 1 0 1 0 1 0 1 0 2 0 2 0 0 2 0 1 0 1 0 1
+ 2019 7 1 3 0 2 1 -1 25 0 0 3 2 0 1 2 0 2 0 0 0 0 0 1 0 1 0 1 2 0 1 1 2 3 1 1 0 0 1
+ 2020 7 1 3 0 2 1 -1 25 0 0 2 1 0 0 3 0 1 1 2 0 0 1 2 0 0 1 0 0 0 1 4 0 1 1 1 1 0 2
+ 2021 7 1 3 0 2 1 -1 25 0 1 3 0 2 2 1 0 1 0 0 1 0 0 1 0 0 1 2 2 3 0 1 0 0 1 0 0 1 2
+ 2022 7 1 3 0 2 1 -1 25 0 0 0 1 0 4 1 2 4 2 0 0 1 0 0 0 1 0 0 0 1 3 1 2 0 1 1 0 0 0
+ 2013 7 2 3 0 2 1 -1 25 0 2 0 0 0 3 1 0 0 0 0 0 0 0 1 1 0 2 3 2 2 0 2 0 1 3 0 0 0 2
+ 2016 7 2 3 0 2 1 -1 25 0 2 3 4 0 2 1 3 0 1 1 0 0 0 2 0 0 1 1 2 0 1 0 0 1 0 0 0 0 0
+ 2019 7 2 3 0 2 1 -1 25 2 0 2 4 3 3 1 0 0 1 0 0 0 0 0 1 0 0 2 2 2 0 1 0 1 0 0 0 0 0
+ 2022 7 2 3 0 2 1 -1 25 0 5 1 2 0 0 1 1 0 0 1 0 0 0 1 1 2 1 2 1 1 1 1 0 2 0 0 0 0 1
+-9999  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+#
+1 #_Use_MeanSize-at-Age_obs (0/1)
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+# ageerr codes:  positive means mean length-at-age; negative means mean bodywt_at_age
+#_yr month fleet sex part ageerr ignore datavector(female-male)
+#                                          samplesize(female-male)
+ 2011 7 1 3 0 1 2 31.1001 37.6247 44.5123 49.1577 55.1767 54.8937 57.8764 60.9767 63.407 65.7124 66.0334 66.4773 68.9837 66.6942 71.1547 33.5091 39.7331 44.3663 51.8999 58.3478 56.2476 56.8415 61.0979 62.6753 66.6534 63.2961 68.0089 70.3394 63.6221 64.0669 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2022 7 1 3 0 1 2 34.0009 39.7221 47.1123 45.7161 53.0591 55.0899 59.624 59.0393 60.7682 62.2977 66.8939 69.1026 64.3979 70.8878 70.3055 29.1412 38.5441 43.7912 49.411 53.6436 51.6203 56.6552 58.4004 57.7949 60.8038 60.1757 68.7121 65.4931 71.1546 68.5105 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2011 7 2 3 0 1 2 37.0771 42.2374 45.8496 46.4639 48.0909 53.0877 51.8078 53.4516 59.9566 63.2883 65.2691 61.7347 65.8296 62.7968 65.2321 35.7322 39.2453 41.9728 43.2482 51.1391 51.5264 56.5467 57.8468 59.1265 55.1627 63.584 64.4084 66.1779 67.2428 76.4436 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2022 7 2 3 0 1 2 35.8298 39.9798 43.9049 47.5526 51.4533 55.7381 58.9463 55.3783 64.3036 61.418 59.7164 63.6052 65.2001 65.1717 62.5297 34.9269 40.0321 40.2621 46.1732 49.9184 54.0177 50.9587 56.7169 62.4114 65.9328 63.2853 65.8266 71.7125 69.1363 71.1985 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+-9999  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+#
+0 #_N_environ_variables
+# -2 in yr will subtract mean for that env_var; -1 will subtract mean and divide by stddev (e.g. Z-score)
+#Yr Variable Value
+#
+# Sizefreq data. Defined by method because a fleet can use multiple methods
+0 # N sizefreq methods to read (or -1 for expanded options)
+# 
+0 # do tags (0/1/2); where 2 allows entry of TG_min_recap
+#
+0 #    morphcomp data(0/1) 
+#  Nobs, Nmorphs, mincomp
+#  yr, seas, type, partition, Nsamp, datavector_by_Nmorphs
+#
+0  #  Do dataread for selectivity priors(0/1)
+# Yr, Seas, Fleet,  Age/Size,  Bin,  selex_prior,  prior_sd
+# feature not yet implemented
+#
+999
+
diff --git a/inst/extdata/forecast.ss b/inst/extdata/forecast.ss
new file mode 100644
index 0000000..5db7ba5
--- /dev/null
+++ b/inst/extdata/forecast.ss
@@ -0,0 +1,64 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#C forecast file written by R function SS_writeforecast
+#C rerun model to get more complete formatting in forecast.ss_new
+#C should work with SS version: 3.3
+#C file write time: 2022-06-15 15:53:34
+# for all year entries except rebuilder; enter either: actual year, -999 for styr, 0 for endyr, neg number for rel. endyr
+1 # Benchmarks: 0=skip; 1=calc F_spr,F_btgt,F_msy; 2=calc F_spr,F0.1,F_msy; 3=add F_Blimit; 
+2 # Do_MSY: 1= set to F(SPR); 2=calc F(MSY); 3=set to F(Btgt) or F0.1; 4=set to F(endyr); 5=calc F(MEY) with MSY_unit options
+# if Do_MSY=5, enter MSY_Units; then list fleet_ID, cost/F, price/mt, include_in_Fmey_scaling; # -fleet_ID to fill; -9999 to terminate
+0.5 # SPR target (e.g. 0.40)
+0.4 # Biomass target (e.g. 0.40)
+#_Bmark_years: beg_bio, end_bio, beg_selex, end_selex, beg_relF, end_relF, beg_recr_dist, end_recr_dist, beg_SRparm, end_SRparm (enter actual year, or values of 0 or -integer to be rel. endyr)
+ 0 0 0 0 0 0 0 0 0 0
+#  2022 2022 2022 2022 2022 2022 2022 2022 2022 2022
+# value <0 convert to endyr-value; except -999 converts to start_yr; must be >=start_yr and <=endyr
+1 #Bmark_relF_Basis: 1 = use year range; 2 = set relF same as forecast below
+#
+1 # Forecast: -1=none; 0=simple_1yr; 1=F(SPR); 2=F(MSY) 3=F(Btgt) or F0.1; 4=Ave F (uses first-last relF yrs); 5=input annual F scalar
+# where none and simple require no input after this line; simple sets forecast F same as end year F
+10 # N forecast years 
+0.2 # Fmult (only used for Do_Forecast==5) such that apical_F(f)=Fmult*relF(f)
+#_Fcast_years:  beg_selex, end_selex, beg_relF, end_relF, beg_mean recruits, end_recruits  (enter actual year, or values of 0 or -integer to be rel. endyr)
+ 0 0 0 0 0 0
+#  2022 2022 2022 2022 2022 2022
+0 # Forecast selectivity (0=fcast selex is mean from year range; 1=fcast selectivity from annual time-vary parms)
+1 # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch) 
+# values for top, bottom and buffer exist, but not used when Policy=0
+0.4 # Control rule inflection for constant F (as frac of Bzero, e.g. 0.40); must be > control rule cutoff, or set to -1 to use Bmsy/SSB_unf 
+0.1 # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) 
+0.75 # Buffer:  enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax 
+3 #_N forecast loops (1=OFL only; 2=ABC; 3=get F from forecast ABC catch with allocations applied)
+3 #_First forecast loop with stochastic recruitment
+1 #_Forecast recruitment:  0= spawn_recr; 1=value*spawn_recr_fxn; 2=value*VirginRecr; 3=recent mean from yr range above (need to set phase to -1 in control to get constant recruitment in MCMC)
+1 # value is multiplier of SRR 
+0 #_Forecast loop control #5 (reserved for future bells&whistles) 
+2050  #FirstYear for caps and allocations (should be after years with fixed inputs) 
+0 # stddev of log(realized catch/target catch) in forecast (set value>0.0 to cause active impl_error)
+0 # Do West Coast gfish rebuilder output: 0=no; 1=yes 
+1999 # Rebuilder:  first year catch could have been set to zero (Ydecl)(-1 to set to 1999)
+2002 # Rebuilder:  year for current age structure (Yinit) (-1 to set to endyear+1)
+1 # fleet relative F:  1=use first-last alloc year; 2=read seas, fleet, alloc list below
+# Note that fleet allocation is used directly as average F if Do_Forecast=4 
+2 # basis for fcast catch tuning and for fcast catch caps and allocation  (2=deadbio; 3=retainbio; 5=deadnum; 6=retainnum); NOTE: same units for all fleets
+# Conditional input if relative F choice = 2
+# enter list of:  season,  fleet, relF; if used, terminate with season=-9999
+# 1 1 1
+# -9999 0 0  # terminator for list of relF
+# enter list of: fleet number, max annual catch for fleets with a max; terminate with fleet=-9999
+-9999 -1
+# enter list of area ID and max annual catch; terminate with area=-9999
+-9999 -1
+# enter list of fleet number and allocation group assignment, if any; terminate with fleet=-9999
+-9999 -1
+#_if N allocation groups >0, list year, allocation fraction for each group 
+# list sequentially because read values fill to end of N forecast
+# terminate with -9999 in year field 
+# no allocation groups
+#
+2 # basis for input Fcast catch: -1=read basis with each obs; 2=dead catch; 3=retained catch; 99=input apical_F; NOTE: bio vs num based on fleet's catchunits
+#enter list of Fcast catches or Fa; terminate with line having year=-9999
+#_Yr Seas Fleet Catch(or_F)
+-9999 1 1 0 
+#
+999 # verify end of input 
diff --git a/inst/extdata/starter.ss b/inst/extdata/starter.ss
new file mode 100644
index 0000000..f75b160
--- /dev/null
+++ b/inst/extdata/starter.ss
@@ -0,0 +1,47 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#C starter file written by R function SS_writestarter
+#C rerun model to get more complete formatting in starter.ss_new
+#C should work with SS version: 3.30
+#C file write time: 2022-06-15 15:53:34
+data.ss
+control.ss
+0 # 0=use init values in control file; 1=use ss.par
+0 # run display detail (0,1,2)
+1 # detailed output (0=minimal for data-limited, 1=high (w/ wtatage.ss_new), 2=brief, 3=custom) 
+# custom report options: -100 to start with minimal; -101 to start with all; -number to remove, +number to add, -999 to end
+0 # write 1st iteration details to echoinput.sso file (0,1) 
+0 # write parm values to ParmTrace.sso (0=no,1=good,active; 2=good,all; 3=every_iter,all_parms; 4=every,active)
+0 # write to cumreport.sso (0=no,1=like&timeseries; 2=add survey fits)
+1 # Include prior_like for non-estimated parameters (0,1) 
+1 # Use Soft Boundaries to aid convergence (0,1) (recommended)
+#
+3 # Number of datafiles to produce:  0 turns off all *.ss_new; 1st is data_echo.ss_new, 2nd is data_expval.ss, 3rd and higher are data_boot_**N.ss,
+10 # Turn off estimation for parameters entering after this phase
+#
+0 # MCeval burn interval
+1 # MCeval thin interval
+0 # jitter initial parm value by this fraction
+-1 # min yr for sdreport outputs (-1 for styr); #_2009
+-2 # max yr for sdreport outputs (-1 for endyr+1; -2 for endyr+Nforecastyrs); #_2032
+0 # N individual STD years 
+#vector of year values 
+
+0.0001 # final convergence criteria (e.g. 1.0e-04) 
+0 # retrospective year relative to end year (e.g. -4)
+1 # min age for calc of summary biomass
+1 # Depletion basis:  denom is: 0=skip; 1=X*SPBvirgin; 2=X*SPBmsy; 3=X*SPB_styr; 4=X*SPB_endyr; 5=X*dyn_Bzero;  values>=11 invoke N multiyr (up to 9!) with 10's digit; >100 invokes log(ratio)
+1 # Fraction (X) for Depletion denominator (e.g. 0.4)
+1 # SPR_report_basis:  0=skip; 1=(1-SPR)/(1-SPR_tgt); 2=(1-SPR)/(1-SPR_MSY); 3=(1-SPR)/(1-SPR_Btarget); 4=rawSPR
+3 # F_reporting_units: 0=skip; 1=exploitation(Bio); 2=exploitation(Num); 3=sum(Apical_F's); 4=true F for range of ages; 5=unweighted avg. F for range of ages
+#COND 10 15 #_min and max age over which average F will be calculated with F_reporting=4 or 5
+0 # F_std_basis: 0=raw_annual_F; 1=F/Fspr; 2=F/Fmsy; 3=F/Fbtgt; where F means annual_F; values >=11 invoke N multiyr (up to 9!) with 10's digit; >100 invokes log(ratio)
+0 # MCMC output detail: integer part (0=default; 1=adds obj func components; 2= +write_report_for_each_mceval); and decimal part (added to SR_LN(R0) on first call to mcmc)
+0 # ALK tolerance ***disabled in code (example 0.0001)
+-1 # random number seed for bootstrap data (-1 to use long(time) as seed): # 1677620645
+3.30 # check value for end of file and for version control
diff --git a/model_recipes/test/control.ss b/model_recipes/test/control.ss
new file mode 100644
index 0000000..c56e819
--- /dev/null
+++ b/model_recipes/test/control.ss
@@ -0,0 +1,314 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#C growth parameters are estimated
+#C spawner-recruitment bias adjustment Not tuned For optimality
+#_data_and_control_files: data.ss // control.ss
+0  # 0 means do not read wtatage.ss; 1 means read and use wtatage.ss and also read and use growth parameters
+1  #_N_Growth_Patterns (Growth Patterns, Morphs, Bio Patterns, GP are terms used interchangeably in SS3)
+1 #_N_platoons_Within_GrowthPattern 
+#_Cond 1 #_Platoon_within/between_stdev_ratio (no read if N_platoons=1)
+#_Cond  1 #vector_platoon_dist_(-1_in_first_val_gives_normal_approx)
+#
+4 # recr_dist_method for parameters:  2=main effects for GP, Area, Settle timing; 3=each Settle entity; 4=none (only when N_GP*Nsettle*pop==1)
+1 # not yet implemented; Future usage: Spawner-Recruitment: 1=global; 2=by area
+1 #  number of recruitment settlement assignments 
+0 # unused option
+#GPattern month  area  age (for each settlement assignment)
+ 1 1 1 0
+#
+#_Cond 0 # N_movement_definitions goes here if Nareas > 1
+#_Cond 1.0 # first age that moves (real age at begin of season, not integer) also cond on do_migration>0
+#_Cond 1 1 1 2 4 10 # example move definition for seas=1, morph=1, source=1 dest=2, age1=4, age2=10
+#
+0 #_Nblock_Patterns
+#_Cond 0 #_blocks_per_pattern 
+# begin and end years of blocks
+#
+# controls for all timevary parameters 
+1 #_time-vary parm bound check (1=warn relative to base parm bounds; 3=no bound check); Also see env (3) and dev (5) options to constrain with base bounds
+#
+# AUTOGEN
+ 0 0 0 0 0 # autogen: 1st element for biology, 2nd for SR, 3rd for Q, 4th reserved, 5th for selex
+# where: 0 = autogen time-varying parms of this category; 1 = read each time-varying parm line; 2 = read then autogen if parm min==-12345
+#
+#_Available timevary codes
+#_Block types: 0: P_block=P_base*exp(TVP); 1: P_block=P_base+TVP; 2: P_block=TVP; 3: P_block=P_block(-1) + TVP
+#_Block_trends: -1: trend bounded by base parm min-max and parms in transformed units (beware); -2: endtrend and infl_year direct values; -3: end and infl as fraction of base range
+#_EnvLinks:  1: P(y)=P_base*exp(TVP*env(y));  2: P(y)=P_base+TVP*env(y);  3: P(y)=f(TVP,env_Zscore) w/ logit to stay in min-max;  4: P(y)=2.0/(1.0+exp(-TVP1*env(y) - TVP2))
+#_DevLinks:  1: P(y)*=exp(dev(y)*dev_se;  2: P(y)+=dev(y)*dev_se;  3: random walk;  4: zero-reverting random walk with rho;  5: like 4 with logit transform to stay in base min-max
+#_DevLinks(more):  21-25 keep last dev for rest of years
+#
+#_Prior_codes:  0=none; 6=normal; 1=symmetric beta; 2=CASAL's beta; 3=lognormal; 4=lognormal with biascorr; 5=gamma
+#
+# setup for M, growth, wt-len, maturity, fecundity, (hermaphro), recr_distr, cohort_grow, (movement), (age error), (catch_mult), sex ratio 
+#_NATMORT
+0 #_natM_type:_0=1Parm; 1=N_breakpoints;_2=Lorenzen;_3=agespecific;_4=agespec_withseasinterpolate;_5=BETA:_Maunder_link_to_maturity;_6=Lorenzen_range
+  #_no additional input for selected M option; read 1P per morph
+#
+1 # GrowthModel: 1=vonBert with L1&L2; 2=Richards with L1&L2; 3=age_specific_K_incr; 4=age_specific_K_decr; 5=age_specific_K_each; 6=NA; 7=NA; 8=growth cessation
+0 #_Age(post-settlement)_for_L1;linear growth below this
+25 #_Growth_Age_for_L2 (999 to use as Linf)
+-999 #_exponential decay for growth above maxage (value should approx initial Z; -999 replicates 3.24; -998 to not allow growth above maxage)
+0  #_placeholder for future growth feature
+#
+0 #_SD_add_to_LAA (set to 0.1 for SS2 V1.x compatibility)
+0 #_CV_Growth_Pattern:  0 CV=f(LAA); 1 CV=F(A); 2 SD=F(LAA); 3 SD=F(A); 4 logSD=F(A)
+#
+1 #_maturity_option:  1=length logistic; 2=age logistic; 3=read age-maturity matrix by growth_pattern; 4=read age-fecundity; 5=disabled; 6=read length-maturity
+1 #_First_Mature_Age
+1 #_fecundity_at_length option:(1)eggs=Wt*(a+b*Wt);(2)eggs=a*L^b;(3)eggs=a*Wt^b; (4)eggs=a+b*L; (5)eggs=a+b*W
+0 #_hermaphroditism option:  0=none; 1=female-to-male age-specific fxn; -1=male-to-female age-specific fxn
+1 #_parameter_offset_approach for M, G, CV_G:  1- direct, no offset**; 2- male=fem_parm*exp(male_parm); 3: male=female*exp(parm) then old=young*exp(parm)
+#_** in option 1, any male parameter with value = 0.0 and phase <0 is set equal to female parameter
+#
+#_growth_parms
+#_ LO HI INIT PRIOR PR_SD PR_type PHASE env_var&link dev_link dev_minyr dev_maxyr dev_PH Block Block_Fxn
+# Sex: 1  BioPattern: 1  NatMort
+ 0.05 0.15 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # NatM_uniform_Fem_GP_1
+# Sex: 1  BioPattern: 1  Growth
+ -10 45 22.769 36 10 0 2 0 0 0 0 0 0 0 # L_at_Amin_Fem_GP_1
+ 40 90 71.8072 70 10 0 4 0 0 0 0 0 0 0 # L_at_Amax_Fem_GP_1
+ 0.05 0.25 0.142165 0.15 0.8 0 4 0 0 0 0 0 0 0 # VonBert_K_Fem_GP_1
+ 0.05 0.25 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # CV_young_Fem_GP_1
+ 0.05 0.25 0.1 0.1 0.8 0 -3 0 0 0 0 0 0 0 # CV_old_Fem_GP_1
+# Sex: 1  BioPattern: 1  WtLen
+ -3 3 2.44e-06 2.44e-06 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_1_Fem_GP_1
+ -3 4 3.34694 3.34694 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_2_Fem_GP_1
+# Sex: 1  BioPattern: 1  Maturity&Fecundity
+ 50 60 55 55 0.8 0 -3 0 0 0 0 0 0 0 # Mat50%_Fem_GP_1
+ -3 3 -0.25 -0.25 0.8 0 -3 0 0 0 0 0 0 0 # Mat_slope_Fem_GP_1
+ -3 3 1 1 0.8 0 -3 0 0 0 0 0 0 0 # Eggs/kg_inter_Fem_GP_1
+ -3 3 0 0 0.8 0 -3 0 0 0 0 0 0 0 # Eggs/kg_slope_wt_Fem_GP_1
+# Sex: 2  BioPattern: 1  NatMort
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # NatM_uniform_Mal_GP_1
+# Sex: 2  BioPattern: 1  Growth
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # L_at_Amin_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # L_at_Amax_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # VonBert_K_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # CV_young_Mal_GP_1
+ -3 3 0 0 99 0 -3 0 0 0 0 0 0 0 # CV_old_Mal_GP_1
+# Sex: 2  BioPattern: 1  WtLen
+ -3 3 2.44e-06 2.44e-06 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_1_Mal_GP_1
+ -3 4 3.34694 3.34694 0.8 0 -3 0 0 0 0 0 0 0 # Wtlen_2_Mal_GP_1
+# Hermaphroditism
+#  Recruitment Distribution 
+#  Cohort growth dev base
+ 0.1 10 1 1 1 0 -1 0 0 0 0 0 0 0 # CohortGrowDev
+#  Movement
+#  Age Error from parameters
+#  catch multiplier
+#  fraction female, by GP
+ 1e-06 0.999999 0.5 0.5 0.5 0 -99 0 0 0 0 0 0 0 # FracFemale_GP_1
+#  M2 parameter for each predator fleet
+#
+#_no timevary MG parameters
+#
+#_seasonal_effects_on_biology_parms
+ 0 0 0 0 0 0 0 0 0 0 #_femwtlen1,femwtlen2,mat1,mat2,fec1,fec2,Malewtlen1,malewtlen2,L1,K
+#_ LO HI INIT PRIOR PR_SD PR_type PHASE
+#_Cond -2 2 0 0 -1 99 -2 #_placeholder when no seasonal MG parameters
+#
+3 #_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm
+0  # 0/1 to use steepness in initial equ recruitment calculation
+0  #  future feature:  0/1 to make realized sigmaR a function of SR curvature
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn #  parm_name
+             6            12        8.9274          10.3            10             0          1          0          0          0          0          0          0          0 # SR_LN(R0)
+           0.2             1           0.8           0.8             1             0         -4          0          0          0          0          0          0          0 # SR_BH_steep
+             0             2           0.6           0.8           0.8             0         -4          0          0          0          0          0          0          0 # SR_sigmaR
+            -5             5             0             0             1             0         -4          0          0          0          0          0          0          0 # SR_regime
+             0             0             0             0             0             0        -99          0          0          0          0          0          0          0 # SR_autocorr
+#_no timevary SR parameters
+1 #do_recdev:  0=none; 1=devvector (R=F(SSB)+dev); 2=deviations (R=F(SSB)+dev); 3=deviations (R=R0*dev; dev2=R-f(SSB)); 4=like 3 with sum(dev2) adding penalty
+2011 # first year of main recr_devs; early devs can preceed this era
+2022 # last year of main recr_devs; forecast devs start in following year
+2 #_recdev phase 
+1 # (0/1) to read 13 advanced options
+ 0 #_recdev_early_start (0=none; neg value makes relative to recdev_start)
+ -4 #_recdev_early_phase
+ 0 #_forecast_recruitment phase (incl. late recr) (0 value resets to maxphase+1)
+ 1 #_lambda for Fcast_recr_like occurring before endyr+1
+ 1900 #_last_yr_nobias_adj_in_MPD; begin of ramp
+ 1900 #_first_yr_fullbias_adj_in_MPD; begin of plateau
+ 2001 #_last_yr_fullbias_adj_in_MPD
+ 2002 #_end_yr_for_ramp_in_MPD (can be in forecast to shape ramp, but SS3 sets bias_adj to 0.0 for fcast yrs)
+ 1 #_max_bias_adj_in_MPD (typical ~0.8; -3 sets all years to 0.0; -2 sets all non-forecast yrs w/ estimated recdevs to 1.0; -1 sets biasadj=1.0 for all yrs w/ recdevs)
+ 0 #_period of cycles in recruitment (N parms read below)
+ -5 #min rec_dev
+ 5 #max rec_dev
+ 0 #_read_recdevs
+#_end of advanced SR options
+#
+#_placeholder for full parameter lines for recruitment cycles
+# read specified recr devs
+#_Yr Input_value
+#
+# all recruitment deviations
+#  2011R 2012R 2013R 2014R 2015R 2016R 2017R 2018R 2019R 2020R 2021R 2022R 2023F 2024F 2025F 2026F 2027F 2028F 2029F 2030F 2031F 2032F
+#  -0.0600525 -0.0128487 0.134787 0.114101 -0.083254 0.145568 -0.12974 -0.183081 0.0244276 -0.0897981 -0.0781232 0.218013 0 0 0 0 0 0 0 0 0 0
+#
+#Fishing Mortality info 
+0.3 # F ballpark value in units of annual_F
+-2022 # F ballpark year (neg value to disable)
+3 # F_Method:  1=Pope midseason rate; 2=F as parameter; 3=F as hybrid; 4=fleet-specific parm/hybrid (#4 is superset of #2 and #3 and is recommended)
+2.95 # max F (methods 2-4) or harvest fraction (method 1)
+4  # N iterations for tuning in hybrid mode; recommend 3 (faster) to 5 (more precise if many fleets)
+#
+#_initial_F_parms; for each fleet x season that has init_catch; nest season in fleet; count = 1
+#_for unconstrained init_F, use an arbitrary initial catch and set lambda=0 for its logL
+#_ LO HI INIT PRIOR PR_SD  PR_type  PHASE
+ 0.05 0.15 0.1 0.1 0.8 0 -3 # InitF_seas_1_flt_1FISHERY
+#
+# F rates by fleet x season
+# Yr:  2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032
+# seas:  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+# FISHERY 0.220229 0.239068 0.263759 0.287798 0.316758 0.350726 0.404085 0.405179 0.390484 0.37618 0.33635 0.226617 0.0405108 0.0466812 0.051486 0.0553062 0.0583986 0.0609127 0.0629524 0.0646111 0.0659699 0.067099
+#
+#_Q_setup for fleets with cpue or survey data
+#_1:  fleet number
+#_2:  link type: (1=simple q, 1 parm; 2=mirror simple q, 1 mirrored parm; 3=q and power, 2 parm; 4=mirror with offset, 2 parm)
+#_3:  extra input for link, i.e. mirror fleet# or dev index number
+#_4:  0/1 to select extra sd parameter
+#_5:  0/1 for biasadj or not
+#_6:  0/1 to float
+#_   fleet      link link_info  extra_se   biasadj     float  #  fleetname
+         2         1         0         1         0         0  #  SURVEY1
+         3         1         0         0         0         0  #  SURVEY2
+-9999 0 0 0 0 0
+#
+#_Q_parms(if_any);Qunits_are_ln(q)
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn  #  parm_name
+            -7             5      0.476759             0             1             0          1          0          0          0          0          0          0          0  #  LnQ_base_SURVEY1(2)
+             0           0.5             0          0.05             1             0         -4          0          0          0          0          0          0          0  #  Q_extraSD_SURVEY1(2)
+            -7             5      -6.30407             0             1             0          1          0          0          0          0          0          0          0  #  LnQ_base_SURVEY2(3)
+#_no timevary Q parameters
+#
+#_size_selex_patterns
+#Pattern:_0;  parm=0; selex=1.0 for all sizes
+#Pattern:_1;  parm=2; logistic; with 95% width specification
+#Pattern:_5;  parm=2; mirror another size selex; PARMS pick the min-max bin to mirror
+#Pattern:_11; parm=2; selex=1.0  for specified min-max population length bin range
+#Pattern:_15; parm=0; mirror another age or length selex
+#Pattern:_6;  parm=2+special; non-parm len selex
+#Pattern:_43; parm=2+special+2;  like 6, with 2 additional param for scaling (average over bin range)
+#Pattern:_8;  parm=8; double_logistic with smooth transitions and constant above Linf option
+#Pattern:_9;  parm=6; simple 4-parm double logistic with starting length; parm 5 is first length; parm 6=1 does desc as offset
+#Pattern:_21; parm=2+special; non-parm len selex, read as pairs of size, then selex
+#Pattern:_22; parm=4; double_normal as in CASAL
+#Pattern:_23; parm=6; double_normal where final value is directly equal to sp(6) so can be >1.0
+#Pattern:_24; parm=6; double_normal with sel(minL) and sel(maxL), using joiners
+#Pattern:_2;  parm=6; double_normal with sel(minL) and sel(maxL), using joiners, back compatibile version of 24 with 3.30.18 and older
+#Pattern:_25; parm=3; exponential-logistic in length
+#Pattern:_27; parm=special+3; cubic spline in length; parm1==1 resets knots; parm1==2 resets all 
+#Pattern:_42; parm=special+3+2; cubic spline; like 27, with 2 additional param for scaling (average over bin range)
+#_discard_options:_0=none;_1=define_retention;_2=retention&mortality;_3=all_discarded_dead;_4=define_dome-shaped_retention
+#_Pattern Discard Male Special
+ 1 0 0 0 # 1 FISHERY
+ 1 0 0 0 # 2 SURVEY1
+ 0 0 0 0 # 3 SURVEY2
+#
+#_age_selex_patterns
+#Pattern:_0; parm=0; selex=1.0 for ages 0 to maxage
+#Pattern:_10; parm=0; selex=1.0 for ages 1 to maxage
+#Pattern:_11; parm=2; selex=1.0  for specified min-max age
+#Pattern:_12; parm=2; age logistic
+#Pattern:_13; parm=8; age double logistic. Recommend using pattern 18 instead.
+#Pattern:_14; parm=nages+1; age empirical
+#Pattern:_15; parm=0; mirror another age or length selex
+#Pattern:_16; parm=2; Coleraine - Gaussian
+#Pattern:_17; parm=nages+1; empirical as random walk  N parameters to read can be overridden by setting special to non-zero
+#Pattern:_41; parm=2+nages+1; // like 17, with 2 additional param for scaling (average over bin range)
+#Pattern:_18; parm=8; double logistic - smooth transition
+#Pattern:_19; parm=6; simple 4-parm double logistic with starting age
+#Pattern:_20; parm=6; double_normal,using joiners
+#Pattern:_26; parm=3; exponential-logistic in age
+#Pattern:_27; parm=3+special; cubic spline in age; parm1==1 resets knots; parm1==2 resets all 
+#Pattern:_42; parm=2+special+3; // cubic spline; with 2 additional param for scaling (average over bin range)
+#Age patterns entered with value >100 create Min_selage from first digit and pattern from remainder
+#_Pattern Discard Male Special
+ 0 0 0 0 # 1 FISHERY
+ 0 0 0 0 # 2 SURVEY1
+ 11 0 0 0 # 3 SURVEY2
+#
+#_          LO            HI          INIT         PRIOR         PR_SD       PR_type      PHASE    env-var    use_dev   dev_mnyr   dev_mxyr     dev_PH      Block    Blk_Fxn  #  parm_name
+# 1   FISHERY LenSelex
+            19            80       56.5138            50            99             0          2          0          0          0          0          0          0          0  #  Size_inflection_FISHERY(1)
+          0.01            60       20.1793            15            99             0          3          0          0          0          0          0          0          0  #  Size_95%width_FISHERY(1)
+# 2   SURVEY1 LenSelex
+            19            70       36.0233            30            99             0          2          0          0          0          0          0          0          0  #  Size_inflection_SURVEY1(2)
+          0.01            60       5.33502            10            99             0          3          0          0          0          0          0          0          0  #  Size_95%width_SURVEY1(2)
+# 3   SURVEY2 LenSelex
+# 1   FISHERY AgeSelex
+# 2   SURVEY1 AgeSelex
+# 3   SURVEY2 AgeSelex
+             0            20             0             5            99             0        -99          0          0          0          0          0          0          0  #  minage@sel=1_SURVEY2(3)
+             0            20             0             6            99             0        -99          0          0          0          0          0          0          0  #  maxage@sel=1_SURVEY2(3)
+#_No_Dirichlet parameters
+#_no timevary selex parameters
+#
+0   #  use 2D_AR1 selectivity(0/1)
+#_no 2D_AR1 selex offset used
+#
+# Tag loss and Tag reporting parameters go next
+0  # TG_custom:  0=no read and autogen if tag data exist; 1=read
+#_Cond -6 6 1 1 2 0.01 -4 0 0 0 0 0 0 0  #_placeholder if no parameters
+#
+# no timevary parameters
+#
+#
+# Input variance adjustments factors: 
+ #_1=add_to_survey_CV
+ #_2=add_to_discard_stddev
+ #_3=add_to_bodywt_CV
+ #_4=mult_by_lencomp_N
+ #_5=mult_by_agecomp_N
+ #_6=mult_by_size-at-age_N
+ #_7=mult_by_generalized_sizecomp
+#_Factor  Fleet  Value
+ -9999   1    0  # terminator
+#
+4 #_maxlambdaphase
+1 #_sd_offset; must be 1 if any growthCV, sigmaR, or survey extraSD is an estimated parameter
+# read 3 changes to default Lambdas (default value is 1.0)
+# Like_comp codes:  1=surv; 2=disc; 3=mnwt; 4=length; 5=age; 6=SizeFreq; 7=sizeage; 8=catch; 9=init_equ_catch; 
+# 10=recrdev; 11=parm_prior; 12=parm_dev; 13=CrashPen; 14=Morphcomp; 15=Tag-comp; 16=Tag-negbin; 17=F_ballpark; 18=initEQregime
+#like_comp fleet  phase  value  sizefreq_method
+ 1 2 2 1 1
+ 4 2 2 1 1
+ 4 2 3 1 1
+-9999  1  1  1  1  #  terminator
+#
+# lambdas (for info only; columns are phases)
+#  0 0 0 0 #_CPUE/survey:_1
+#  1 1 1 1 #_CPUE/survey:_2
+#  1 1 1 1 #_CPUE/survey:_3
+#  1 1 1 1 #_lencomp:_1
+#  1 1 1 1 #_lencomp:_2
+#  0 0 0 0 #_lencomp:_3
+#  1 1 1 1 #_agecomp:_1
+#  1 1 1 1 #_agecomp:_2
+#  0 0 0 0 #_agecomp:_3
+#  1 1 1 1 #_size-age:_1
+#  1 1 1 1 #_size-age:_2
+#  0 0 0 0 #_size-age:_3
+#  1 1 1 1 #_init_equ_catch1
+#  1 1 1 1 #_init_equ_catch2
+#  1 1 1 1 #_init_equ_catch3
+#  1 1 1 1 #_recruitments
+#  1 1 1 1 #_parameter-priors
+#  1 1 1 1 #_parameter-dev-vectors
+#  1 1 1 1 #_crashPenLambda
+#  0 0 0 0 # F_ballpark_lambda
+1 # (0/1/2) read specs for more stddev reporting: 0 = skip, 1 = read specs for reporting stdev for selectivity, size, and numbers, 2 = add options for M,Dyn. Bzero, SmryBio
+ 1 1 -1 5 # Selectivity: (1) 0 to skip or fleet, (2) 1=len/2=age/3=combined, (3) year, (4) N selex bins; NOTE: combined reports in age bins
+ 1 5 # Growth: (1) 0 to skip or growth pattern, (2) growth ages; NOTE: does each sex
+ 1 -1 5 # Numbers-at-age: (1) 0 or area(-1 for all), (2) year, (3) N ages;  NOTE: sums across morphs
+ 5 15 25 35 43 # vector with selex std bins (-1 in first bin to self-generate)
+ 1 2 14 20 20 # vector with growth std ages picks (-1 in first bin to self-generate)
+ 1 2 14 20 20 # vector with NatAge std ages (-1 in first bin to self-generate)
+999
+
diff --git a/model_recipes/test/data.ss b/model_recipes/test/data.ss
new file mode 100644
index 0000000..07c785a
--- /dev/null
+++ b/model_recipes/test/data.ss
@@ -0,0 +1,213 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#_Start_time: Tue Feb 28 13:44:05 2023
+#_echo_input_data
+#C data file for simple example
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+2011 #_StartYr
+2022 #_EndYr
+1 #_Nseas
+ 12 #_months/season
+2 #_Nsubseasons (even number, minimum is 2)
+1 #_spawn_month
+2 #_Ngenders: 1, 2, -1  (use -1 for 1 sex setup with SSB multiplied by female_frac parameter)
+20 #_Nages=accumulator age, first age is always age 0
+1 #_Nareas
+3 #_Nfleets (including surveys)
+#_fleet_type: 1=catch fleet; 2=bycatch only fleet; 3=survey; 4=predator(M2) 
+#_sample_timing: -1 for fishing fleet to use season-long catch-at-age for observations, or 1 to use observation month;  (always 1 for surveys)
+#_fleet_area:  area the fleet/survey operates in 
+#_units of catch:  1=bio; 2=num (ignored for surveys; their units read later)
+#_catch_mult: 0=no; 1=yes
+#_rows are fleets
+#_fleet_type fishery_timing area catch_units need_catch_mult fleetname
+ 1 -1 1 1 0 FISHERY  # 1
+ 3 1 1 2 0 SURVEY1  # 2
+ 3 1 1 2 0 SURVEY2  # 3
+#Bycatch_fleet_input_goes_next
+#a:  fleet index
+#b:  1=include dead bycatch in total dead catch for F0.1 and MSY optimizations and forecast ABC; 2=omit from total catch for these purposes (but still include the mortality)
+#c:  1=Fmult scales with other fleets; 2=bycatch F constant at input value; 3=bycatch F from range of years
+#d:  F or first year of range
+#e:  last year of range
+#f:  not used
+# a   b   c   d   e   f 
+#_Catch data: yr, seas, fleet, catch, catch_se
+#_catch_se:  standard error of log(catch)
+#_NOTE:  catch data is ignored for survey fleets
+-999 1 1 4153.75 0.2
+2011 1 1 10015.9 0.01
+2012 1 1 10026.4 0.01
+2013 1 1 10124.8 0.01
+2014 1 1 10036.8 0.01
+2015 1 1 9969.51 0.01
+2016 1 1 9893.53 0.01
+2017 1 1 10107 0.01
+2018 1 1 9003.8 0.01
+2019 1 1 7882.21 0.01
+2020 1 1 7053.95 0.01
+2021 1 1 6000.76 0.01
+2022 1 1 4024.62 0.01
+-9999 0 0 0 0
+#
+#_CPUE_and_surveyabundance_and_index_observations
+#_Units: 0=numbers; 1=biomass; 2=F; 30=spawnbio; 31=exp(recdev); 36=recdev; 32=spawnbio*recdev; 33=recruitment; 34=depletion(&see Qsetup); 35=parm_dev(&see Qsetup)
+#_Errtype:  -1=normal; 0=lognormal; >0=T
+#_SD_Report: 0=no sdreport; 1=enable sdreport
+#_note that link functions are specified in Q_setup section of control file
+#_Fleet Units Errtype SD_Report
+1 1 0 0 # FISHERY
+2 1 0 1 # SURVEY1
+3 0 0 0 # SURVEY2
+#_yr month fleet obs stderr
+2013 7 2 144745 0.3 #_ SURVEY1
+2016 7 2 63760.3 0.3 #_ SURVEY1
+2019 7 2 59242.9 0.3 #_ SURVEY1
+2022 7 2 49649.7 0.3 #_ SURVEY1
+2011 7 3 11.5668 0.7 #_ SURVEY2
+2012 7 3 13.9955 0.7 #_ SURVEY2
+2013 7 3 12.5783 0.7 #_ SURVEY2
+2014 7 3 16.7479 0.7 #_ SURVEY2
+2015 7 3 7.7595 0.7 #_ SURVEY2
+2016 7 3 9.36206 0.7 #_ SURVEY2
+2017 7 3 16.9079 0.7 #_ SURVEY2
+2018 7 3 6.90196 0.7 #_ SURVEY2
+2019 7 3 14.6227 0.7 #_ SURVEY2
+2020 7 3 7.4737 0.7 #_ SURVEY2
+2021 7 3 7.60085 0.7 #_ SURVEY2
+2022 7 3 15.8286 0.7 #_ SURVEY2
+-9999 1 1 1 1 # terminator for survey observations 
+#
+0 #_N_fleets_with_discard
+#_discard_units (1=same_as_catchunits(bio/num); 2=fraction; 3=numbers)
+#_discard_errtype:  >0 for DF of T-dist(read CV below); 0 for normal with CV; -1 for normal with se; -2 for lognormal; -3 for trunc normal with CV
+# note: only enter units and errtype for fleets with discard 
+# note: discard data is the total for an entire season, so input of month here must be to a month in that season
+#_Fleet units errtype
+# -9999 0 0 0.0 0.0 # terminator for discard data 
+#
+0 #_use meanbodysize_data (0/1)
+#_COND_0 #_DF_for_meanbodysize_T-distribution_like
+# note:  type=1 for mean length; type=2 for mean body weight 
+#_yr month fleet part type obs stderr
+#  -9999 0 0 0 0 0 0 # terminator for mean body size data 
+#
+# set up population length bin structure (note - irrelevant if not using size data and using empirical wtatage
+2 # length bin method: 1=use databins; 2=generate from binwidth,min,max below; 3=read vector
+2 # binwidth for population size comp 
+10 # minimum size in the population (lower edge of first bin and size at age 0.00) 
+94 # maximum size in the population (lower edge of last bin) 
+1 # use length composition data (0/1/2) where 2 invokes new comp_control format
+#_mintailcomp: upper and lower distribution for females and males separately are accumulated until exceeding this level.
+#_addtocomp:  after accumulation of tails; this value added to all bins
+#_combM+F: males and females treated as combined gender below this bin number 
+#_compressbins: accumulate upper tail by this number of bins; acts simultaneous with mintailcomp; set=0 for no forced accumulation
+#_Comp_Error:  0=multinomial, 1=dirichlet using Theta*n, 2=dirichlet using beta, 3=MV_Tweedie
+#_ParmSelect:  consecutive index for dirichlet or MV_Tweedie
+#_minsamplesize: minimum sample size; set to 1 to match 3.24, minimum value is 0.001
+#
+#_Using old format for composition controls
+#_mintailcomp addtocomp combM+F CompressBins CompError ParmSelect minsamplesize
+-1 0.0001 0 0 0 0 0.1 #_fleet:1_FISHERY
+-1 0.0001 0 0 0 0 0.1 #_fleet:2_SURVEY1
+-1 0.0001 0 0 0 0 0.1 #_fleet:3_SURVEY2
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+25 #_N_LengthBins; then enter lower edge of each length bin
+ 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 68 72 76 80 90
+#_yr month fleet sex part Nsamp datavector(female-male)
+ 2011 7 1 3 0 50 0 0 0 0 0 0 0 1 0 1 1 2 0 0 0 0 3 5 1 3 5 1 2 0 1 0 0 1 0 0 0 0 0 0 2 2 1 2 0 3 0 0 0 3 5 2 3 0 0 0
+ 2012 7 1 3 0 50 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 3 1 4 1 5 1 3 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 2 2 3 4 4 4 2 0 1 0 0
+ 2013 7 1 3 0 50 0 0 0 0 0 0 0 0 0 0 2 0 3 0 0 2 0 4 3 1 2 2 0 1 0 0 0 0 0 1 1 1 1 0 2 0 0 1 2 1 2 0 2 6 6 1 1 1 1 0
+ 2014 7 1 3 0 50 0 0 0 0 1 0 0 0 0 1 5 0 1 1 2 3 2 0 2 3 3 1 1 0 1 0 0 1 0 0 0 1 1 0 2 0 0 1 2 1 1 2 1 1 5 1 1 1 1 0
+ 2015 7 1 3 0 50 0 0 0 0 0 1 0 0 0 0 1 0 0 2 0 0 1 0 2 4 4 2 2 0 0 0 0 1 1 0 0 0 0 1 2 0 3 1 1 2 1 1 3 4 2 4 1 2 1 0
+ 2016 7 1 3 0 50 1 0 0 0 0 0 0 0 0 2 0 2 1 2 5 3 1 0 4 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 2 0 0 2 1 2 3 3 1 2 2 0 0
+ 2017 7 1 3 0 50 0 0 0 0 1 0 0 0 0 0 1 0 1 2 3 2 1 3 2 6 2 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 1 2 0 1 1 2 3 3 3 2 2 0 2 0
+ 2018 7 1 3 0 50 0 0 0 0 0 0 1 0 0 1 0 3 2 1 7 3 2 2 0 4 4 1 0 1 0 0 0 0 0 1 0 0 1 0 1 2 0 1 1 2 1 0 0 1 3 2 2 0 0 0
+ 2019 7 1 3 0 50 1 0 0 0 0 0 0 0 1 2 0 0 1 4 5 1 6 1 1 3 4 0 1 0 0 0 0 0 0 0 1 1 0 0 3 0 0 0 1 1 1 0 1 2 3 4 1 0 0 0
+ 2020 7 1 3 0 50 0 0 0 0 1 2 0 0 1 2 1 1 1 2 4 1 2 1 0 6 0 1 0 0 0 0 0 0 1 0 1 0 1 1 1 0 1 0 3 4 0 2 4 1 0 1 1 1 1 0
+ 2021 7 1 3 0 50 0 0 0 1 1 0 0 1 1 2 2 0 1 1 1 2 2 2 1 2 3 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 3 4 2 2 3 1 2 3 3 0 0 0 0
+ 2022 7 1 3 0 50 0 0 0 0 0 0 1 0 0 2 2 0 1 2 2 1 0 2 1 4 0 1 1 0 0 1 0 0 0 0 0 0 1 1 0 2 1 1 3 3 2 3 1 2 2 2 4 1 0 0
+ 2013 7 2 3 0 50 0 0 0 0 0 2 0 3 1 4 2 2 0 4 1 2 2 0 1 3 2 0 0 2 0 0 0 0 1 1 0 0 1 1 1 2 2 0 1 3 0 2 0 0 2 1 1 0 0 0
+ 2016 7 2 3 0 50 0 0 0 0 2 0 1 0 0 2 0 2 2 2 1 0 1 4 1 2 3 1 0 0 0 0 0 0 0 1 1 2 1 4 2 1 0 4 0 2 1 0 2 1 0 3 0 1 0 0
+ 2019 7 2 3 0 50 0 0 0 0 1 1 0 4 1 1 4 0 0 3 1 2 2 1 0 2 0 0 0 0 0 1 0 1 0 0 2 1 2 3 1 3 4 3 1 0 2 1 0 1 0 0 1 0 0 0
+ 2022 7 2 3 0 50 0 0 0 1 3 0 2 0 0 0 0 1 6 2 2 0 4 0 1 4 1 0 1 0 0 0 0 0 0 1 1 0 4 0 1 1 1 3 1 2 2 1 0 1 1 2 0 0 0 0
+-9999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+#
+15 #_N_age_bins
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+2 #_N_ageerror_definitions
+ 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5
+ 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
+ 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5
+ 0.5 0.65 0.67 0.7 0.73 0.76 0.8 0.84 0.88 0.92 0.97 1.03 1.09 1.16 1.23 1.32 1.41 1.51 1.62 1.75 1.89
+#_mintailcomp: upper and lower distribution for females and males separately are accumulated until exceeding this level.
+#_addtocomp:  after accumulation of tails; this value added to all bins
+#_combM+F: males and females treated as combined gender below this bin number 
+#_compressbins: accumulate upper tail by this number of bins; acts simultaneous with mintailcomp; set=0 for no forced accumulation
+#_Comp_Error:  0=multinomial, 1=dirichlet using Theta*n, 2=dirichlet using beta, 3=MV_Tweedie
+#_ParmSelect:  consecutive index for dirichlet or MV_Tweedie
+#_minsamplesize: minimum sample size; set to 1 to match 3.24, minimum value is 0.001
+#
+#_mintailcomp addtocomp combM+F CompressBins CompError ParmSelect minsamplesize
+-1 0.0001 0 0 0 0 0.1 #_fleet:1_FISHERY
+-1 0.0001 0 0 0 0 0.1 #_fleet:2_SURVEY1
+-1 0.0001 0 0 0 0 0.1 #_fleet:3_SURVEY2
+1 #_Lbin_method_for_Age_Data: 1=poplenbins; 2=datalenbins; 3=lengths
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+#_yr month fleet sex part ageerr Lbin_lo Lbin_hi Nsamp datavector(female-male)
+ 2011 7 1 3 0 2 1 -1 25 0 1 0 0 0 2 2 0 0 1 0 0 0 0 4 0 0 2 0 2 1 1 1 1 0 2 0 0 1 4
+ 2012 7 1 3 0 2 1 -1 25 0 0 2 0 1 1 1 2 1 1 0 0 1 0 2 0 1 1 1 0 1 0 2 1 0 2 1 1 1 1
+ 2013 7 1 3 0 2 1 -1 25 1 0 2 1 2 1 1 1 0 1 1 0 0 0 3 0 1 0 0 1 1 3 0 1 0 0 1 0 2 1
+ 2014 7 1 3 0 2 1 -1 25 1 1 1 2 1 1 1 0 0 1 2 0 1 0 3 0 0 1 2 2 0 0 1 0 1 0 1 0 0 2
+ 2015 7 1 3 0 2 1 -1 25 0 0 0 1 0 2 1 0 2 0 2 1 0 1 3 0 0 0 3 1 0 1 0 3 0 0 0 1 0 3
+ 2016 7 1 3 0 2 1 -1 25 1 1 1 1 2 2 3 2 1 1 0 0 2 1 1 0 0 1 0 1 0 0 0 0 0 1 1 0 0 2
+ 2017 7 1 3 0 2 1 -1 25 0 1 2 2 0 2 1 2 0 0 0 0 1 1 1 0 0 1 0 1 1 2 2 0 2 1 0 0 1 1
+ 2018 7 1 3 0 2 1 -1 25 0 1 1 3 3 0 2 0 1 1 1 0 1 0 1 0 1 0 2 0 2 0 0 2 0 1 0 1 0 1
+ 2019 7 1 3 0 2 1 -1 25 0 0 3 2 0 1 2 0 2 0 0 0 0 0 1 0 1 0 1 2 0 1 1 2 3 1 1 0 0 1
+ 2020 7 1 3 0 2 1 -1 25 0 0 2 1 0 0 3 0 1 1 2 0 0 1 2 0 0 1 0 0 0 1 4 0 1 1 1 1 0 2
+ 2021 7 1 3 0 2 1 -1 25 0 1 3 0 2 2 1 0 1 0 0 1 0 0 1 0 0 1 2 2 3 0 1 0 0 1 0 0 1 2
+ 2022 7 1 3 0 2 1 -1 25 0 0 0 1 0 4 1 2 4 2 0 0 1 0 0 0 1 0 0 0 1 3 1 2 0 1 1 0 0 0
+ 2013 7 2 3 0 2 1 -1 25 0 2 0 0 0 3 1 0 0 0 0 0 0 0 1 1 0 2 3 2 2 0 2 0 1 3 0 0 0 2
+ 2016 7 2 3 0 2 1 -1 25 0 2 3 4 0 2 1 3 0 1 1 0 0 0 2 0 0 1 1 2 0 1 0 0 1 0 0 0 0 0
+ 2019 7 2 3 0 2 1 -1 25 2 0 2 4 3 3 1 0 0 1 0 0 0 0 0 1 0 0 2 2 2 0 1 0 1 0 0 0 0 0
+ 2022 7 2 3 0 2 1 -1 25 0 5 1 2 0 0 1 1 0 0 1 0 0 0 1 1 2 1 2 1 1 1 1 0 2 0 0 0 0 1
+-9999  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+#
+1 #_Use_MeanSize-at-Age_obs (0/1)
+# sex codes:  0=combined; 1=use female only; 2=use male only; 3=use both as joint sexxlength distribution
+# partition codes:  (0=combined; 1=discard; 2=retained
+# ageerr codes:  positive means mean length-at-age; negative means mean bodywt_at_age
+#_yr month fleet sex part ageerr ignore datavector(female-male)
+#                                          samplesize(female-male)
+ 2011 7 1 3 0 1 2 31.1001 37.6247 44.5123 49.1577 55.1767 54.8937 57.8764 60.9767 63.407 65.7124 66.0334 66.4773 68.9837 66.6942 71.1547 33.5091 39.7331 44.3663 51.8999 58.3478 56.2476 56.8415 61.0979 62.6753 66.6534 63.2961 68.0089 70.3394 63.6221 64.0669 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2022 7 1 3 0 1 2 34.0009 39.7221 47.1123 45.7161 53.0591 55.0899 59.624 59.0393 60.7682 62.2977 66.8939 69.1026 64.3979 70.8878 70.3055 29.1412 38.5441 43.7912 49.411 53.6436 51.6203 56.6552 58.4004 57.7949 60.8038 60.1757 68.7121 65.4931 71.1546 68.5105 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2011 7 2 3 0 1 2 37.0771 42.2374 45.8496 46.4639 48.0909 53.0877 51.8078 53.4516 59.9566 63.2883 65.2691 61.7347 65.8296 62.7968 65.2321 35.7322 39.2453 41.9728 43.2482 51.1391 51.5264 56.5467 57.8468 59.1265 55.1627 63.584 64.4084 66.1779 67.2428 76.4436 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+ 2022 7 2 3 0 1 2 35.8298 39.9798 43.9049 47.5526 51.4533 55.7381 58.9463 55.3783 64.3036 61.418 59.7164 63.6052 65.2001 65.1717 62.5297 34.9269 40.0321 40.2621 46.1732 49.9184 54.0177 50.9587 56.7169 62.4114 65.9328 63.2853 65.8266 71.7125 69.1363 71.1985 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+-9999  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+#
+0 #_N_environ_variables
+# -2 in yr will subtract mean for that env_var; -1 will subtract mean and divide by stddev (e.g. Z-score)
+#Yr Variable Value
+#
+# Sizefreq data. Defined by method because a fleet can use multiple methods
+0 # N sizefreq methods to read (or -1 for expanded options)
+# 
+0 # do tags (0/1/2); where 2 allows entry of TG_min_recap
+#
+0 #    morphcomp data(0/1) 
+#  Nobs, Nmorphs, mincomp
+#  yr, seas, type, partition, Nsamp, datavector_by_Nmorphs
+#
+0  #  Do dataread for selectivity priors(0/1)
+# Yr, Seas, Fleet,  Age/Size,  Bin,  selex_prior,  prior_sd
+# feature not yet implemented
+#
+999
+
diff --git a/model_recipes/test/forecast.ss b/model_recipes/test/forecast.ss
new file mode 100644
index 0000000..5db7ba5
--- /dev/null
+++ b/model_recipes/test/forecast.ss
@@ -0,0 +1,64 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#C forecast file written by R function SS_writeforecast
+#C rerun model to get more complete formatting in forecast.ss_new
+#C should work with SS version: 3.3
+#C file write time: 2022-06-15 15:53:34
+# for all year entries except rebuilder; enter either: actual year, -999 for styr, 0 for endyr, neg number for rel. endyr
+1 # Benchmarks: 0=skip; 1=calc F_spr,F_btgt,F_msy; 2=calc F_spr,F0.1,F_msy; 3=add F_Blimit; 
+2 # Do_MSY: 1= set to F(SPR); 2=calc F(MSY); 3=set to F(Btgt) or F0.1; 4=set to F(endyr); 5=calc F(MEY) with MSY_unit options
+# if Do_MSY=5, enter MSY_Units; then list fleet_ID, cost/F, price/mt, include_in_Fmey_scaling; # -fleet_ID to fill; -9999 to terminate
+0.5 # SPR target (e.g. 0.40)
+0.4 # Biomass target (e.g. 0.40)
+#_Bmark_years: beg_bio, end_bio, beg_selex, end_selex, beg_relF, end_relF, beg_recr_dist, end_recr_dist, beg_SRparm, end_SRparm (enter actual year, or values of 0 or -integer to be rel. endyr)
+ 0 0 0 0 0 0 0 0 0 0
+#  2022 2022 2022 2022 2022 2022 2022 2022 2022 2022
+# value <0 convert to endyr-value; except -999 converts to start_yr; must be >=start_yr and <=endyr
+1 #Bmark_relF_Basis: 1 = use year range; 2 = set relF same as forecast below
+#
+1 # Forecast: -1=none; 0=simple_1yr; 1=F(SPR); 2=F(MSY) 3=F(Btgt) or F0.1; 4=Ave F (uses first-last relF yrs); 5=input annual F scalar
+# where none and simple require no input after this line; simple sets forecast F same as end year F
+10 # N forecast years 
+0.2 # Fmult (only used for Do_Forecast==5) such that apical_F(f)=Fmult*relF(f)
+#_Fcast_years:  beg_selex, end_selex, beg_relF, end_relF, beg_mean recruits, end_recruits  (enter actual year, or values of 0 or -integer to be rel. endyr)
+ 0 0 0 0 0 0
+#  2022 2022 2022 2022 2022 2022
+0 # Forecast selectivity (0=fcast selex is mean from year range; 1=fcast selectivity from annual time-vary parms)
+1 # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch) 
+# values for top, bottom and buffer exist, but not used when Policy=0
+0.4 # Control rule inflection for constant F (as frac of Bzero, e.g. 0.40); must be > control rule cutoff, or set to -1 to use Bmsy/SSB_unf 
+0.1 # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) 
+0.75 # Buffer:  enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax 
+3 #_N forecast loops (1=OFL only; 2=ABC; 3=get F from forecast ABC catch with allocations applied)
+3 #_First forecast loop with stochastic recruitment
+1 #_Forecast recruitment:  0= spawn_recr; 1=value*spawn_recr_fxn; 2=value*VirginRecr; 3=recent mean from yr range above (need to set phase to -1 in control to get constant recruitment in MCMC)
+1 # value is multiplier of SRR 
+0 #_Forecast loop control #5 (reserved for future bells&whistles) 
+2050  #FirstYear for caps and allocations (should be after years with fixed inputs) 
+0 # stddev of log(realized catch/target catch) in forecast (set value>0.0 to cause active impl_error)
+0 # Do West Coast gfish rebuilder output: 0=no; 1=yes 
+1999 # Rebuilder:  first year catch could have been set to zero (Ydecl)(-1 to set to 1999)
+2002 # Rebuilder:  year for current age structure (Yinit) (-1 to set to endyear+1)
+1 # fleet relative F:  1=use first-last alloc year; 2=read seas, fleet, alloc list below
+# Note that fleet allocation is used directly as average F if Do_Forecast=4 
+2 # basis for fcast catch tuning and for fcast catch caps and allocation  (2=deadbio; 3=retainbio; 5=deadnum; 6=retainnum); NOTE: same units for all fleets
+# Conditional input if relative F choice = 2
+# enter list of:  season,  fleet, relF; if used, terminate with season=-9999
+# 1 1 1
+# -9999 0 0  # terminator for list of relF
+# enter list of: fleet number, max annual catch for fleets with a max; terminate with fleet=-9999
+-9999 -1
+# enter list of area ID and max annual catch; terminate with area=-9999
+-9999 -1
+# enter list of fleet number and allocation group assignment, if any; terminate with fleet=-9999
+-9999 -1
+#_if N allocation groups >0, list year, allocation fraction for each group 
+# list sequentially because read values fill to end of N forecast
+# terminate with -9999 in year field 
+# no allocation groups
+#
+2 # basis for input Fcast catch: -1=read basis with each obs; 2=dead catch; 3=retained catch; 99=input apical_F; NOTE: bio vs num based on fleet's catchunits
+#enter list of Fcast catches or Fa; terminate with line having year=-9999
+#_Yr Seas Fleet Catch(or_F)
+-9999 1 1 0 
+#
+999 # verify end of input 
diff --git a/model_recipes/test/starter.ss b/model_recipes/test/starter.ss
new file mode 100644
index 0000000..f75b160
--- /dev/null
+++ b/model_recipes/test/starter.ss
@@ -0,0 +1,47 @@
+#V3.30.21.00;_safe;_compile_date:_Feb 10 2023;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.1
+#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.
+#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.
+#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov
+#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis
+#_Source_code_at:_https://github.com/nmfs-stock-synthesis/stock-synthesis
+
+#C starter file written by R function SS_writestarter
+#C rerun model to get more complete formatting in starter.ss_new
+#C should work with SS version: 3.30
+#C file write time: 2022-06-15 15:53:34
+data.ss
+control.ss
+0 # 0=use init values in control file; 1=use ss.par
+0 # run display detail (0,1,2)
+1 # detailed output (0=minimal for data-limited, 1=high (w/ wtatage.ss_new), 2=brief, 3=custom) 
+# custom report options: -100 to start with minimal; -101 to start with all; -number to remove, +number to add, -999 to end
+0 # write 1st iteration details to echoinput.sso file (0,1) 
+0 # write parm values to ParmTrace.sso (0=no,1=good,active; 2=good,all; 3=every_iter,all_parms; 4=every,active)
+0 # write to cumreport.sso (0=no,1=like&timeseries; 2=add survey fits)
+1 # Include prior_like for non-estimated parameters (0,1) 
+1 # Use Soft Boundaries to aid convergence (0,1) (recommended)
+#
+3 # Number of datafiles to produce:  0 turns off all *.ss_new; 1st is data_echo.ss_new, 2nd is data_expval.ss, 3rd and higher are data_boot_**N.ss,
+10 # Turn off estimation for parameters entering after this phase
+#
+0 # MCeval burn interval
+1 # MCeval thin interval
+0 # jitter initial parm value by this fraction
+-1 # min yr for sdreport outputs (-1 for styr); #_2009
+-2 # max yr for sdreport outputs (-1 for endyr+1; -2 for endyr+Nforecastyrs); #_2032
+0 # N individual STD years 
+#vector of year values 
+
+0.0001 # final convergence criteria (e.g. 1.0e-04) 
+0 # retrospective year relative to end year (e.g. -4)
+1 # min age for calc of summary biomass
+1 # Depletion basis:  denom is: 0=skip; 1=X*SPBvirgin; 2=X*SPBmsy; 3=X*SPB_styr; 4=X*SPB_endyr; 5=X*dyn_Bzero;  values>=11 invoke N multiyr (up to 9!) with 10's digit; >100 invokes log(ratio)
+1 # Fraction (X) for Depletion denominator (e.g. 0.4)
+1 # SPR_report_basis:  0=skip; 1=(1-SPR)/(1-SPR_tgt); 2=(1-SPR)/(1-SPR_MSY); 3=(1-SPR)/(1-SPR_Btarget); 4=rawSPR
+3 # F_reporting_units: 0=skip; 1=exploitation(Bio); 2=exploitation(Num); 3=sum(Apical_F's); 4=true F for range of ages; 5=unweighted avg. F for range of ages
+#COND 10 15 #_min and max age over which average F will be calculated with F_reporting=4 or 5
+0 # F_std_basis: 0=raw_annual_F; 1=F/Fspr; 2=F/Fmsy; 3=F/Fbtgt; where F means annual_F; values >=11 invoke N multiyr (up to 9!) with 10's digit; >100 invokes log(ratio)
+0 # MCMC output detail: integer part (0=default; 1=adds obj func components; 2= +write_report_for_each_mceval); and decimal part (added to SR_LN(R0) on first call to mcmc)
+0 # ALK tolerance ***disabled in code (example 0.0001)
+-1 # random number seed for bootstrap data (-1 to use long(time) as seed): # 1677620645
+3.30 # check value for end of file and for version control
diff --git a/tests/testthat.R b/tests/testthat.R
index 3e21509..54013a7 100644
--- a/tests/testthat.R
+++ b/tests/testthat.R
@@ -1,6 +1,25 @@
 library(testthat)
 library(ss3diags)
+library(r4ss)
 
-test_example_path <- system.file("data", package = "ss3diags")
+test_example_path <- system.file("extdata", "simple_small", package = "r4ss")
+files_path <- system.file("extdata", package = "ss3diags")
+run_tmp <- file.path(tempdir(check = TRUE), "test-runs")
+dir.create(run_tmp, showWarnings = FALSE)
+file.copy(from = list.files(files_path, full.names = TRUE), to = run_tmp)
+
+## Run retrospectives
+if (file.exists(file.path(files_path, "ss")) | file.exists(file.path(files_path, "ss.exe"))) {
+  r4ss::retro(dir = run_tmp, oldsubdir = "", newsubdir = "retrospectives", years = 0:-3, show_in_console = FALSE)
+  # Creating retrospective object here so that multiple test files can access it without having to re-run retrospective for each test
+  retroModels <- r4ss::SSgetoutput(
+    dirvec = file.path(
+      run_tmp, "retrospectives",
+      paste0("retro", 0:-3)
+    )
+  )
+  retrosum.simple <- r4ss::SSsummarize(retroModels)
+}
+on.exit(unlink(tmp_path, recursive = TRUE))
 
 test_check("ss3diags")
diff --git a/tests/testthat/SSmase_test.R b/tests/testthat/SSmase_test.R
index 2425bd8..893b116 100644
--- a/tests/testthat/SSmase_test.R
+++ b/tests/testthat/SSmase_test.R
@@ -1,29 +1,31 @@
-library(ss3diags)
+## File is deprecated by MOshima 6/12/23
+
+# library(ss3diags)
 # Tests
 # sma = r4ss::SSsummarize(retro.sma)
 # Check joint MASE for indices
-SSmase(retro.sma)
+# SSmase(retro.sma)
 # select indices
-SSmase(retro.sma, indexselect = c(4))
+# SSmase(retro.sma, indexselect = c(4))
 # check for length comps
 # smaL = ss3diags::SSretroComps(retro.sma)
-SSmase(retro.sma, quants = "len", indexselect = 2:4)
+# SSmase(retro.sma, quants = "len", indexselect = 2:4)
 
 # check manually
-mase <- SSmase(retro.sma, residuals = T)
+# mase <- SSmase(retro.sma, residuals = T)
 # Joint mase from index 1, 3,4
-mae <- apply(
-  mase$Residuals[c("Pred.Res", "Naive.Res")], 2,
-  function(x) mean(abs(x))
-)
+# mae <- apply(
+#  mase$Residuals[c("Pred.Res", "Naive.Res")], 2,
+#  function(x) mean(abs(x))
+# )
 # MASE
-mae[1] / mae[2]
+# mae[1] / mae[2]
 # compare
-SSmase(retro.sma)
+# SSmase(retro.sma)
 
 # Check hake
 # phk= r4ss::SSsummarize(retro.phk)
-SSmase(retro.phk, verbose = T, indexselect = 1)
+# SSmase(retro.phk, verbose = T, indexselect = 1)
 # check for age comps
-phkA <- ss3diags::SSretroComps(retro.phk)
-SSmase(retro.phk, quants = "age")
+# phkA <- ss3diags::SSretroComps(retro.phk)
+# SSmase(retro.phk, quants = "age")
diff --git a/tests/testthat/test-HCxval.R b/tests/testthat/test-HCxval.R
index 5e36b12..d86dc18 100644
--- a/tests/testthat/test-HCxval.R
+++ b/tests/testthat/test-HCxval.R
@@ -1,7 +1,7 @@
 ## Test script for hindcast cross validation and MASE
 
-retroSimple <- ss3diags::retroSimple
-retrosum.simple <- r4ss::SSsummarize(retroSimple)
+# retroSimple <- ss3diags::retroSimple
+# retrosum.simple <- r4ss::SSsummarize(retroSimple)
 
 path <- file.path(tempdir(), "test_runs")
 dir.create(path, showWarnings = FALSE)
@@ -9,6 +9,10 @@ dir.create(path, showWarnings = FALSE)
 ## Hindcast plotting of indices
 
 test_that("Hindcast plot is created for simple model indices", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   SSplotHCxval(retrosum.simple,
     add = T,
     verbose = F,
@@ -37,15 +41,24 @@ test_that("Hindcast plot is created for simple model indices", {
 #### note: function was too complicated to replicate in the test script so used values directly from running the function. If code for calculations changes, the values will change and it will error or if the SS3 input files change, the values will be different and it will error.
 
 test_that("MASE table gives expected values for simple model", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
+
   mase <- SSplotHCxval(retrosum.simple, add = T, verbose = F)
 
-  expect_match(mase$Index[1], "Survey")
-  expect_equal(round(mase$MASE[1], 7), 0.5836339)
-  expect_equal(round(mase$MAE.PR[1], 7), 0.1856207)
+  expect_match(mase$Index[1], "SURVEY1")
+  expect_equal(round(mase$MASE[1], 7), 0.494371)
+  expect_equal(round(mase$MAE.PR[1], 7), 0.0873324)
 })
 
 
 test_that("SSretroComps returns the correct comp data for simple model", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   retro_comps <- SSretroComps(retroSimple)
 
   expect_equal(retro_comps$n, 6)
@@ -58,7 +71,11 @@ test_that("SSretroComps returns the correct comp data for simple model", {
 
 
 test_that("SSmase base.adj changes", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   ssmase <- SSmase(retrosum.simple, MAE.base.adj = 0.15)
 
-  expect_equal(round(ssmase$MASE.adj[1], 7), 0.8503379)
+  expect_equal(round(ssmase$MASE.adj[1], 7), 0.4943710)
 })
diff --git a/tests/testthat/test-plotJABBA.R b/tests/testthat/test-plotJABBA.R
index 8dde672..78644f3 100644
--- a/tests/testthat/test-plotJABBA.R
+++ b/tests/testthat/test-plotJABBA.R
@@ -1,13 +1,9 @@
 ## Test script for JABBA residual plot
-
-simple <- ss3diags::simple
-
+simple <- r4ss::SS_output(dir = test_example_path, verbose = FALSE, printstats = FALSE)
 
 path <- file.path(tempdir(), "test_runs")
 dir.create(path, showWarnings = FALSE)
 
-
-
 ## Simple
 test_that("file of simple_cpue_jabbaresiduals plot exists", {
   SSplotJABBAres(simple,
@@ -26,7 +22,6 @@ test_that("file of simple_len_jabbaresiduals plot exists", {
     use_png = TRUE,
     print_plot = TRUE,
     subplots = "len",
-    # indexselect = 2,
     plotdir = path,
     filenameprefix = "simple_len2_"
   )
@@ -35,29 +30,29 @@ test_that("file of simple_len_jabbaresiduals plot exists", {
 })
 
 ## CAAL uncomment when con option is finished in function
-test_that("file of simple_con_jabbaresiduals plot exists", {
-  #
-  SSplotJABBAres(simple,
-    use_png = TRUE,
-    print_plot = TRUE,
-    subplots = "con",
-    plotdir = path,
-    filenameprefix = "simple_con_"
-  )
-
-  expect_true(file.exists(file.path(path, "simple_con_jabbaresidual.png")))
-  #
-})
+# test_that("file of simple_con_jabbaresiduals plot exists", {
+#
+# SSplotJABBAres(simple,
+# use_png = TRUE,
+# print_plot = TRUE,
+# subplots = "con",
+# plotdir = path,
+# filenameprefix = "simple_con_"
+# )
+
+# expect_true(file.exists(file.path(path, "simple_con_jabbaresidual.png")))
+#
+# })
 
 test_that("calculate RMSE for CPUE index", {
   rmse <- SSrmse(simple, quants = "cpue")
 
-  expect_equal(rmse$RMSE$RMSE.perc[1], 20.4)
+  expect_equal(rmse$RMSE$RMSE.perc[1], 21.6)
 })
 
 
 test_that("Calculate combined RMSE for length comp data", {
   rmse <- SSrmse(simple, quants = "len")$RMSE
 
-  expect_equal(rmse$RMSE.perc[3], 4.3)
+  expect_equal(rmse$RMSE.perc[3], 2.5)
 })
diff --git a/tests/testthat/test-retro.R b/tests/testthat/test-retro.R
index e00ec9f..0272a80 100644
--- a/tests/testthat/test-retro.R
+++ b/tests/testthat/test-retro.R
@@ -3,11 +3,13 @@
 path <- file.path(tempdir(), "test_runs")
 dir.create(path, showWarnings = FALSE)
 
-retroSimple <- ss3diags::retroSimple
-retrosum.simple <- r4ss::SSsummarize(retroSimple)
 
 ## SSB
 test_that("Retrospective plot is created for SSB (w/ deprecated parameters)", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   SSplotRetro(retrosum.simple,
     subplots = "SSB",
     use_png = TRUE,
@@ -21,6 +23,10 @@ test_that("Retrospective plot is created for SSB (w/ deprecated parameters)", {
 
 ## F
 test_that("Retrospective plot is created for F (w/ deprecated parameters)", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   SSplotRetro(retrosum.simple,
     subplots = "F",
     use_png = TRUE,
@@ -35,6 +41,10 @@ test_that("Retrospective plot is created for F (w/ deprecated parameters)", {
 ## Mohn's Rho table using SShcbias function. It's the same code used in SSplotRetro to make the table.
 ## SSB
 test_that("Mohn's Rho table is created for SSB", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   n <- retrosum.simple$n
   endyrs <- retrosum.simple$endyrs - seq(0, n - 1, 1)
   exp <- retrosum.simple$SpawnBio
@@ -60,6 +70,10 @@ test_that("Mohn's Rho table is created for SSB", {
 
 ## F
 test_that("Mohn's Rho table is created for F", {
+  skip_if(
+    !exists("retrosum.simple"),
+    message = "skipping test that requires SS3 executable"
+  )
   n <- retrosum.simple$n
   endyrs <- retrosum.simple$endyrs - seq(0, n - 1, 1)
   exp <- retrosum.simple$Fvalue
diff --git a/tests/testthat/test-runs-test.R b/tests/testthat/test-runs-test.R
index b8d1323..094ad92 100644
--- a/tests/testthat/test-runs-test.R
+++ b/tests/testthat/test-runs-test.R
@@ -1,4 +1,4 @@
-simple <- ss3diags::simple
+simple <- r4ss::SS_output(dir = test_example_path, verbose = FALSE, printstats = FALSE)
 
 #### Simple Model ###############################################################
 
@@ -7,7 +7,7 @@ simple <- ss3diags::simple
 
 ## SSrunstest function
 test_that("runs test works with simple model", {
-  test.resids <- simple$cpue[which(simple$cpue$Fleet_name == "Survey"), c("Fleet_name", "Yr", "Obs", "Exp")]
+  test.resids <- simple$cpue[which(simple$cpue$Fleet_name == "SURVEY2"), c("Fleet_name", "Yr", "Obs", "Exp")]
   test.resids$residuals <- log(test.resids$Obs) - log(test.resids$Exp)
 
   mu <- 0
@@ -30,13 +30,13 @@ test_that("runs test works with simple model", {
   run_cpue <- SSrunstest(simple, quants = "cpue")
 
   ## testing structure of dataframe
-  expect_match(run_cpue$Index[1], "Survey")
+  expect_match(run_cpue$Index[1], "SURVEY1")
   expect_equal(nrow(run_cpue), n.cpue)
   ## testing values in the first row
-  expect_match(run_cpue$test[1], "Failed")
-  expect_equal(run_cpue$runs.p[1], test.p)
-  expect_equal(run_cpue$sigma3.lo[1], lcl)
-  expect_equal(run_cpue$sigma3.hi[1], ucl)
+  expect_match(run_cpue$test[2], "Passed")
+  expect_equal(run_cpue$runs.p[2], test.p)
+  expect_equal(run_cpue$sigma3.lo[2], lcl)
+  expect_equal(run_cpue$sigma3.hi[2], ucl)
 
   ## for length comp
   len.test.resids <- simple$lendbase[which(simple$lendbase$Fleet == 1), ]
@@ -116,7 +116,7 @@ test_that("runs test works with simple model", {
   run_fish <- SSrunstest(simple, quants = "len")
 
   ## testing structure of dataframe
-  expect_match(run_fish$Index[1], "Fishery")
+  expect_match(run_fish$Index[1], "FISHERY")
   expect_equal(nrow(run_fish), n.fish)
   ## testing values in first row
   expect_equal(run_fish$runs.p[1], test.p)
@@ -125,13 +125,13 @@ test_that("runs test works with simple model", {
 
   ## specifying cpue index
   run_cpue <- SSrunstest(simple, quants = "len", indexselect = 2)
-  expect_match(run_cpue$Index, "Survey")
+  expect_match(run_cpue$Index, "SURVEY1")
 
-  # CAAL
-  run_con <- SSrunstest(simple, quants = "con")
+  # CAAL since switching to simple_small, there is no conditional-age-at-length data.
+  # run_con <- SSrunstest(simple, quants = "con")
 
-  expect_match(run_con$Index, "Fishery")
-  expect_match(run_con$test, "Passed")
+  # expect_match(run_con$Index, "Fishery")
+  # expect_match(run_con$test, "Passed")
 })
 
 
@@ -148,26 +148,26 @@ dir.create(path, showWarnings = FALSE)
 ## Simple
 test_that("file of simple_cpue_residruns plot exists", {
   SSplotRunstest(simple,
-    png = TRUE,
-    print = T,
+    use_png = TRUE,
+    print_plot = T,
     subplots = "cpue",
     indexselect = 1,
     plotdir = path,
     filenameprefix = "simple_"
   )
 
-  expect_true(file.exists(file.path(path, "simple_residruns_Survey.png")))
+  expect_true(file.exists(file.path(path, "simple_residruns_SURVEY1.png")))
 })
 
 
-test_that("file of simple_con_residruns plot exists", {
+test_that("file of simple_len_residruns plot exists", {
   SSplotRunstest(simple,
-    png = TRUE,
-    print = T,
-    subplots = "con",
+    use_png = TRUE,
+    print_plot = T,
+    subplots = "len",
     plotdir = path,
-    filenameprefix = "simple_con_"
+    filenameprefix = "simple_len_"
   )
 
-  expect_true(file.exists(file.path(path, "simple_con_residruns_Fishery.png")))
+  expect_true(file.exists(file.path(path, "simple_len_residruns_FISHERY.png")))
 })
diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R
index 66dc0c1..f8670d7 100644
--- a/tests/testthat/test-uncertainty.R
+++ b/tests/testthat/test-uncertainty.R
@@ -1,18 +1,18 @@
 # Test script for evaluating model uncertainty
+simple <- r4ss::SS_output(dir = test_example_path, verbose = FALSE, printstats = FALSE)
 
-simple <- ss3diags::simple
 
 mvln.msy <- SSdeltaMVLN(simple, Fref = "MSY")
 
 test_that("mle table values are as expected for Simple with Fref = MSY", {
-  expect_equal(round(mvln.msy$mle$harvest[1], 7), 0.8845945)
-  expect_equal(round(mvln.msy$mle$stock[1], 7), 2.0921525)
-  expect_equal(round(mvln.msy$mle$F[1], 7), 0.1009570)
+  expect_equal(round(mvln.msy$mle$harvest[1], 6), 1.210293)
+  expect_equal(round(mvln.msy$mle$stock[1], 6), 1.0998)
+  expect_equal(round(mvln.msy$mle$F[1], 6), 0.239068)
 })
 
 
 test_that("Kobe plot values are correct", {
-  kb <- mvln.msy$kb[which(mvln.msy$kb$year == 100), ] # last year of model
+  kb <- mvln.msy$kb[which(mvln.msy$kb$year == 2022), ] # last year of model
 
   kb$green <- ifelse(kb$stock > 1 & kb$harvest < 1, 1, 0)
   kb$red <- ifelse(kb$stock < 1 & kb$harvest > 1, 1, 0)