diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 00000000..4ed9401a
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,67 @@
+name: FSPS Regression Tests
+
+on:
+  push:
+    branches: [ master ]
+  pull_request:
+    branches: [ master ]
+
+jobs:
+  test-matrix:
+    runs-on: ubuntu-latest
+    
+    env:
+      SPS_HOME: ${{ github.workspace }}
+      FSPS_TEST_RTOL: 1.0E-5  # Set global relative tolerance here
+
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          # Define the configurations to test.
+          # 'opts': The library options passed to test_runner.
+          # 'ref_suffix': The suffix of the reference file in tests/data/
+          
+          # 1. Standard (Default)
+          - name: MILES + MIST (default)
+            opts: "--isoc mist --spec miles"
+            ref_suffix: "MILES-1_MIST-1"
+
+          # 2. BaSeL+Padova
+          - name: BaSEL + Padova
+            opts: "--isoc pdva --spec basel"
+            ref_suffix: "MILES-0_MIST-0_BASEL-1_PADOVA-1"
+
+          # 3. Dust Model (THEMIS)
+          - name: THEMIS (dust)
+            opts: "--isoc mist --spec miles --dust themis"
+            ref_suffix: "MILES-1_MIST-1_THEMIS-1_DL07-0"
+
+          # 4. C3K
+          - name: C3K + MIST
+            opts: "--isoc mist --spec c3k_afe+0.0"
+            ref_suffix: "MILES-0_C3K-1_MIST-1"
+
+          # 5. BPASS (Different logic path)
+          - name: BPASS
+            opts: "--isoc bpss"
+            ref_suffix: "MIST-0_BPASS-1"
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
+    - name: Install GFortran
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y gfortran
+
+    - name: Compile FSPS and Test Runner
+      # We force FFLAGS to include the matrix flags + standard flags + silence underflows
+      run: |
+        make clean
+        make test FFLAGS="-cpp -O3 -fPIC -ffpe-summary=invalid,zero,overflow"
+
+    - name: Run Test
+      run: |
+        ./test_runner tests/data/sps_ref_${{ matrix.ref_suffix }}.bin ${{ matrix.opts }}
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 8cd11d43..a3d4c05f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,3 @@
-*.o
-*.mod
-*.exe
 .DS_Store
 nr
 SPECTRA/CKC14/*
@@ -11,3 +8,17 @@ metalml.f90
 colorteff.f90
 bytype.f90
 Makefile.prv
+
+# Ignore compiler outputs
+build/
+*.o
+*.mod
+*.exe
+
+# Executables produced by Makefile
+autosps
+simple
+lesssimple
+spec_bin
+generate_test_data
+test_runner
\ No newline at end of file
diff --git a/LICENSE b/LICENSE
index 99ffd933..beb57cb1 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2009-2021 Charlie Conroy & contributors.
+Copyright (c) 2009-2026 Charlie Conroy & contributors.
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/Makefile b/Makefile
new file mode 100644
index 00000000..f0ea587c
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,111 @@
+# ===================================
+# Compiler Configuration
+# ===================================
+
+# Default to gfortran, allow for an override (e.g. make FC=ifort)
+FC = gfortran
+
+# Default flags
+FFLAGS ?= -O3 -cpp -fPIC
+
+# Directory configuration
+SRC_DIR := src
+TEST_DIR := tests
+BUILD_DIR := build
+
+# Module directory flags:
+# Gfortran uses -J to specify where to put/find .mod files
+# Intel (ifort) uses -module. We default to Gfortran syntax here.
+# You can override this: make MOD_FLAG=-module
+MOD_FLAG ?= -J
+
+# Combine flags to include the build dir for .mod search
+FCFLAGS := $(FFLAGS) $(MOD_FLAG)$(BUILD_DIR) -I$(BUILD_DIR)
+
+# ===================================
+# Source and Object Definitions
+# ===================================
+
+# Tell make to look for source files in these directories
+VPATH = $(SRC_DIR):$(TEST_DIR)
+
+# The list of programs to build (executables)
+PROGS = simple lesssimple autosps spec_bin
+
+# The common object files required by the programs
+# We wrap them in addprefix to place them inside the build directory
+COMMON_NAMES = sps_vars.o sps_utils.o compsp.o csp_gen.o ssp_gen.o \
+    getmags.o locate.o funcint.o sps_setup.o pz_convol.o \
+    get_tuniv.o intsfwght.o imf.o imf_weight.o add_dust.o \
+    getspec.o sbf.o add_bs.o mod_hb.o add_remnants.o getindx.o \
+    smoothspec.o mod_gb.o add_nebular.o add_xrb.o write_isochrone.o \
+    sfhstat.o linterp.o tsum.o add_agb_dust.o linterparr.o \
+    ztinterp.o vacairconv.o igm_absorb.o get_lumdist.o attn_curve.o \
+    sfh_weight.o sfhlimit.o sfhinfo.o setup_tabular_sfh.o agn_dust.o
+
+COMMON_OBJS = $(addprefix $(BUILD_DIR)/, $(COMMON_NAMES))
+
+# ===================================
+# Rules
+# ===================================
+
+.PHONY: all clean test
+
+all: $(PROGS)
+
+# --- Compilation Rules ---
+
+# Ensure build directory exists before compiling
+$(BUILD_DIR):
+	@mkdir -p $(BUILD_DIR)
+
+# Pattern rule: Compile any .f90 found in VPATH to .o in BUILD_DIR
+$(BUILD_DIR)/%.o: %.f90 | $(BUILD_DIR)
+	$(FC) $(FCFLAGS) -c $< -o $@
+
+# Specific dependencies to enforce compilation order
+
+# sps_utils.o specifically depends on sps_vars.o
+$(BUILD_DIR)/sps_utils.o: $(BUILD_DIR)/sps_vars.o
+
+# All other common objects depend on both vars and utils.
+REST_OF_COMMON = $(filter-out $(BUILD_DIR)/sps_vars.o $(BUILD_DIR)/sps_utils.o, $(COMMON_OBJS))
+
+$(REST_OF_COMMON): $(BUILD_DIR)/sps_vars.o $(BUILD_DIR)/sps_utils.o
+
+# Main program objects also wait for modules
+$(BUILD_DIR)/simple.o $(BUILD_DIR)/lesssimple.o $(BUILD_DIR)/autosps.o $(BUILD_DIR)/spec_bin.o: $(BUILD_DIR)/sps_vars.o $(BUILD_DIR)/sps_utils.o
+
+# Dependencies for test objects
+$(BUILD_DIR)/generate_test_data.o: $(BUILD_DIR)/sps_vars.o $(BUILD_DIR)/sps_utils.o
+$(BUILD_DIR)/test_runner.o: $(BUILD_DIR)/sps_vars.o $(BUILD_DIR)/sps_utils.o
+
+# --- Linking Rules ---
+
+autosps: $(BUILD_DIR)/autosps.o $(COMMON_OBJS)
+	$(FC) $(FCFLAGS) -o $@ $^
+
+simple: $(BUILD_DIR)/simple.o $(COMMON_OBJS)
+	$(FC) $(FCFLAGS) -o $@ $^
+
+lesssimple: $(BUILD_DIR)/lesssimple.o $(COMMON_OBJS)
+	$(FC) $(FCFLAGS) -o $@ $^
+
+spec_bin: $(BUILD_DIR)/spec_bin.o $(BUILD_DIR)/sps_vars.o
+	$(FC) $(FCFLAGS) -o $@ $^
+
+# --- Test Targets ---
+
+generate_test_data: $(BUILD_DIR)/generate_test_data.o $(COMMON_OBJS)
+	$(FC) $(FCFLAGS) -o $@ $^
+
+test_runner: $(BUILD_DIR)/test_runner.o $(COMMON_OBJS)
+	$(FC) $(FCFLAGS) -o $@ $^
+
+test: test_runner
+
+# --- Utilities ---
+
+clean:
+	rm -rf $(BUILD_DIR) $(PROGS) generate_test_data test_runner
+
diff --git a/README.md b/README.md
index 93c4c668..ea55ea20 100644
--- a/README.md
+++ b/README.md
@@ -1,23 +1,24 @@
 FSPS: Flexible Stellar Population Synthesis
 =====
-version 3.2
+![Version Badge](https://img.shields.io/badge/version-v3.2-blue)
 
 References
 ---------
 When using this code please cite the following papers:
- * Conroy, Gunn, & White 2009, ApJ, 699, 486
- * Conroy & Gunn 2010, ApJ, 712, 833
+ * [Conroy, Gunn, & White 2009, ApJ, 699, 486](https://ui.adsabs.harvard.edu/abs/2009ApJ...699..486C)
+ * [Conroy & Gunn 2010, ApJ, 712, 833](https://ui.adsabs.harvard.edu/abs/2010ApJ...712..833C)
 
 Installation
 ----------
-If you have git installed, FSPS can be obtained with the following commands:
-```
+If you have Git installed, FSPS can be obtained with the following commands:
+
+```sh
 cd /path/to/desired/location/
 git clone https://github.com/cconroy20/fsps
 ```
-Otherwise download a gzipped tarball from [here](https://github.com/cconroy20/fsps/releases). Then follow the instructions at [doc/INSTALL](doc/INSTALL).
+Otherwise download a gzipped tarball from [here](https://github.com/cconroy20/fsps/releases). Then follow the instructions at [`doc/INSTALL`](doc/INSTALL).
 
-You should not need to update the git repository until an update is announced (which is why you need to be on the mailing list - see [doc/INSTALL](doc/INSTALL)).  If you've obtained FSPS using git then when an update is announced you will need to simply type ``cd $SPS_HOME; git pull`` and then recompile.  If you have made your own edits to the FSPS files, git will attempt to gracefully merge your local version with the repository version.
+You should not need to update the Git repository until an update is announced (which is why you need to be on the mailing list - see [`doc/INSTALL`](doc/INSTALL)).  If you've obtained FSPS using Git then when an update is announced you will need to simply type `cd $SPS_HOME; git pull` and then recompile (just type `make`).  If you have made your own edits to the FSPS files, Git will attempt to gracefully merge your local version with the repository version.
 
 Documentation
 ------
@@ -28,17 +29,17 @@ Contents
 Below is a brief description of the contents of the directories in the
 fsps root directory:
 
- * `ISOCHRONES`: Contains the isochrone tables for the BaSTI and Padova
+ * [`ISOCHRONES`](ISOCHRONES): Contains the isochrone tables for the BaSTI and Padova
 isochrone sets.  The Geneva isochrones have been pasted onto the BaSTI
 and Padova tables for high masses (M>70Msun), and the low-mass Lyon
 models have been pasted on at low masses.  You should not edit these
 files unless you know what you're doing.
 
- * `OUTPUTS`: Contains the outputs of a few example calls of the routines
+ * [`OUTPUTS`](OUTPUTS): Contains the outputs of a few example calls of the routines
 autosps and simple.  You may wish to use this directory for all
 outputs of the fsps routines.
 
- * `SPECTRA`: Contains the spectral libraries, the spectrum of an A0V star
+ * [`SPECTRA`](SPECTRA): Contains the spectral libraries, the spectrum of an A0V star
 used to set the Vega magnitude zero points, and a spectrum of the Sun.
 The BaSeL spectra (based on the Kurucz models) are in binary format,
 primarily to make the read in time faster and to decrease the size of
@@ -46,28 +47,28 @@ the fsps download.  The Hot_spectra directory contain the libraries
 for O stars, WR stars, and post-AGB stars, from Smith et al. 2002 and
 Rauch 2003, respectively.
 
- * `data`: Contains files that define the set of filters and indices used
+ * `build`: Contains the compiled object files (`.o`) and Fortran modules (`.mod`).
+This directory is automatically created when you run `make`.
+
+ * [`data`](data): Contains files that define the set of filters and indices used
 in FSPS and the tabulated imfs and sfhs if those options are set.  The
 files in this directory are readily user editable.
 
- * `doc`: Contains the manual, revision history, and installation
+ * [`doc`](doc): Contains the manual, revision history, and installation
 instructions.
 
-* `dust`: Contains the dust attenuation curves for the Witt & Gordon
+* [`dust`](dust): Contains the dust attenuation curves for the Witt & Gordon
 (2000) dust model and the dust emission spectra from the Draine & Li
 2007 grain model.  Also contains the circumstellar dust models from 
 Villaume et al. 2015 and the AGN dusty torus models of Nenkova et al. 2008.
 
-* `nebular`: Contains the Cloudy lookup tables for nebular emission 
+* [`nebular`](nebular): Contains the Cloudy lookup tables for nebular emission 
 (both continuum and line emission) computed by Nell Byler.
 
-* `pro`: Contains IDL files for reading in the .mag, .indx, and .spec
+* [`pro`](pro): Contains IDL files for reading in the .mag, .indx, and .spec
 output files
 
-* `src`: Contains the source files and routines from Numerical Recipes.
-
-
-
-
-
+* [`src`](src): Contains the source files and routines from Numerical Recipes.
 
+* [`tests`](tests): Contains the regression test suite and scripts for generating
+reference comparison data.
diff --git a/doc/INSTALL b/doc/INSTALL
index c779aefe..df0f8bf7 100644
--- a/doc/INSTALL
+++ b/doc/INSTALL
@@ -1,18 +1,22 @@
 SETUP INSTRUCTIONS
 
-1.   edit your .cshrc (or .bashrc) file to something like:
-     	  setenv SPS_HOME /home/user/fsps/
-	  (this should point to the directory containing the src directory)
+1.  Edit your .cshrc (or .bashrc/zshrc) file to set the SPS_HOME variable:
+            export SPS_HOME="/home/user/fsps/"
+       (this must point to the root directory containing the src/ directory)
 
-2.  edit the makefile in the src directory so that it is using your
-    fortran compiler one should then be able to compile everything by
-    simply typing "make"
+2.  Compile the code. The Makefile is located in the root directory.
+    By default, it uses gfortran. You can compile by simply typing:
+            make
+
+    If you wish to use a different compiler (e.g., ifort), you can specify it
+    on the command line without editing the Makefile:
+            make FC=ifort
 
 3.  email cconroy@cfa.harvard.edu so that I can keep you posted on updates.  
       -> this step is actually important!  Every code has bugs.  If I
-	 don't know that you're using this code, I can't tell you
-	 that I found a catastrophic bug that invalidates all of your
-	 results.
+      don't know that you're using this code, I can't tell you
+      that I found a catastrophic bug that invalidates all of your
+      results.
 
 4.  please email me if have any problems or find any bugs.
 
@@ -23,4 +27,10 @@ NOTES:
 1.  It has been reported that the code crashes if compiled with
     earlier versions of gfortran (specifically v4.2.1).  The code has
     been tested and compiles sucessfully with gfortran v4.4 and later.
-    (FYI, I currently compile FSPS with gfortran v6.4.0).
\ No newline at end of file
+    (FYI, I currently compile FSPS with gfortran v6.4.0).
+
+2.  As of 2026, the code allows for runtime selection of isochrones and
+    spectral libraries. You do NOT need to edit sps_vars.f90 and 
+    recompile to switch between MIST, Padova, MILES, BaSeL, etc.
+    These options can now be passed directly to the setup routines or 
+    selected interactively in the `autosps` driver.
\ No newline at end of file
diff --git a/doc/REVISION_HISTORY b/doc/REVISION_HISTORY
index 825c5ed3..2b440e89 100644
--- a/doc/REVISION_HISTORY
+++ b/doc/REVISION_HISTORY
@@ -149,4 +149,18 @@
 
 04/2/21
  - Release of v3.2
- - added THEMIS dust emission model
\ No newline at end of file
+ - added THEMIS dust emission model
+
+01/28/26
+ - Major modernization of the codebase and build system.
+ - Key Feature: Isochrone and spectral libraries can now be selected at 
+   runtime. You no longer need to edit `sps_vars.f90` and recompile to switch 
+   between MIST, Padova, MILES, BaSeL, etc.
+ - Implemented dynamic memory allocation throughout the driver programs 
+   (`simple`, `lesssimple`, `autosps`), removing compile-time static array limits.
+ - `autosps` is now interactive and will prompt the user to select their desired
+   libraries at runtime.
+ - The `Makefile` has been moved to the root directory and now uses an 
+   out-of-source build system (all artifacts are output to `build/`).
+ - Fixed legacy compiler warnings regarding FORMAT strings.
+ - Updated CI/CD infrastructure for regression testing.
\ No newline at end of file
diff --git a/src/Makefile b/src/Makefile
deleted file mode 100644
index d783a455..00000000
--- a/src/Makefile
+++ /dev/null
@@ -1,75 +0,0 @@
-# NB: this code is fairly memory intensive.  If you run into memory errors
-# when building this code, you may have to increase the default memory 
-# allocation.
-
-#======================
-# Fortran 90 compiler
-# (Uncomment only one)
-#======================
-# GNU
-F90 = gfortran
-#---------------------
-# Intel
-#F90 = ifort
-#---------------------
-# PGI
-#F90 = pgf90
-# PGI (franklin.nersc.gov)
-#F90 = ftn
-#---------------------
-# Pathscale
-#F90 = pathf90
-#---------------------
-# IBM
-#F90 = xlf90_r
-#---------------------
-
-#========================
-# Compiler Optimizations
-#  (Uncomment only one)
-#========================
-F90FLAGS := -O3 -cpp -fPIC $(F90FLAGS) # prepend so a user can set -Ofast from outside the Makefile
-#---------------------
-# Intel (Itanium 2 processor)
-#F90FLAGS = -O3 -mcpu=itanium2 -cpp
-# Intel (EM64T/AMD64 processor)
-#F90FLAGS = -O3 -funroll-loops -xP -cpp
-#---------------------
-
-PROGS = simple lesssimple autosps spec_bin
-
-COMMON = sps_vars.o sps_utils.o compsp.o csp_gen.o ssp_gen.o \
-	getmags.o locate.o funcint.o sps_setup.o pz_convol.o \
-	get_tuniv.o intsfwght.o imf.o imf_weight.o add_dust.o \
-	getspec.o sbf.o add_bs.o mod_hb.o add_remnants.o getindx.o \
-	smoothspec.o mod_gb.o add_nebular.o add_xrb.o write_isochrone.o \
-	sfhstat.o linterp.o tsum.o add_agb_dust.o linterparr.o \
-	ztinterp.o vacairconv.o igm_absorb.o get_lumdist.o attn_curve.o \
-	sfh_weight.o sfhlimit.o sfhinfo.o setup_tabular_sfh.o agn_dust.o
-
-all : $(PROGS)
-
-clean :
-	rm -rf *.o *.mod *.MOD *~
-
-autosps : autosps.o $(COMMON)
-	$(F90) $(F90FLAGS) -o autosps.exe autosps.o $(COMMON)
-
-simple : simple.o $(COMMON)
-	$(F90) $(F90FLAGS) -o simple.exe simple.o $(COMMON)
-
-lesssimple : lesssimple.o $(COMMON)
-	$(F90) $(F90FLAGS) -o lesssimple.exe lesssimple.o $(COMMON)
-
-spec_bin : spec_bin.o sps_vars.o
-	$(F90) $(F90FLAGS) -o spec_bin.exe spec_bin.o sps_vars.o
-
-sps_vars.o : sps_vars.f90
-	$(F90) $(F90FLAGS) sps_vars.f90 -c
-
-sps_utils.o : sps_utils.f90 sps_vars.o
-	$(F90) $(F90FLAGS) sps_utils.f90 -c
-
-%.o : %.f90 sps_vars.o sps_utils.o
-	$(F90) $(F90FLAGS) -o $@ -c $<
-
diff --git a/src/autosps.f90 b/src/autosps.f90
index 210563fc..8614e0f2 100644
--- a/src/autosps.f90
+++ b/src/autosps.f90
@@ -6,23 +6,59 @@ PROGRAM AUTOSPS
   IMPLICIT NONE
 
   INTEGER :: z
-  REAL(SP), DIMENSION(ntfull,nspec)  :: spec_ssp
-  REAL(SP), DIMENSION(ntfull)        :: mass_ssp,lbol_ssp
+
+  REAL(SP), ALLOCATABLE :: spec_ssp(:,:)
+  REAL(SP), ALLOCATABLE :: mass_ssp(:),lbol_ssp(:)
+  TYPE(COMPSPOUT), ALLOCATABLE :: ocompsp(:)
+
   CHARACTER(100) :: file1='',aux
   CHARACTER(3)  :: str
   TYPE(PARAMS)  :: pset
-  TYPE(COMPSPOUT), DIMENSION(ntfull) :: ocompsp
+
+  ! Variables for library selection
+  CHARACTER(10) :: iso_in, spec_in
 
   !---------------------------------------------------------------!
   !---------------------------------------------------------------!
   
-  IF (isoc_type.NE.'pdva') THEN
-     WRITE(*,*) 'ERROR: autosps only works with the "Padova" isochrones'
-     WRITE(*,*) '       edit sps_vars.f90 to turn these isochrones on'
-     RETURN
+  WRITE(6,*) '=================================================='
+  WRITE(6,*) '             FSPS INTERACTIVE MODE                '
+  WRITE(6,*) '=================================================='
+
+  ! --- Ask for libraries at runtime ---
+  WRITE(6,*) 'Choose Isochrones [def: mist]:'
+  WRITE(6,*) '(mist, pdva, parsec, basti, bpass, etc.)'
+  READ(5,'(A)') aux
+  IF (LEN_TRIM(aux).EQ.0) THEN
+     iso_in = 'mist'
+  ELSE
+     READ(aux,'(A)') iso_in
+  ENDIF
+  
+  WRITE(6,*) 'Choose Spectral Library [def: miles]:'
+  WRITE(6,*) '(miles, basel, bpass, c3k_afe+0.0, etc.)'
+  READ(5,'(A)') aux
+  IF (LEN_TRIM(aux).EQ.0) THEN
+     spec_in = 'miles'
+  ELSE
+     READ(aux,'(A)') spec_in
   ENDIF
 
+  ! --- Initialize Environment Immediately ---
+  WRITE(6,*) 'Initializing libraries...'
+  ! We load ALL metallicities (-1) so we can query nz and zlegend
+  CALL SPS_SETUP(-1, TRIM(iso_in), TRIM(spec_in))
+
+  ! --- Allocate Memory ---
+  IF (.NOT. ALLOCATED(spec_ssp)) THEN
+      ALLOCATE(spec_ssp(ntfull, nspec))
+      ALLOCATE(mass_ssp(ntfull))
+      ALLOCATE(lbol_ssp(ntfull))
+      ALLOCATE(ocompsp(ntfull))
+  END IF
+
   !set IMF
+  WRITE(6,*)
   WRITE(6,*)  'enter IMF [0-5; def:0]:'
   WRITE(6,*) ' (0=Salpeter, 1=Chabrier 2003, 2=Kroupa 2001, '//&
        '3=van Dokkum 2008, 4=Dave 2008, 5=tabulated)'
@@ -37,21 +73,6 @@ PROGRAM AUTOSPS
      STOP
   ENDIF
   WRITE(6,'(" ---> Using IMF",1x,I1)') imf_type
-  
-
-  !setup directory and metallicity array
-  CALL GETENV('SPS_HOME',SPS_HOME)
-  IF (LEN_TRIM(SPS_HOME).EQ.0) THEN
-     WRITE(*,*) 'SETUP_SPS ERROR: spsdir environment variable not set!'
-     STOP
-  ENDIF
-  OPEN(90,FILE=TRIM(SPS_HOME)//'/ISOCHRONES/Padova/Padova2007/zlegend.dat',&
-       STATUS='OLD',ACTION='READ')
-  DO z=1,nz
-     READ(90,'(F6.4)') zlegend(z)
-  ENDDO
-  CLOSE(90)
-
 
   !set SFH
   WRITE(6,*)
@@ -65,7 +86,7 @@ PROGRAM AUTOSPS
   ENDIF
   IF (pset%sfh.EQ.0) WRITE(6,'(" ---> Computing an SSP")') 
   IF (pset%sfh.EQ.1) WRITE(6,'(" ---> Computing a CSP")') 
-  IF (pset%sfh.EQ.2) WRITE(6,'(" ---> Computing a tabulated SFH")') 
+  IF (pset%sfh.EQ.2) WRITE(6,'(" ---> Computing a tabulated SFH")')
 
   IF (pset%sfh.EQ.1) THEN
      WRITE(6,*)
@@ -93,10 +114,11 @@ PROGRAM AUTOSPS
   IF (pset%sfh.NE.2) THEN
      !set metallicity
      WRITE(6,*)
-     WRITE(6,*)  'enter metallicity [1-22; def:20]:'
+     ! Dynamic prompt using the loaded nz
+     WRITE(6,'("enter metallicity index [1-",I0,"; def:",I0,"]:")') nz, nz
      READ(5,'(A)')  aux
      IF (len(trim(aux)).EQ.0) THEN
-        pset%zmet = 20
+        pset%zmet = nz  ! Default to the last index (usually safest/solar-ish)
      ELSE
         READ(aux,'(I2)') pset%zmet
      ENDIF
@@ -139,17 +161,22 @@ PROGRAM AUTOSPS
   WRITE(6,'(" ---> Running model.......")')
   
   IF (pset%sfh.EQ.2) THEN
-     CALL SPS_SETUP(-1) 
+     ! We already called SPS_SETUP(-1) at the top, so variables are ready.
      DO z=1,nz
         pset%zmet=z
         CALL SSP_GEN(pset,mass_ssp_zz(:,z),lbol_ssp_zz(:,z),spec_ssp_zz(:,:,z))
      ENDDO
      CALL COMPSP(3,nz,file1,mass_ssp_zz,lbol_ssp_zz,spec_ssp_zz,pset,ocompsp)
   ELSE
-     CALL SPS_SETUP(pset%zmet)
+     ! We already called SPS_SETUP(-1), so speclib is populated.
      CALL SSP_GEN(pset,mass_ssp,lbol_ssp,spec_ssp)
      CALL COMPSP(3,1,file1,mass_ssp,lbol_ssp,spec_ssp,pset,ocompsp)
   ENDIF
 
+  ! Clean up
+  IF (ALLOCATED(spec_ssp)) DEALLOCATE(spec_ssp)
+  IF (ALLOCATED(mass_ssp)) DEALLOCATE(mass_ssp)
+  IF (ALLOCATED(lbol_ssp)) DEALLOCATE(lbol_ssp)
+  IF (ALLOCATED(ocompsp))  DEALLOCATE(ocompsp)
 
-END PROGRAM AUTOSPS 
+END PROGRAM AUTOSPS
diff --git a/src/compsp.f90 b/src/compsp.f90
index ec23b234..cef99a23 100644
--- a/src/compsp.f90
+++ b/src/compsp.f90
@@ -411,12 +411,12 @@ SUBROUTINE COMPSP_SETUP_OUTPUT(write_compsp,pset,outfile,imin,imax)
    ENDIF
 
    !formats
-30 FORMAT('#   SFH: tabulated input, dust=(',F6.2,','F6.2,')')
+30 FORMAT('#   SFH: tabulated input, dust=(',F6.2,',',F6.2,')')
 31 FORMAT('#   log(age) log(mass) Log(lbol) log(SFR) spectra')
 32 FORMAT('#   log(age) log(mass) Log(lbol) log(SFR) mags (see FILTER_LIST)')
 33 FORMAT('#   SFH: Tage=',F6.2,' Gyr, log(tau/Gyr)= ',F6.3,&
         ', const= ',F6.3,', fb= ',F6.3,', tb= ',F6.2,&
-        ' Gyr, sf_start= 'F6.3,', dust=(',F6.2,','F6.2,')')
+        ' Gyr, sf_start= ',F6.3,', dust=(',F6.2,',',F6.2,')')
 34 FORMAT('#   log(age) indices (see allindices.dat)')
 
 
diff --git a/src/csp_gen.f90 b/src/csp_gen.f90
index 16ef70ec..ec200d27 100644
--- a/src/csp_gen.f90
+++ b/src/csp_gen.f90
@@ -59,7 +59,7 @@ subroutine csp_gen(mass_ssp, lbol_ssp, spec_ssp, &
   real(SP), dimension(nspec) :: lw_age, temp_spec !,csp1, csp2
   real(SP), dimension(nemline) :: ncsp1, ncsp2, nlw_age, temp_lin
   real(SP), dimension(ntfull, nzin) :: total_weights
-  real(SP), dimension(ntfull) :: w1=0., w2=0.
+  real(SP), dimension(ntfull) :: w1, w2
   integer :: i, j, k, imin, imax, i_tesc
   type(SFHPARAMS) :: sfhpars
   real(SP) :: m1, m2, frac_linear, mfrac, sfr, fburst
@@ -67,6 +67,8 @@ subroutine csp_gen(mass_ssp, lbol_ssp, spec_ssp, &
   real(SP) :: lbol_age, mass_age ! for mass and lbol weighted ages
 
   ! ------- Setup ----------
+  w1 = 0.
+  w2 = 0.
 
   ! Build a structure containing useful units, numbers, and switches for the
   ! weight calculations. The units of the parameters in the `sfhparams`
diff --git a/src/lesssimple.f90 b/src/lesssimple.f90
index 00affdba..c16b797e 100644
--- a/src/lesssimple.f90
+++ b/src/lesssimple.f90
@@ -10,14 +10,14 @@ PROGRAM LESSSIMPLE
   !    variables not explicitly defined here are defined in sps_vars.f90
   INTEGER :: i
   !define variable for SSP spectrum
-  REAL(SP), DIMENSION(ntfull,nspec)  :: spec_pz
+  REAL(SP), ALLOCATABLE :: spec_pz(:,:)
   !define variables for Mass and Lbol info
-  REAL(SP), DIMENSION(ntfull)    :: mass_pz,lbol_pz
+  REAL(SP), ALLOCATABLE :: mass_pz(:),lbol_pz(:)
   CHARACTER(100) :: file2=''
   !structure containing all necessary parameters
   TYPE(PARAMS) :: pset
   !define structure for CSP spectrum
-  TYPE(COMPSPOUT), DIMENSION(ntfull) :: ocompsp
+  TYPE(COMPSPOUT), ALLOCATABLE :: ocompsp(:)
   REAL(SP) :: zave
 
   !---------------------------------------------------------------!
@@ -32,6 +32,14 @@ PROGRAM LESSSIMPLE
   !here we have to read in all the librarries
   CALL SPS_SETUP(-1)
 
+  ! Allocate memory now that SPS_SETUP has defined ntfull/nspec
+  IF (.NOT. ALLOCATED(spec_pz)) THEN
+      ALLOCATE(spec_pz(ntfull, nspec))
+      ALLOCATE(mass_pz(ntfull))
+      ALLOCATE(lbol_pz(ntfull))
+      ALLOCATE(ocompsp(ntfull))
+  END IF
+
   !compute all SSPs (i.e. at all Zs)
   !nz and the various *ssp_zz arrays are stored 
   !in the common block set up in sps_vars.f90
@@ -56,5 +64,10 @@ PROGRAM LESSSIMPLE
   CALL COMPSP(1,nz,file2,mass_ssp_zz,lbol_ssp_zz,&
           spec_ssp_zz,pset,ocompsp)
 
+  ! Clean up memory before exiting
+  IF (ALLOCATED(spec_pz)) DEALLOCATE(spec_pz)
+  IF (ALLOCATED(mass_pz)) DEALLOCATE(mass_pz)
+  IF (ALLOCATED(lbol_pz)) DEALLOCATE(lbol_pz)
+  IF (ALLOCATED(ocompsp))  DEALLOCATE(ocompsp)
 
 END PROGRAM LESSSIMPLE
diff --git a/src/sfh_weight.f90 b/src/sfh_weight.f90
index 0c5bcc5a..6861b565 100644
--- a/src/sfh_weight.f90
+++ b/src/sfh_weight.f90
@@ -36,7 +36,7 @@ function sfh_weight(sfh, imin, imax)
   integer :: i, istart
   real(SP), dimension(2) :: tlim
   real(SP) :: dt, delta_time, log_tb
-  real(SP), dimension(ntfull) :: tmp_wght=0. !left=0., right=0.
+  real(SP), dimension(ntfull) :: tmp_wght !left=0., right=0.
 
 
   ! Check if this is an SSP.  If so, do simple weights and return.
diff --git a/src/simple.f90 b/src/simple.f90
index 380afd4f..12cfdc95 100644
--- a/src/simple.f90
+++ b/src/simple.f90
@@ -10,14 +10,14 @@ PROGRAM SIMPLE
 
   INTEGER :: i
   !define variable for SSP spectrum
-  REAL(SP), DIMENSION(ntfull,nspec)  :: spec_ssp
+  REAL(SP), ALLOCATABLE :: spec_ssp(:,:)
   !define variables for Mass and Lbol info
-  REAL(SP), DIMENSION(ntfull)    :: mass_ssp,lbol_ssp
+  REAL(SP), ALLOCATABLE :: mass_ssp(:), lbol_ssp(:)
   CHARACTER(100) :: file1='', file2=''
   !structure containing all necessary parameters
   TYPE(PARAMS) :: pset
   !define structure for CSP spectrum
-  TYPE(COMPSPOUT), DIMENSION(ntfull) :: ocompsp
+  TYPE(COMPSPOUT), ALLOCATABLE :: ocompsp(:)
   REAL(SP) :: ssfr6,ssfr7,ssfr8,ave_age
 
   !---------------------------------------------------------------!
@@ -34,6 +34,14 @@ PROGRAM SIMPLE
 
   CALL SPS_SETUP(pset%zmet) !read in the isochrones and spectral libraries
 
+  ! Allocate memory now that SPS_SETUP has defined ntfull/nspec
+  IF (.NOT. ALLOCATED(spec_ssp)) THEN
+      ALLOCATE(spec_ssp(ntfull, nspec))
+      ALLOCATE(mass_ssp(ntfull))
+      ALLOCATE(lbol_ssp(ntfull))
+      ALLOCATE(ocompsp(ntfull))
+  END IF
+
   !define the parameter set.  These are the default values, specified 
   !in sps_vars.f90, but are explicitly included here for transparency
   pset%sfh   = 0     !set SFH to "SSP"
@@ -86,5 +94,11 @@ PROGRAM SIMPLE
   !results are returned in the variables ssfr6,...,ave_age
   CALL SFHSTAT(pset,ocompsp(1),ssfr6,ssfr7,ssfr8,ave_age)
 
+  ! Clean up memory before exiting
+  IF (ALLOCATED(spec_ssp)) DEALLOCATE(spec_ssp)
+  IF (ALLOCATED(mass_ssp)) DEALLOCATE(mass_ssp)
+  IF (ALLOCATED(lbol_ssp)) DEALLOCATE(lbol_ssp)
+  IF (ALLOCATED(ocompsp))  DEALLOCATE(ocompsp)
+
 
 END PROGRAM SIMPLE
diff --git a/src/spec_bin.f90 b/src/spec_bin.f90
index d6359622..abe9c09b 100644
--- a/src/spec_bin.f90
+++ b/src/spec_bin.f90
@@ -5,23 +5,59 @@ PROGRAM SPEC_BIN
 
   USE sps_vars
   IMPLICIT NONE
-  INTEGER  :: z,dumi1,i,j
+  INTEGER  :: z,dumi1,i,j,status
   REAL(SP) :: dumr1,d2,d3
   CHARACTER(6) :: zstype
+  CHARACTER(100) :: arg_spec_type
 
   !----------------------------------------------------------------!
 
   CALL GETENV('SPS_HOME',SPS_HOME)
   
+  ! Read spec_type from command line or default to 'miles'
+  CALL GET_COMMAND_ARGUMENT(1, arg_spec_type, STATUS=status)
+  IF (status /= 0 .OR. LEN_TRIM(arg_spec_type) == 0) THEN
+     spec_type = 'miles'
+     WRITE(*,*) 'No spectral library specified, defaulting to: ', spec_type
+  ELSE
+     spec_type = TRIM(arg_spec_type)
+     WRITE(*,*) 'Using spectral library: ', spec_type
+  END IF
+
+  ! Set dimensions based on spec_type
+  IF (spec_type.EQ.'miles') THEN
+     nzinit=5
+     nspec=5994
+  ELSE IF (spec_type.EQ.'basel') THEN
+     nzinit=6
+     nspec=1963
+  ELSE IF (spec_type(1:5).EQ.'ckc14') THEN
+     nzinit=11
+     nspec=11149
+  ELSE IF (spec_type(1:3).EQ.'c3k') THEN
+     nzinit=11
+     nspec=11149
+  ELSE
+     WRITE(*,*) 'SPEC_BIN ERROR: Unknown spec_type: ', spec_type
+     STOP
+  END IF
+
+  ! Allocate arrays
+  ALLOCATE(zlegendinit(nzinit))
+  ALLOCATE(speclib(nspec,nzinit,ndim_logt,ndim_logg))
+  
   IF (spec_type.EQ.'basel') THEN
      OPEN(90,FILE=TRIM(SPS_HOME)//'/SPECTRA/BaSeL3.1/zlegend.dat',&
           STATUS='OLD',ACTION='READ')
   ELSE IF (spec_type.EQ.'miles') THEN
      OPEN(90,FILE=TRIM(SPS_HOME)//'/SPECTRA/MILES/zlegend.dat',&
           STATUS='OLD',ACTION='READ')
-  ELSE IF (spec_type.EQ.'ckc14') THEN
+  ELSE IF (spec_type.EQ.'ckc14'.OR.spec_type(1:5).EQ.'ckc14') THEN
      OPEN(90,FILE=TRIM(SPS_HOME)//'/SPECTRA/CKC14/zlegend.dat',&
           STATUS='OLD',ACTION='READ')
+  ELSE IF (spec_type(1:3).EQ.'c3k') THEN
+     OPEN(90,FILE=TRIM(SPS_HOME)//'/SPECTRA/C3K/zlegend.dat',&
+          STATUS='OLD',ACTION='READ')
   ENDIF
   DO z=1,nzinit
      READ(90,'(F6.4)') zlegendinit(z)
@@ -58,6 +94,15 @@ PROGRAM SPEC_BIN
              //zstype//'.spectra.bin',FORM='UNFORMATTED',&
              STATUS='REPLACE',access='direct',&
              recl=nspec*ndim_logg*ndim_logt*4)
+             
+     ELSE IF (spec_type(1:3).EQ.'c3k') THEN
+        OPEN(92,FILE=TRIM(SPS_HOME)//'/SPECTRA/C3K/'//spec_type//'_z'&
+             //zstype//'.spectra',FORM='FORMATTED',&
+             STATUS='OLD',ACTION='READ')
+        OPEN(93,FILE=TRIM(SPS_HOME)//'/SPECTRA/C3K/'//spec_type//'_z'&
+             //zstype//'.spectra.bin',FORM='UNFORMATTED',&
+             STATUS='REPLACE',access='direct',&
+             recl=nspec*ndim_logg*ndim_logt*4)
   
     ENDIF
 
diff --git a/src/sps_setup.f90 b/src/sps_setup.f90
index 592d3898..948d9890 100644
--- a/src/sps_setup.f90
+++ b/src/sps_setup.f90
@@ -1,4 +1,4 @@
-SUBROUTINE SPS_SETUP(zin)
+SUBROUTINE SPS_SETUP(zin, isoc_type_in, spec_type_in, dust_type_in)
 
   !read in isochrones and spectral libraries for all metallicities.
   !read in band-pass info and the spectrum for Vega.
@@ -10,9 +10,14 @@ SUBROUTINE SPS_SETUP(zin)
   !in a much faster setup.
 
   USE sps_vars
-  USE sps_utils
+  USE sps_utils, ONLY: locate, linterparr, linterp, tsum, get_tuniv, &
+       get_lumdist, airtovac, sps_takedown
   IMPLICIT NONE
   INTEGER, INTENT(in) :: zin
+  CHARACTER(LEN=*), INTENT(in), OPTIONAL :: isoc_type_in
+  CHARACTER(LEN=*), INTENT(in), OPTIONAL :: spec_type_in
+  CHARACTER(LEN=*), INTENT(in), OPTIONAL :: dust_type_in
+
   INTEGER :: stat=1,n,i,j,m,jj,k,i1,i2,stat2=1
   INTEGER, PARAMETER :: ntlam=1221,nspec_agb=6146,nspec_aringer=9032
   INTEGER, PARAMETER :: nlamwr=1963,nspec_pagb=9281
@@ -22,9 +27,10 @@ SUBROUTINE SPS_SETUP(zin)
   CHARACTER(6) :: zstype
   CHARACTER(5) :: zstype5
   REAL(SP) :: dumr1,d1,d2,logage,x,a,zero=0.0,d,one=1.0,dz,dlam
-  CHARACTER(5), DIMENSION(nz) :: zlegend_str=''
-  CHARACTER(5), DIMENSION(nz_xrb) :: zz_str_xrb=''
-  REAL(SP), DIMENSION(nspec) :: tspec=0.
+  
+  CHARACTER(5), ALLOCATABLE :: zlegend_str(:)
+  CHARACTER(5), ALLOCATABLE :: zz_str_xrb(:)
+  REAL(SP), ALLOCATABLE :: tspec(:)
   REAL(SP), DIMENSION(ntlam) :: tvega_lam=0.,tvega_spec=0.
   REAL(SP), DIMENSION(ntlam) :: tsun_lam=0.,tsun_spec=0.
   REAL(SP), DIMENSION(nlamwr) :: tlamwr=0.,tspecwr=0.
@@ -48,28 +54,278 @@ SUBROUTINE SPS_SETUP(zin)
   REAL(SP), DIMENSION(nspec_aringer,n_agb_car) :: aringer_specinit=0.
   REAL(SP), DIMENSION(nagndust_spec)           :: agndust_lam=0.
   REAL(SP), DIMENSION(nagndust_spec,nagndust)  :: agndust_specinit=0.
-  !REAL(KIND(1.0)), DIMENSION(nspec,nzinit,ndim_logt,ndim_logg) :: speclibinit=0.
-  REAL(KIND(1.0)), allocatable :: speclibinit(:,:,:,:)
-  REAL(SP), DIMENSION(nspec,nzwmb,ndim_wmb_logt,ndim_wmb_logg) :: wmbsi=0.
+  REAL(KIND(1.0)), ALLOCATABLE :: speclibinit(:,:,:,:)
+  REAL(SP), ALLOCATABLE :: wmbsi(:,:,:,:)
   REAL(SP), DIMENSION(nzwmb)     :: zwmb=0.
   REAL(SP), DIMENSION(nspec_wmb) :: wmb_lam=0.
   REAL(SP), DIMENSION(nspec_wmb,ndim_wmb_logt,ndim_wmb_logg) :: wmb_specinit=0.
   REAL(SP), DIMENSION(ntabmax)   :: lsflam=0.,lsfsig=0.
   REAL(SP), DIMENSION(30) :: g03lam=0., g03smc=0.
-  REAL(SP), DIMENSION(nspec_xrb) :: tspec_xrb
+  REAL(SP), ALLOCATABLE :: tspec_xrb(:)
 
   !---------------------------------------------------------------!
   !---------------------------------------------------------------!
 
+  CALL SPS_TAKEDOWN()
+
   IF (verbose.EQ.1) THEN
      WRITE(*,*)
      WRITE(*,*) '    Setting up SPS...'
   ENDIF
 
+  ! Initialize Library Variables
+  IF (PRESENT(isoc_type_in)) THEN
+     isoc_type = isoc_type_in
+  ELSE
+     isoc_type = 'mist'
+  END IF
+  
+  IF (PRESENT(spec_type_in)) THEN
+     spec_type = spec_type_in
+  ELSE
+     spec_type = 'miles'
+  END IF
+
+  IF (TRIM(isoc_type) == 'bpss' .AND. TRIM(spec_type) /= 'bpass') THEN
+     WRITE(*,*) 'Notice: BPASS isochrones selected; forcing spec_type="bpass"'
+     spec_type = 'bpass'
+  END IF
+
+  ! Set isochrone dimensions
+  SELECT CASE (isoc_type)
+  CASE ('mist')
+     zsol = 0.0142
+     nt=107
+     nz=12
+  CASE ('pdva')
+     zsol = 0.019
+     nt=94
+     nz=22
+  CASE ('prsc')
+     zsol = 0.01524
+     nt=93
+     nz=15
+  CASE ('bsti')
+     zsol = 0.020
+     nt=94
+     nz=10
+  CASE ('gnva')
+     zsol = 0.020
+     nt=51
+     nz=5
+  CASE ('bpss')
+     zsol = 0.020
+     nt=43
+     nz=12
+  CASE DEFAULT
+      WRITE(*,*) 'SPS_SETUP ERROR: Unknown isoc_type: ', isoc_type
+      STOP
+  END SELECT
+  
+  ! Set spectral library dimensions
+  SELECT CASE (spec_type)
+  CASE ('miles')
+     zsol_spec = 0.019
+     nzinit=5
+     nspec=5994
+  CASE ('basel')
+     zsol_spec = 0.020
+     nzinit=6
+     nspec=1963
+  CASE ('bpass')
+     zsol_spec = 0.020
+     nzinit=1
+     nspec=15000
+  CASE DEFAULT
+     ! Check for C3K
+     IF (spec_type(1:3).EQ.'c3k') THEN
+        zsol_spec = 0.0134
+        nzinit=11
+        nspec=11149
+     ELSE
+        WRITE(*,*) 'SPS_SETUP ERROR: Unknown spec_type: ', spec_type
+        STOP
+     END IF
+  END SELECT
+
+  ! Set other dimensions
+  nbands = 159
+  nindx = 30
+  ntfull = time_res_incr * nt
+  
+  ! Set XRB dimensions
+  nspec_xrb=15000
+  nt_xrb=10
+  nz_xrb=11
+
+  ! Set Dust Emission Model
+  IF (PRESENT(dust_type_in)) THEN
+     IF (dust_type_in == 'themis' .OR. dust_type_in == 'THEMIS') THEN
+        str_dustem = 'THEMIS'
+        ndim_dustem = 576
+        numin_dustem = 37
+        nqpah_dustem = 11
+     ELSE
+        ! Default to DL07
+        str_dustem = 'DL07'
+        ndim_dustem = 1001
+        numin_dustem = 22
+        nqpah_dustem = 7
+     END IF
+  ELSE
+     str_dustem = 'DL07'
+     ndim_dustem = 1001
+     numin_dustem = 22
+     nqpah_dustem = 7
+  END IF
+
   IF (zin.GT.nz) THEN
      WRITE(*,*) 'SPS_SETUP ERROR: zin GT nz', zin,nz
      STOP
   ENDIF
+  
+
+  ! Allocate arrays
+  ALLOCATE(indexdefined(7,nindx))
+  ALLOCATE(wgdust(nspec,18,6,2))
+  ALLOCATE(g03smcextn(nspec))
+  ALLOCATE(bands(nspec,nbands))
+  ALLOCATE(magsun(nbands),magvega(nbands),filter_leff(nbands))
+  ALLOCATE(vega_spec(nspec),sun_spec(nspec))
+  ALLOCATE(spec_lambda(nspec),spec_nu(nspec))
+  ALLOCATE(spec_res(nspec))
+  ALLOCATE(speclib(nspec,nz,ndim_logt,ndim_logg))
+  ALLOCATE(wmb_spec(nspec,nz,ndim_wmb_logt,ndim_wmb_logg))
+  ALLOCATE(agb_spec_o(nspec,n_agb_o))
+  ALLOCATE(agb_logt_o(nz,n_agb_o))
+  ALLOCATE(agb_spec_c(nspec,n_agb_c))
+  ALLOCATE(agb_logt_c(n_agb_c))
+  ALLOCATE(agb_spec_car(nspec,n_agb_car))
+  ALLOCATE(pagb_spec(nspec,ndim_pagb,2))
+  ALLOCATE(wrn_spec(nspec,ndim_wr,nz),wrc_spec(nspec,ndim_wr,nz))
+  
+  ! Dust models
+  ALLOCATE(qpaharr(nqpah_dustem))
+  ALLOCATE(uminarr(numin_dustem))
+  ALLOCATE(lambda_dustem(ndim_dustem))
+  ALLOCATE(dustem_dustem(ndim_dustem,numin_dustem*2))
+  ALLOCATE(dustem2_dustem(nspec,nqpah_dustem,numin_dustem*2))
+  
+  ALLOCATE(flux_dagb(nspec,2,nteff_dagb,ntau_dagb))
+  ALLOCATE(nebem_cont(nspec,nebnz,nebnage,nebnip),xnebem_cont(nspec,nebnz,nebnage,nebnip))
+  ALLOCATE(neb_res_min(nspec))
+  ALLOCATE(gaussnebarr(nspec,nemline))
+  ALLOCATE(agndust_spec(nspec,nagndust))
+  ALLOCATE(mact_isoc(nz,nt,nm),logl_isoc(nz,nt,nm),logt_isoc(nz,nt,nm),logg_isoc(nz,nt,nm))
+  ALLOCATE(ffco_isoc(nz,nt,nm),phase_isoc(nz,nt,nm),mini_isoc(nz,nt,nm),lmdot_isoc(nz,nt,nm))
+  ALLOCATE(nmass_isoc(nz,nt))
+  ALLOCATE(timestep_isoc(nz,nt))
+  ALLOCATE(zlegend(nz))
+  ALLOCATE(zlegendinit(nzinit))
+  ALLOCATE(spec_ssp_zz(nspec,ntfull,nz))
+  ALLOCATE(mass_ssp_zz(ntfull,nz),lbol_ssp_zz(ntfull,nz))
+  ALLOCATE(time_full(ntfull))
+  ALLOCATE(weight_ssp(ntfull,nz))
+  ALLOCATE(spec_young(nspec),spec_old(nspec))
+  ALLOCATE(bpass_spec_ssp(nspec,nt,nz))
+  ALLOCATE(bpass_mass_ssp(nt,nz))
+  ALLOCATE(lam_xrb(nspec_xrb))
+  ALLOCATE(spec_xrb(nspec,nt_xrb,nz_xrb))
+  ALLOCATE(ages_xrb(nt_xrb))
+  ALLOCATE(zmet_xrb(nz_xrb))
+  
+  ALLOCATE(lsfinfo%lsf(nspec))
+  
+  ! Allocate local arrays
+  ALLOCATE(tspec(nspec))
+  ALLOCATE(wmbsi(nspec,nzwmb,ndim_wmb_logt,ndim_wmb_logg))
+  ALLOCATE(tspec_xrb(nspec_xrb))
+  ALLOCATE(zlegend_str(nz))
+  ALLOCATE(zz_str_xrb(nz_xrb))
+
+  ! Initialize new arrays to 0.0 or default values
+  bands = 0.0
+  magsun = 0.0
+  magvega = 0.0
+  filter_leff = 0.0
+  vega_spec = 0.0
+  sun_spec = 0.0
+  spec_lambda = 0.0
+  spec_nu = 0.0
+  spec_res = 0.0
+  speclib = 0.0
+  wmb_spec = 0.0
+  agb_spec_o = 0.0
+  agb_logt_o = 0.0
+  agb_spec_c = 0.0
+  agb_logt_c = 0.0
+  agb_spec_car = 0.0
+  pagb_spec = 0.0
+  wrn_spec = 0.0
+  wrc_spec = 0.0
+  
+  ! Dust initialization
+  dustem2_dustem = 0.0
+  
+  IF (TRIM(str_dustem) == 'THEMIS') THEN
+     qpaharr = (/0.02,0.06,0.10,0.14,0.17,0.20,0.24,0.28,0.32,0.36,0.40/)/2.2*100
+     uminarr = (/0.1,0.12,0.15,0.17,0.2,0.25,0.3,0.35,0.4,0.5,0.6,0.7,0.8,1.0,&
+       1.2,1.5,1.7, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0,&
+       12.0, 15.0, 17.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 80.0/)
+  ELSE
+     ! DL07
+     qpaharr = (/0.47,1.12,1.77,2.50,3.19,3.90,4.58/)
+     uminarr = (/0.1,0.15,0.2,0.3,0.4,0.5,0.7,0.8,1.0,1.2,1.5,2.0,&
+       2.5,3.0,4.0,5.0,7.0,8.0,12.0,15.0,20.0,25.0/)
+  END IF
+  
+  lambda_dustem = 0.0
+  dustem_dustem = 0.0
+  
+  flux_dagb = 0.0
+  nebem_cont = 0.0
+  xnebem_cont = 0.0
+  neb_res_min = 0.0
+  gaussnebarr = 0.0
+  agndust_spec = 0.0
+  mact_isoc = 0.0
+  logl_isoc = 0.0
+  logt_isoc = 0.0
+  logg_isoc = 0.0
+  ffco_isoc = 0.0
+  phase_isoc = 0.0
+  mini_isoc = 0.0
+  lmdot_isoc = 0.0
+  nmass_isoc = 0
+  timestep_isoc = 0.0
+  zlegend = -99.0
+  zlegendinit = -99.0
+  spec_ssp_zz = 0.0
+  mass_ssp_zz = 0.0
+  lbol_ssp_zz = 0.0
+  time_full = 0.0
+  weight_ssp = 0.0
+  spec_young = 0.0
+  spec_old = 0.0
+  bpass_spec_ssp = 0.0
+  bpass_mass_ssp = 0.0
+  lam_xrb = 0.0
+  spec_xrb = 0.0
+  ages_xrb = 0.0
+  zmet_xrb = 0.0
+  lsfinfo%lsf = 0.0
+  
+  tspec = 0.0
+  wmbsi = 0.0
+  tspec_xrb = 0.0
+  zlegend_str = ''
+  zz_str_xrb = ''
+  indexdefined = 0.0
+  mwdindex = 0
+  wgdust = 0.0
+  g03smcextn = 0.0
+  sfh_tab = 0.0
+  ntabsfh = 0
 
   !clean out all the common block arrays
   mini_isoc     = 0.
@@ -786,7 +1042,7 @@ SUBROUTINE SPS_SETUP(zin)
      !now interpolate the dust spectra onto the master wavelength array
      DO j=1,numin_dustem*2
         !the dust models only extend to 1um
-        jj = locate(spec_lambda/1E4,one)
+        jj = MAX(locate(spec_lambda/1E4,one), 1)
         dustem2_dustem(jj:,k,j) = linterparr(lambda_dustem,&
              dustem_dustem(:,j),spec_lambda(jj:))
      ENDDO
@@ -821,7 +1077,7 @@ SUBROUTINE SPS_SETUP(zin)
            STOP
         ENDIF
         !interpolate the dust spectra onto the master wavelength array
-        jj = locate(spec_lambda,lambda_dagb(1))
+        jj = MAX(locate(spec_lambda,lambda_dagb(1)), 1)
         flux_dagb(jj:,1,i,j) = linterparr(lambda_dagb(1:nlam),&
              fluxin_dagb(1:nlam),spec_lambda(jj:))
      ENDDO
@@ -850,7 +1106,7 @@ SUBROUTINE SPS_SETUP(zin)
            STOP
         ENDIF
         !interpolate the dust spectra onto the master wavelength array
-        jj = locate(spec_lambda,lambda_dagb(1))
+        jj = MAX(locate(spec_lambda,lambda_dagb(1)), 1)
         flux_dagb(jj:,2,i,j) = linterparr(lambda_dagb(1:nlam),&
              fluxin_dagb(1:nlam),spec_lambda(jj:))
      ENDDO
@@ -882,8 +1138,8 @@ SUBROUTINE SPS_SETUP(zin)
      READ(99,*) agndust_lam(i),agndust_specinit(i,:)
   ENDDO
 
-  i1 = locate(spec_lambda,agndust_lam(1))
-  i2 = locate(spec_lambda,agndust_lam(nagndust_spec))
+  i1 = MAX(locate(spec_lambda,agndust_lam(1)), 1)
+  i2 = MAX(locate(spec_lambda,agndust_lam(nagndust_spec)), 1)
   DO i=1,nagndust
      agndust_spec(i1:i2,i) = 10**linterparr(LOG10(agndust_lam),&
           LOG10(agndust_specinit(:,i)+tiny30),LOG10(spec_lambda(i1:i2)))-tiny30
@@ -1115,7 +1371,7 @@ SUBROUTINE SPS_SETUP(zin)
 
   !interpolate the Vega spectrum onto the wavelength grid
   !and convert to fnu
-  jj = locate(vega_spec,tvega_lam(ntlam))
+  jj = locate(spec_lambda,tvega_lam(ntlam))
   vega_spec(:jj) = 10**linterparr(LOG10(tvega_lam),&
           LOG10(tvega_spec+tiny_number),LOG10(spec_lambda(:jj)))
   vega_spec = vega_spec*spec_lambda**2
@@ -1136,7 +1392,7 @@ SUBROUTINE SPS_SETUP(zin)
   CLOSE(98)
 
   !interpolate the Solar spectrum onto the wavelength grid
-  jj = locate(sun_spec,tsun_lam(ntlam))
+  jj = locate(spec_lambda,tsun_lam(ntlam))
   sun_spec(:jj) = 10**linterparr(LOG10(tsun_lam),&
           LOG10(tsun_spec+tiny_number),LOG10(spec_lambda(:jj)))
   sun_spec(jj+1:) = tiny_number
diff --git a/src/sps_utils.f90 b/src/sps_utils.f90
index 0f646bb9..6ef8bccc 100644
--- a/src/sps_utils.f90
+++ b/src/sps_utils.f90
@@ -1,5 +1,15 @@
 MODULE SPS_UTILS
 
+  INTERFACE
+     SUBROUTINE SPS_SETUP(zin, isoc_type_in, spec_type_in, dust_type_in)
+       USE sps_vars
+       INTEGER, INTENT(in) :: zin
+       CHARACTER(LEN=*), INTENT(in), OPTIONAL :: isoc_type_in
+       CHARACTER(LEN=*), INTENT(in), OPTIONAL :: spec_type_in
+       CHARACTER(LEN=*), INTENT(in), OPTIONAL :: dust_type_in
+     END SUBROUTINE SPS_SETUP
+  END INTERFACE
+
   INTERFACE
      SUBROUTINE ADD_AGB_DUST(weight,tspec,mact,logt,logl,logg,&
           zz,tco,lmdot)
@@ -380,4 +390,105 @@ SUBROUTINE ZTINTERP(zpos,spec,lbol,mass,tpos,zpow)
      END SUBROUTINE ZTINTERP
   END INTERFACE
 
+CONTAINS
+
+  SUBROUTINE SPS_TAKEDOWN()
+    USE sps_vars
+    IMPLICIT NONE
+
+    ! Deallocate all arrays
+    IF (ALLOCATED(isoc_type)) DEALLOCATE(isoc_type)
+    IF (ALLOCATED(spec_type)) DEALLOCATE(spec_type)
+    IF (ALLOCATED(indexdefined)) DEALLOCATE(indexdefined)
+    IF (ALLOCATED(wgdust)) DEALLOCATE(wgdust)
+    IF (ALLOCATED(g03smcextn)) DEALLOCATE(g03smcextn)
+    IF (ALLOCATED(bands)) DEALLOCATE(bands)
+    IF (ALLOCATED(magsun)) DEALLOCATE(magsun)
+    IF (ALLOCATED(magvega)) DEALLOCATE(magvega)
+    IF (ALLOCATED(filter_leff)) DEALLOCATE(filter_leff)
+    IF (ALLOCATED(vega_spec)) DEALLOCATE(vega_spec)
+    IF (ALLOCATED(sun_spec)) DEALLOCATE(sun_spec)
+    IF (ALLOCATED(spec_lambda)) DEALLOCATE(spec_lambda)
+    IF (ALLOCATED(spec_nu)) DEALLOCATE(spec_nu)
+    IF (ALLOCATED(spec_res)) DEALLOCATE(spec_res)
+    IF (ALLOCATED(speclib)) DEALLOCATE(speclib)
+    IF (ALLOCATED(wmb_spec)) DEALLOCATE(wmb_spec)
+    IF (ALLOCATED(agb_spec_o)) DEALLOCATE(agb_spec_o)
+    IF (ALLOCATED(agb_logt_o)) DEALLOCATE(agb_logt_o)
+    IF (ALLOCATED(agb_spec_c)) DEALLOCATE(agb_spec_c)
+    IF (ALLOCATED(agb_logt_c)) DEALLOCATE(agb_logt_c)
+    IF (ALLOCATED(agb_spec_car)) DEALLOCATE(agb_spec_car)
+    IF (ALLOCATED(pagb_spec)) DEALLOCATE(pagb_spec)
+    IF (ALLOCATED(wrn_spec)) DEALLOCATE(wrn_spec)
+    IF (ALLOCATED(wrc_spec)) DEALLOCATE(wrc_spec)
+    IF (ALLOCATED(qpaharr)) DEALLOCATE(qpaharr)
+    IF (ALLOCATED(uminarr)) DEALLOCATE(uminarr)
+    IF (ALLOCATED(lambda_dustem)) DEALLOCATE(lambda_dustem)
+    IF (ALLOCATED(dustem_dustem)) DEALLOCATE(dustem_dustem)
+    IF (ALLOCATED(dustem2_dustem)) DEALLOCATE(dustem2_dustem)
+    IF (ALLOCATED(flux_dagb)) DEALLOCATE(flux_dagb)
+    IF (ALLOCATED(nebem_cont)) DEALLOCATE(nebem_cont)
+    IF (ALLOCATED(xnebem_cont)) DEALLOCATE(xnebem_cont)
+    IF (ALLOCATED(neb_res_min)) DEALLOCATE(neb_res_min)
+    IF (ALLOCATED(gaussnebarr)) DEALLOCATE(gaussnebarr)
+    IF (ALLOCATED(agndust_spec)) DEALLOCATE(agndust_spec)
+    IF (ALLOCATED(mact_isoc)) DEALLOCATE(mact_isoc)
+    IF (ALLOCATED(logl_isoc)) DEALLOCATE(logl_isoc)
+    IF (ALLOCATED(logt_isoc)) DEALLOCATE(logt_isoc)
+    IF (ALLOCATED(logg_isoc)) DEALLOCATE(logg_isoc)
+    IF (ALLOCATED(ffco_isoc)) DEALLOCATE(ffco_isoc)
+    IF (ALLOCATED(phase_isoc)) DEALLOCATE(phase_isoc)
+    IF (ALLOCATED(mini_isoc)) DEALLOCATE(mini_isoc)
+    IF (ALLOCATED(lmdot_isoc)) DEALLOCATE(lmdot_isoc)
+    IF (ALLOCATED(nmass_isoc)) DEALLOCATE(nmass_isoc)
+    IF (ALLOCATED(timestep_isoc)) DEALLOCATE(timestep_isoc)
+    IF (ALLOCATED(zlegend)) DEALLOCATE(zlegend)
+    IF (ALLOCATED(zlegendinit)) DEALLOCATE(zlegendinit)
+    IF (ALLOCATED(spec_ssp_zz)) DEALLOCATE(spec_ssp_zz)
+    IF (ALLOCATED(mass_ssp_zz)) DEALLOCATE(mass_ssp_zz)
+    IF (ALLOCATED(lbol_ssp_zz)) DEALLOCATE(lbol_ssp_zz)
+    IF (ALLOCATED(time_full)) DEALLOCATE(time_full)
+    IF (ALLOCATED(weight_ssp)) DEALLOCATE(weight_ssp)
+    IF (ALLOCATED(spec_young)) DEALLOCATE(spec_young)
+    IF (ALLOCATED(spec_old)) DEALLOCATE(spec_old)
+    IF (ALLOCATED(bpass_spec_ssp)) DEALLOCATE(bpass_spec_ssp)
+    IF (ALLOCATED(bpass_mass_ssp)) DEALLOCATE(bpass_mass_ssp)
+    IF (ALLOCATED(lam_xrb)) DEALLOCATE(lam_xrb)
+    IF (ALLOCATED(spec_xrb)) DEALLOCATE(spec_xrb)
+    IF (ALLOCATED(ages_xrb)) DEALLOCATE(ages_xrb)
+    IF (ALLOCATED(zmet_xrb)) DEALLOCATE(zmet_xrb)
+
+    IF (ALLOCATED(lsfinfo%lsf)) DEALLOCATE(lsfinfo%lsf)
+
+    IF (ALLOCATED(powell_data%mags)) DEALLOCATE(powell_data%mags)
+    IF (ALLOCATED(powell_data%magerr)) DEALLOCATE(powell_data%magerr)
+    IF (ALLOCATED(powell_data%spec)) DEALLOCATE(powell_data%spec)
+    IF (ALLOCATED(powell_data%specerr)) DEALLOCATE(powell_data%specerr)
+
+    IF (ALLOCATED(sedfit_data%mags)) DEALLOCATE(sedfit_data%mags)
+    IF (ALLOCATED(sedfit_data%magerr)) DEALLOCATE(sedfit_data%magerr)
+    IF (ALLOCATED(sedfit_data%spec)) DEALLOCATE(sedfit_data%spec)
+    IF (ALLOCATED(sedfit_data%specerr)) DEALLOCATE(sedfit_data%specerr)
+
+    ! Reset dimensions
+    nspec = 0
+    nt = 0
+    nz = 0
+    nbands = 0
+    nindx = 0
+    ntfull = 0
+    nzinit = 0
+    nspec_xrb = 0
+    nt_xrb = 0
+    nz_xrb = 0
+    ndim_dustem = 0
+    numin_dustem = 0
+    nqpah_dustem = 0
+    n_user_imf = 0
+
+    ! Reset flag
+    check_sps_setup = 0
+
+  END SUBROUTINE SPS_TAKEDOWN
+
 END MODULE SPS_UTILS
diff --git a/src/sps_vars.f90 b/src/sps_vars.f90
index 8f8bbc08..9f969a3a 100644
--- a/src/sps_vars.f90
+++ b/src/sps_vars.f90
@@ -5,55 +5,6 @@ MODULE SPS_VARS
   IMPLICIT NONE
   SAVE
 
-!-------set the spectral library------!
-#ifndef MILES
-#define MILES 1
-#endif
-
-#ifndef BASEL
-#define BASEL 0
-#endif
-
-#ifndef C3K
-#define C3K 0
-#endif
-
-!------set the isochrone library------!
-#ifndef MIST
-#define MIST 1
-#endif
-
-#ifndef PADOVA
-#define PADOVA 0
-#endif
-
-#ifndef PARSEC
-#define PARSEC 0
-#endif
-
-#ifndef BASTI
-#define BASTI 0
-#endif
-
-#ifndef GENEVA
-#define GENEVA 0
-#endif
-
-!note that in the case of BPASS the SSPs are already pre-computed
-!so the spectral library, IMF, etc. is fixed in this case.
-#ifndef BPASS
-#define BPASS 0
-#endif
-
-!------set the dust emission model------!
-#ifndef DL07
-#define DL07 1
-#endif
-
-#ifndef THEMIS
-#define THEMIS 0
-#endif
-
   !--------------------------------------------------------------!
   !--------------------------------------------------------------!
 
@@ -211,66 +162,21 @@ MODULE SPS_VARS
   INTEGER :: nebemlineinspec=1
 
   !------------Pre-compiler defintions------------!
+  ! Replaced by runtime variables
 
   !flag indicating type of isochrones to use
   !and number of metallicities in the set
-#if (BASTI)
-  REAL(SP), PARAMETER :: zsol = 0.020
-  CHARACTER(4), PARAMETER :: isoc_type = 'bsti'
-  INTEGER, PARAMETER :: nt=94
-  INTEGER, PARAMETER :: nz=10
-#elif (GENEVA)
-  REAL(SP), PARAMETER :: zsol = 0.020
-  CHARACTER(4), PARAMETER :: isoc_type = 'gnva'
-  INTEGER, PARAMETER :: nt=51
-  INTEGER, PARAMETER :: nz=5
-#elif (MIST)
-  REAL(SP), PARAMETER :: zsol = 0.0142
-  CHARACTER(4), PARAMETER :: isoc_type = 'mist'
-  INTEGER, PARAMETER :: nt=107
-  INTEGER, PARAMETER :: nz=12
-#elif (PARSEC)
-  REAL(SP), PARAMETER :: zsol = 0.01524
-  CHARACTER(4), PARAMETER :: isoc_type = 'prsc'
-  INTEGER, PARAMETER :: nt=93
-  INTEGER, PARAMETER :: nz=15
-#elif (PADOVA)
-  REAL(SP), PARAMETER :: zsol = 0.019
-  CHARACTER(4), PARAMETER :: isoc_type = 'pdva'
-  INTEGER, PARAMETER :: nt=94
-  INTEGER, PARAMETER :: nz=22
-#elif (BPASS)
-  REAL(SP), PARAMETER :: zsol = 0.020
-  CHARACTER(4), PARAMETER :: isoc_type = 'bpss'
-  INTEGER, PARAMETER :: nt=43
-  INTEGER, PARAMETER :: nz=12
-#endif
+  REAL(SP) :: zsol = 0.0
+  CHARACTER(LEN=:), ALLOCATABLE :: isoc_type
+  INTEGER :: nt=0
+  INTEGER :: nz=0
 
   !flag indicating type of spectral library to use
   !and number of elements per stellar spectrum
-#if (BPASS)
-  REAL(SP), PARAMETER :: zsol_spec = 0.020
-  CHARACTER(5), PARAMETER :: spec_type = 'bpass'
-  INTEGER, PARAMETER :: nzinit=1
-  INTEGER, PARAMETER :: nspec=15000
-#else
-#if (MILES)
-  REAL(SP), PARAMETER :: zsol_spec = 0.019
-  CHARACTER(5), PARAMETER :: spec_type = 'miles'
-  INTEGER, PARAMETER :: nzinit=5
-  INTEGER, PARAMETER :: nspec=5994
-#elif (C3K)
-  REAL(SP), PARAMETER :: zsol_spec = 0.0134
-  CHARACTER(11), PARAMETER :: spec_type = 'c3k_afe+0.0'
-  INTEGER, PARAMETER :: nzinit=11
-  INTEGER, PARAMETER :: nspec=11149
-#elif (BASEL)
-  REAL(SP), PARAMETER :: zsol_spec = 0.020
-  CHARACTER(5), PARAMETER :: spec_type = 'basel'
-  INTEGER, PARAMETER :: nzinit=6
-  INTEGER, PARAMETER :: nspec=1963
-#endif
-#endif
+  REAL(SP) :: zsol_spec = 0.0
+  CHARACTER(LEN=:), ALLOCATABLE :: spec_type
+  INTEGER :: nzinit=0
+  INTEGER :: nspec=0
 
   !flag indicating the type of normalization used in the BaSeL library
   !pdva = normalized to Padova isochrones
@@ -282,9 +188,9 @@ MODULE SPS_VARS
 
   !You must change the number of bands here if
   !filters are added to allfilters.dat
-  INTEGER, PARAMETER :: nbands=159
+  INTEGER :: nbands=0
   !number of indices defined in allindices.dat
-  INTEGER, PARAMETER :: nindx=30
+  INTEGER :: nindx=0
 
   !The following parameters should never be changed
   !unless you are changing the libraries
@@ -316,9 +222,9 @@ MODULE SPS_VARS
   !number of spectral points in the input library
   INTEGER, PARAMETER :: nagndust_spec=125
 
-  INTEGER, PARAMETER :: nspec_xrb=15000
-  INTEGER, PARAMETER :: nt_xrb=10
-  INTEGER, PARAMETER :: nz_xrb=11
+  INTEGER :: nspec_xrb=0
+  INTEGER :: nt_xrb=0
+  INTEGER :: nz_xrb=0
 
   !------------IMF-related Constants--------------!
 
@@ -407,20 +313,20 @@ MODULE SPS_VARS
   INTEGER :: whlam5000,whlylim
 
   !this specifies the size of the full time grid
-  INTEGER, PARAMETER :: ntfull = time_res_incr*nt
+  INTEGER :: ntfull = 0
 
   !array of index definitions
-  REAL(SP), DIMENSION(7,nindx) :: indexdefined=0.
+  REAL(SP), ALLOCATABLE :: indexdefined(:,:)
 
   !array holding MW extinction curve indices
   INTEGER, DIMENSION(6) :: mwdindex=0
 
   !array holding Witt & Gordon dust models
   !wgdust(lam,tau,model,homo/clump)
-  REAL(SP), DIMENSION(nspec,18,6,2) :: wgdust=0.
+  REAL(SP), ALLOCATABLE :: wgdust(:,:,:,:)
 
   !array holding the Gordon et al. (2003) SMC extinction
-  REAL(SP), DIMENSION(nspec) :: g03smcextn=0.
+  REAL(SP), ALLOCATABLE :: g03smcextn(:)
 
   !Index for P(Z) distribution.  1=closed box;
   !P(Z) = z^zpow*exp(-z/pmetals)  (see pz_convol.f90)
@@ -435,124 +341,109 @@ MODULE SPS_VARS
   INTEGER :: ntabsfh=0
 
   !array of bandpass filters
-  REAL(SP), DIMENSION(nspec,nbands) :: bands
+  REAL(SP), ALLOCATABLE :: bands(:,:)
   !magnitude of the Sun in all filters
-  REAL(SP), DIMENSION(nbands) :: magsun,magvega,filter_leff
+  REAL(SP), ALLOCATABLE :: magsun(:),magvega(:),filter_leff(:)
   !Vega-like star spectrum for Vega magnitude zero-point
   !spectrum of Sun, for absolute mags of Sun
-  REAL(SP), DIMENSION(nspec)  :: vega_spec=0.,sun_spec=0.
+  REAL(SP), ALLOCATABLE  :: vega_spec(:),sun_spec(:)
   !common wavelength and frequench arrays
-  REAL(SP), DIMENSION(nspec)  :: spec_lambda=0.,spec_nu=0.0
+  REAL(SP), ALLOCATABLE  :: spec_lambda(:),spec_nu(:)
   !common wavelength and frequency arrays for dummy resolution files
-  REAL(SP), DIMENSION(nspec)  :: spec_res=0.
+  REAL(SP), ALLOCATABLE  :: spec_res(:)
 
   !arrays for stellar spectral information in HR diagram
   REAL(SP), DIMENSION(ndim_logt) :: speclib_logt=0.
   REAL(SP), DIMENSION(ndim_logg) :: speclib_logg=0.
-  REAL(KIND(1.0)), DIMENSION(nspec,nz,ndim_logt,ndim_logg) :: speclib=0.
+  REAL(KIND(1.0)), ALLOCATABLE :: speclib(:,:,:,:)
 
   !arrays for the WMBasic grid
   REAL(SP), DIMENSION(ndim_wmb_logt) :: wmb_logt=0.
   REAL(SP), DIMENSION(ndim_wmb_logg) :: wmb_logg=0.
-  REAL(KIND(1.0)), DIMENSION(nspec,nz,ndim_wmb_logt,ndim_wmb_logg) :: wmb_spec=0.
+  REAL(KIND(1.0)), ALLOCATABLE :: wmb_spec(:,:,:,:)
 
   !AGB library (Lancon & Mouhcine 2002)
-  REAL(SP), DIMENSION(nspec,n_agb_o) :: agb_spec_o=0.
-  REAL(SP), DIMENSION(nz,n_agb_o)    :: agb_logt_o=0.
-  REAL(SP), DIMENSION(nspec,n_agb_c) :: agb_spec_c=0.
-  REAL(SP), DIMENSION(n_agb_c)       :: agb_logt_c=0.
+  REAL(SP), ALLOCATABLE :: agb_spec_o(:,:)
+  REAL(SP), ALLOCATABLE :: agb_logt_o(:,:)
+  REAL(SP), ALLOCATABLE :: agb_spec_c(:,:)
+  REAL(SP), ALLOCATABLE :: agb_logt_c(:)
   !C-rich library (Aringer et al. 2009)
   REAL(SP), DIMENSION(n_agb_car)       :: agb_logt_car=0.
-  REAL(SP), DIMENSION(nspec,n_agb_car) :: agb_spec_car=0.
+  REAL(SP), ALLOCATABLE :: agb_spec_car(:,:)
 
   !post-AGB library (Rauch 2003)
-  REAL(SP), DIMENSION(nspec,ndim_pagb,2) :: pagb_spec=0.
+  REAL(SP), ALLOCATABLE :: pagb_spec(:,:,:)
   REAL(SP), DIMENSION(ndim_pagb)         :: pagb_logt=0.
 
   !WR library (Smith et al. 2002)
-  REAL(SP), DIMENSION(nspec,ndim_wr,nz) :: wrn_spec=0.,wrc_spec=0.
+  REAL(SP), ALLOCATABLE :: wrn_spec(:,:,:),wrc_spec(:,:,:)
   REAL(SP), DIMENSION(ndim_wr)          :: wrn_logt=0.,wrc_logt=0.
 
-#if (DL07)
-  !dust emission model (Draine & Li 2007)
-  INTEGER, PARAMETER :: ndim_dustem=1001
-  INTEGER, PARAMETER :: numin_dustem=22, nqpah_dustem=7
-  CHARACTER(6), PARAMETER :: str_dustem='DL07'
-  REAL(SP), DIMENSION(nqpah_dustem), PARAMETER :: &
-       qpaharr = (/0.47,1.12,1.77,2.50,3.19,3.90,4.58/)
-  REAL(SP), DIMENSION(numin_dustem) :: uminarr = &
-       (/0.1,0.15,0.2,0.3,0.4,0.5,0.7,0.8,1.0,1.2,1.5,2.0,&
-       2.5,3.0,4.0,5.0,7.0,8.0,12.0,15.0,20.0,25.0/)
-#elif (THEMIS)
-  !dust emission model (THEMIS; Jones et al. 2013, 2017)
-  INTEGER, PARAMETER :: ndim_dustem=576
-  INTEGER, PARAMETER :: numin_dustem=37, nqpah_dustem=11
-  CHARACTER(6), PARAMETER :: str_dustem='THEMIS'
-  REAL(SP), DIMENSION(nqpah_dustem), PARAMETER :: &
-       qpaharr = (/0.02,0.06,0.10,0.14,0.17,0.20,0.24,0.28,0.32,0.36,0.40/)/2.2*100
-  REAL(SP), DIMENSION(numin_dustem) :: uminarr = &
-       (/0.1,0.12,0.15,0.17,0.2,0.25,0.3,0.35,0.4,0.5,0.6,0.7,0.8,1.0,&
-       1.2,1.5,1.7, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0,&
-       12.0, 15.0, 17.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 80.0/)
-#endif
-
-  REAL(SP), DIMENSION(ndim_dustem)                :: lambda_dustem=0.
-  REAL(SP), DIMENSION(ndim_dustem,numin_dustem*2) :: dustem_dustem=0.
-  REAL(SP), DIMENSION(nspec,nqpah_dustem,numin_dustem*2) :: dustem2_dustem=0.
+  !dust emission model (Draine & Li 2007 or THEMIS)
+  INTEGER :: ndim_dustem=0
+  INTEGER :: numin_dustem=0, nqpah_dustem=0
+  CHARACTER(6) :: str_dustem='DL07'
+  
+  REAL(SP), ALLOCATABLE :: qpaharr(:)
+  REAL(SP), ALLOCATABLE :: uminarr(:)
+
+  REAL(SP), ALLOCATABLE :: lambda_dustem(:)
+  REAL(SP), ALLOCATABLE :: dustem_dustem(:,:)
+  REAL(SP), ALLOCATABLE :: dustem2_dustem(:,:,:)
 
   !circumstellar AGB dust model (Villaume et al. 2015)
-  REAL(SP), DIMENSION(nspec,2,nteff_dagb,ntau_dagb) :: flux_dagb=0.
+  REAL(SP), ALLOCATABLE :: flux_dagb(:,:,:,:)
   REAL(SP), DIMENSION(2,ntau_dagb)                  :: tau1_dagb=0.
   REAL(SP), DIMENSION(2,nteff_dagb)                 :: teff_dagb=0.
 
   !nebular emission model
   REAL(SP), DIMENSION(nemline) :: nebem_line_pos=0.
   REAL(SP), DIMENSION(nemline,nebnz,nebnage,nebnip) :: nebem_line=0.,xnebem_line=0.
-  REAL(SP), DIMENSION(nspec,nebnz,nebnage,nebnip) :: nebem_cont=0.,xnebem_cont=0.
+  REAL(SP), ALLOCATABLE :: nebem_cont(:,:,:,:),xnebem_cont(:,:,:,:)
   REAL(SP), DIMENSION(nebnz)   :: nebem_logz=0.
   REAL(SP), DIMENSION(nebnage) :: nebem_age=0.
   REAL(SP), DIMENSION(nebnip)  :: nebem_logu=0.
   !minimum resolution for nebular lines, based
   !on the resolution of the spectral libraries.
-  REAL(SP), DIMENSION(nspec)   :: neb_res_min=0.0
-  REAL(SP), DIMENSION(nspec,nemline) :: gaussnebarr=0.0
+  REAL(SP), ALLOCATABLE   :: neb_res_min(:)
+  REAL(SP), ALLOCATABLE :: gaussnebarr(:,:)
 
   !arrays for AGN dust
   REAL(SP), DIMENSION(nagndust)       :: agndust_tau=0.
-  REAL(SP), DIMENSION(nspec,nagndust) :: agndust_spec=0.
+  REAL(SP), ALLOCATABLE :: agndust_spec(:,:)
 
   !arrays for the isochrone data
-  REAL(SP), DIMENSION(nz,nt,nm) :: mact_isoc=0.,logl_isoc=0.,&
-       logt_isoc=0.,logg_isoc=0.,ffco_isoc=0.,phase_isoc=0.,&
-       mini_isoc=0.,lmdot_isoc=0.
+  REAL(SP), ALLOCATABLE :: mact_isoc(:,:,:),logl_isoc(:,:,:),&
+       logt_isoc(:,:,:),logg_isoc(:,:,:),ffco_isoc(:,:,:),phase_isoc(:,:,:),&
+       mini_isoc(:,:,:),lmdot_isoc(:,:,:)
 
   !arrays holding the number of mass elements for each isochrone,
   !the age of each isochrone, and the metallicity of each isochrone
-  INTEGER, DIMENSION(nz,nt)  :: nmass_isoc=0
-  REAL(SP), DIMENSION(nz,nt) :: timestep_isoc=0.
-  REAL(SP), DIMENSION(nz)    :: zlegend=-99.
-  REAL(SP), DIMENSION(nzinit):: zlegendinit=-99.
+  INTEGER, ALLOCATABLE  :: nmass_isoc(:,:)
+  REAL(SP), ALLOCATABLE :: timestep_isoc(:,:)
+  REAL(SP), ALLOCATABLE    :: zlegend(:)
+  REAL(SP), ALLOCATABLE :: zlegendinit(:)
 
   !arrays for the full Z-dep SSP spectra
-  REAL(SP), DIMENSION(nspec,ntfull,nz) :: spec_ssp_zz=0.
-  REAL(SP), DIMENSION(ntfull,nz)       :: mass_ssp_zz=0.,lbol_ssp_zz=0.
-  REAL(SP), DIMENSION(ntfull)          :: time_full=0.
+  REAL(SP), ALLOCATABLE :: spec_ssp_zz(:,:,:)
+  REAL(SP), ALLOCATABLE :: mass_ssp_zz(:,:),lbol_ssp_zz(:,:)
+  REAL(SP), ALLOCATABLE :: time_full(:)
 
   !array for ssp weights
-  REAL(SP), DIMENSION(ntfull,nz)       :: weight_ssp=0.
+  REAL(SP), ALLOCATABLE :: weight_ssp(:,:)
 
   !array for young and old ages
-  REAL(SP), DIMENSION(nspec) :: spec_young=0.,spec_old=0.
+  REAL(SP), ALLOCATABLE :: spec_young(:),spec_old(:)
 
   !array for full BPASS SSPs
-  REAL(SP), DIMENSION(nspec,nt,nz) :: bpass_spec_ssp=0.
-  REAL(SP), DIMENSION(nt,nz)       :: bpass_mass_ssp=0.
+  REAL(SP), ALLOCATABLE :: bpass_spec_ssp(:,:,:)
+  REAL(SP), ALLOCATABLE :: bpass_mass_ssp(:,:)
 
   !arrays for X-ray binaries
-  REAL(SP), DIMENSION(nspec_xrb) :: lam_xrb=0.
-  REAL(SP), DIMENSION(nspec,nt_xrb,nz_xrb) :: spec_xrb=0.
-  REAL(SP), DIMENSION(nt_xrb) :: ages_xrb=0.0
-  REAL(SP), DIMENSION(nz_xrb) :: zmet_xrb=0.0
+  REAL(SP), ALLOCATABLE :: lam_xrb(:)
+  REAL(SP), ALLOCATABLE :: spec_xrb(:,:,:)
+  REAL(SP), ALLOCATABLE :: ages_xrb(:)
+  REAL(SP), ALLOCATABLE :: zmet_xrb(:)
   
   !------------Define TYPE structures-------------!
 
@@ -569,18 +460,18 @@ MODULE SPS_VARS
           max_wave_smooth=1E4,gas_logu=-2.0,gas_logz=0.,igm_factor=1.0,&
           fagn=0.0,agn_tau=10.0,frac_xrb=1.0,dust3=0.
      INTEGER :: zmet=1,sfh=0,wgp1=1,wgp2=1,wgp3=1,evtype=-1
-     INTEGER, DIMENSION(nbands) :: mag_compute=1
-     INTEGER, DIMENSION(nt) :: ssp_gen_age=1
+     INTEGER, ALLOCATABLE :: mag_compute(:)
+     INTEGER, ALLOCATABLE :: ssp_gen_age(:)
      CHARACTER(50) :: imf_filename='', sfh_filename=''
   END TYPE PARAMS
 
   !structure for the output of the compsp routine
   TYPE COMPSPOUT
      REAL(SP) :: age=0.,mass_csp=0.,lbol_csp=0.,sfr=0.,mdust=0.,mformed=0.
-     REAL(SP), DIMENSION(nbands)  :: mags=0.
-     REAL(SP), DIMENSION(nspec)   :: spec=0.
-     REAL(SP), DIMENSION(nindx)   :: indx=0.
-     REAL(SP), DIMENSION(nemline) :: emlines=0.
+     REAL(SP), ALLOCATABLE  :: mags(:)
+     REAL(SP), ALLOCATABLE   :: spec(:)
+     REAL(SP), ALLOCATABLE   :: indx(:)
+     REAL(SP), ALLOCATABLE :: emlines(:)
   END TYPE COMPSPOUT
 
   ! A structure to hold SFH params converted to intrinsic units
@@ -591,7 +482,7 @@ MODULE SPS_VARS
   END TYPE SFHPARAMS
 
   TYPE TLSF
-     REAL(SP), DIMENSION(nspec) :: lsf=0.
+     REAL(SP), ALLOCATABLE :: lsf(:)
      REAL(SP) :: minlam=0.,maxlam=0.
   END TYPE TLSF
 
@@ -603,8 +494,8 @@ MODULE SPS_VARS
   !structure for observational data
   TYPE OBSDAT
      REAL(SP)                    :: zred=0.,logsmass=0.
-     REAL(SP), DIMENSION(nbands) :: mags=0.,magerr=0.
-     REAL(SP), DIMENSION(nspec)  :: spec=0.,specerr=99.
+     REAL(SP), ALLOCATABLE :: mags(:),magerr(:)
+     REAL(SP), ALLOCATABLE  :: spec(:),specerr(:)
   END TYPE OBSDAT
 
   !structure for using P(z) in chi2
@@ -616,12 +507,4 @@ MODULE SPS_VARS
   !used for Powell minimization
   TYPE(OBSDAT) :: powell_data, sedfit_data
 
-  !used for creating a pre-tabulated grid of CSPs as a function
-  !of tau and metallicity
-  !INTEGER, PARAMETER :: ntaugrid=20
-  !REAL(SP), DIMENSION(ntaugrid) :: taugrid=0.0
-  !REAL, DIMENSION(nspec,ntfull,ntaugrid,nz) :: csp_grid=0.0
-  !INTEGER :: csp_grid_flag=0
-
-
 END MODULE SPS_VARS
diff --git a/tests/data/sps_ref_MILES-0_C3K-1_MIST-1.bin b/tests/data/sps_ref_MILES-0_C3K-1_MIST-1.bin
new file mode 100644
index 00000000..cb800db4
Binary files /dev/null and b/tests/data/sps_ref_MILES-0_C3K-1_MIST-1.bin differ
diff --git a/tests/data/sps_ref_MILES-0_MIST-0_BASEL-1_PADOVA-1.bin b/tests/data/sps_ref_MILES-0_MIST-0_BASEL-1_PADOVA-1.bin
new file mode 100644
index 00000000..96841d43
Binary files /dev/null and b/tests/data/sps_ref_MILES-0_MIST-0_BASEL-1_PADOVA-1.bin differ
diff --git a/tests/data/sps_ref_MILES-1_MIST-1.bin b/tests/data/sps_ref_MILES-1_MIST-1.bin
new file mode 100644
index 00000000..b80280cf
Binary files /dev/null and b/tests/data/sps_ref_MILES-1_MIST-1.bin differ
diff --git a/tests/data/sps_ref_MILES-1_MIST-1_THEMIS-1_DL07-0.bin b/tests/data/sps_ref_MILES-1_MIST-1_THEMIS-1_DL07-0.bin
new file mode 100644
index 00000000..cd021167
Binary files /dev/null and b/tests/data/sps_ref_MILES-1_MIST-1_THEMIS-1_DL07-0.bin differ
diff --git a/tests/data/sps_ref_MIST-0_BPASS-1.bin b/tests/data/sps_ref_MIST-0_BPASS-1.bin
new file mode 100644
index 00000000..9d8f373e
Binary files /dev/null and b/tests/data/sps_ref_MIST-0_BPASS-1.bin differ
diff --git a/tests/generate_test_data.f90 b/tests/generate_test_data.f90
new file mode 100644
index 00000000..99dac642
--- /dev/null
+++ b/tests/generate_test_data.f90
@@ -0,0 +1,199 @@
+PROGRAM GENERATE_TEST_DATA
+
+  ! Generates reference data for FSPS regression testing.
+  ! Uses allocatable arrays to support multiple compile-time configurations.
+
+  USE sps_vars
+  USE sps_utils
+  IMPLICIT NONE
+
+  ! Variables for SSP generation (allocatable)
+  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: spec_ssp
+  REAL(SP), ALLOCATABLE, DIMENSION(:)   :: mass_ssp, lbol_ssp
+  
+  ! Variables for CSP generation (allocatable)
+  TYPE(COMPSPOUT), ALLOCATABLE, DIMENSION(:) :: ocompsp
+  
+  ! Control variables
+  TYPE(PARAMS) :: pset
+  INTEGER :: i, unit_out, status, arg_count
+  CHARACTER(LEN=50) :: filename_out
+  CHARACTER(LEN=100) :: csp_dummy_file
+  CHARACTER(LEN=255) :: arg_val
+  CHARACTER(LEN=20) :: isoc_arg, spec_arg, dust_arg
+  LOGICAL :: isoc_set, spec_set, dust_set
+
+  ! Exit codes
+  INTEGER, PARAMETER :: EXIT_FAILURE = 1
+
+  ! Setup and initialization
+  
+  ! Name of the reference output file
+  filename_out = 'sps_test_output.bin'
+  unit_out = 40
+
+  csp_dummy_file = 'dummy_csp.out'
+  isoc_set = .FALSE.
+  spec_set = .FALSE.
+  dust_set = .FALSE.
+
+  WRITE(*,*) '--------------------------------------------------'
+  WRITE(*,*) 'FSPS REFERENCE DATA GENERATOR'
+
+  ! Parse command line arguments
+  arg_count = COMMAND_ARGUMENT_COUNT()
+  i = 1
+  DO WHILE (i <= arg_count)
+     CALL GET_COMMAND_ARGUMENT(i, arg_val, STATUS=status)
+     IF (status /= 0) EXIT
+     
+     IF (TRIM(arg_val) == '--isoc') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, isoc_arg, STATUS=status)
+           isoc_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --isoc requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     ELSE IF (TRIM(arg_val) == '--spec') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, spec_arg, STATUS=status)
+           spec_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --spec requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     ELSE IF (TRIM(arg_val) == '--dust') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, dust_arg, STATUS=status)
+           dust_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --dust requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     END IF
+     i = i + 1
+  END DO
+  
+  ! Initialize FSPS parameters (MIST/MILES defaults)
+  ! We call this FIRST so we can be sure parameters like ntfull/nspec are set
+  ! before we allocate (though they are static in the current codebase).
+  imf_type = 1          
+  pset%zmet = 10        
+  
+  ! Use optional arguments if set
+  IF (isoc_set .AND. spec_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), spec_type_in=TRIM(spec_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (isoc_set .AND. spec_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), spec_type_in=TRIM(spec_arg))
+  ELSE IF (isoc_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (spec_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, spec_type_in=TRIM(spec_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (isoc_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg))
+  ELSE IF (spec_set) THEN
+      CALL SPS_SETUP(pset%zmet, spec_type_in=TRIM(spec_arg))
+  ELSE IF (dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, dust_type_in=TRIM(dust_arg))
+  ELSE
+      CALL SPS_SETUP(pset%zmet)
+  END IF
+
+  ! Memory allocation
+  WRITE(*,*) 'Allocating memory (nspec:', nspec, ' ntfull:', ntfull, ')...'
+  ALLOCATE(spec_ssp(ntfull,nspec))
+  ALLOCATE(mass_ssp(ntfull))
+  ALLOCATE(lbol_ssp(ntfull))
+  ALLOCATE(ocompsp(ntfull))
+
+  ! Write header
+  WRITE(*,*) 'Writing to file: ', TRIM(filename_out)
+  OPEN(UNIT=unit_out, FILE=TRIM(filename_out), STATUS='REPLACE', &
+       FORM='UNFORMATTED', ACCESS='STREAM')
+
+  WRITE(*,*) 'Writing Header...'
+  WRITE(unit_out) nspec
+  WRITE(unit_out) ntfull
+  WRITE(unit_out) nbands
+
+  ! Test Case 1: Simple SSP (Solar, Chabrier)
+  WRITE(*,*) 'Running Test Case 1: SSP (Solar, Chabrier)...'
+
+  ! Define SSP Parameters
+  pset%sfh   = 0     ! SSP
+  pset%const = 0.0
+  pset%zred  = 0.0
+  pset%dust1 = 0.0
+  pset%dust2 = 0.0
+  add_neb_emission = 1 
+  
+  ! Allocate pset allocatable components
+  IF (ALLOCATED(pset%mag_compute)) DEALLOCATE(pset%mag_compute)
+  ALLOCATE(pset%mag_compute(nbands))
+  pset%mag_compute = 1
+  
+  IF (ALLOCATED(pset%ssp_gen_age)) DEALLOCATE(pset%ssp_gen_age)
+  ALLOCATE(pset%ssp_gen_age(nt))
+  pset%ssp_gen_age = 1
+
+  ! Compute SSP
+  CALL SSP_GEN(pset, mass_ssp, lbol_ssp, spec_ssp)
+
+  ! Write SSP Data
+  WRITE(*,*) 'Saving SSP results...'
+  WRITE(unit_out) mass_ssp
+  WRITE(unit_out) lbol_ssp
+  WRITE(unit_out) spec_ssp
+
+  ! Test Case 2: CSP (Tau Model, Dusty)
+  WRITE(*,*) 'Running Test Case 2: CSP (Tau=2.0, Dust=1.0)...'
+
+  pset%sfh   = 1     ! Tau model
+  pset%tau   = 2.0   
+  pset%dust1 = 1.0   
+  pset%dust2 = 0.3
+  
+  ! Re-run SSP_GEN
+  CALL SSP_GEN(pset, mass_ssp, lbol_ssp, spec_ssp)
+
+  ! Manually allocate components of ocompsp array elements
+  DO i = 1, ntfull
+     IF (.NOT. ALLOCATED(ocompsp(i)%mags)) ALLOCATE(ocompsp(i)%mags(nbands))
+     IF (.NOT. ALLOCATED(ocompsp(i)%spec)) ALLOCATE(ocompsp(i)%spec(nspec))
+     IF (.NOT. ALLOCATED(ocompsp(i)%indx)) ALLOCATE(ocompsp(i)%indx(nindx))
+     IF (.NOT. ALLOCATED(ocompsp(i)%emlines)) ALLOCATE(ocompsp(i)%emlines(nemline))
+  END DO
+
+  ! Compute CSP
+  CALL COMPSP(3, 1, csp_dummy_file, mass_ssp, lbol_ssp, spec_ssp, pset, ocompsp)
+
+  ! Write CSP Data
+  WRITE(*,*) 'Saving CSP results...'
+  DO i = 1, ntfull
+     WRITE(unit_out) ocompsp(i)%age
+     WRITE(unit_out) ocompsp(i)%mass_csp
+     WRITE(unit_out) ocompsp(i)%lbol_csp
+     WRITE(unit_out) ocompsp(i)%sfr
+     WRITE(unit_out) ocompsp(i)%mdust
+     WRITE(unit_out) ocompsp(i)%mformed
+     WRITE(unit_out) ocompsp(i)%mags      
+     WRITE(unit_out) ocompsp(i)%spec      
+     WRITE(unit_out) ocompsp(i)%indx      
+     WRITE(unit_out) ocompsp(i)%emlines   
+  END DO
+
+  ! Cleanup
+  CLOSE(unit_out)
+  DEALLOCATE(spec_ssp, mass_ssp, lbol_ssp, ocompsp)
+  IF (ALLOCATED(pset%mag_compute)) DEALLOCATE(pset%mag_compute)
+  IF (ALLOCATED(pset%ssp_gen_age)) DEALLOCATE(pset%ssp_gen_age)
+  
+  CALL SPS_TAKEDOWN()
+
+  WRITE(*,*) 'Complete. Data saved.'
+
+END PROGRAM GENERATE_TEST_DATA
diff --git a/tests/generate_test_data.sh b/tests/generate_test_data.sh
new file mode 100755
index 00000000..6786cb03
--- /dev/null
+++ b/tests/generate_test_data.sh
@@ -0,0 +1,72 @@
+#!/bin/bash
+
+# Generate FSPS test reference data for multiple library configurations.
+# This script must be run from the tests/ directory.
+# Usage: ./generate_test_data.sh
+
+set -e
+
+# Validate we're in the correct directory
+if [[ ! -f "generate_test_data.sh" ]]; then
+    echo "ERROR: This script must be run from the tests/ directory"
+    exit 1
+fi
+
+# Define the library combinations
+# Format: "Legacy_File_Suffix|Runtime_Args"
+declare -a configurations=(
+    "MILES-1_MIST-1|--isoc mist --spec miles"
+    "MILES-0_MIST-0_BASEL-1_PADOVA-1|--isoc pdva --spec basel"
+    "MILES-1_MIST-1_THEMIS-1_DL07-0|--isoc mist --spec miles --dust themis"
+    "MILES-0_C3K-1_MIST-1|--isoc mist --spec c3k_afe+0.0"
+    "MIST-0_BPASS-1|--isoc bpss --spec bpass"
+)
+
+# Create the data directory if it doesn't exist
+mkdir -p data
+
+# Move into root directory to run Make
+cd ..
+
+echo "=========================================================="
+echo "Compiling FSPS objects and generator..."
+echo "=========================================================="
+
+# Clean previous build artifacts
+make clean > /dev/null 2>&1
+
+# Build the generator target ONCE
+make generate_test_data
+
+if [ ! -f ./generate_test_data ]; then
+    echo "ERROR: Compilation failed."
+    exit 1
+fi
+
+# Main Loop
+for config in "${configurations[@]}"; do
+    IFS="|" read -r suffix args <<< "$config"
+    
+    echo "=========================================================="
+    echo "Running configuration: $args"
+    echo "Output: tests/data/sps_ref_${suffix}.bin"
+    echo "=========================================================="
+
+    # Run the generator
+    ./generate_test_data $args
+    
+    # Move and rename output to the tests/data folder
+    if [ -f "sps_test_output.bin" ]; then
+        mv sps_test_output.bin "tests/data/sps_ref_${suffix}.bin"
+        echo "Created: tests/data/sps_ref_${suffix}.bin"
+    else
+        echo "ERROR: Output file not generated for $args"
+        exit 1
+    fi
+    
+    echo ""
+done
+
+# Cleanup the executable from src
+rm -f generate_test_data
+echo "Done."
diff --git a/tests/test_runner.f90 b/tests/test_runner.f90
new file mode 100644
index 00000000..80f4eda7
--- /dev/null
+++ b/tests/test_runner.f90
@@ -0,0 +1,367 @@
+PROGRAM TEST_RUNNER
+
+  ! FSPS Regression Test Runner
+  ! Reads reference data from disk, runs current FSPS code with identical
+  ! parameters, and compares outputs with a relative tolerance.
+
+  USE, INTRINSIC :: IEEE_ARITHMETIC
+  USE sps_vars
+  USE sps_utils
+  IMPLICIT NONE
+
+  ! Exit codes
+  INTEGER, PARAMETER :: EXIT_SUCCESS = 0
+  INTEGER, PARAMETER :: EXIT_FAILURE = 1
+  
+  ! Default tolerance (can be overridden by env var FSPS_TEST_RTOL)
+  REAL(SP), PARAMETER :: DEFAULT_RTOL = 1.0E-5
+
+  ! Test arrays (allocatable)
+  ! Reference data (read from disk)
+  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: ref_spec_ssp
+  REAL(SP), ALLOCATABLE, DIMENSION(:)   :: ref_mass_ssp, ref_lbol_ssp
+  TYPE(COMPSPOUT), ALLOCATABLE, DIMENSION(:) :: ref_ocompsp
+
+  ! New data (computed on the fly)
+  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: new_spec_ssp
+  REAL(SP), ALLOCATABLE, DIMENSION(:)   :: new_mass_ssp, new_lbol_ssp
+  TYPE(COMPSPOUT), ALLOCATABLE, DIMENSION(:) :: new_ocompsp
+
+  ! Control variables
+  TYPE(PARAMS) :: pset
+  INTEGER :: i, j, unit_in, status, arg_count
+  CHARACTER(LEN=255) :: filename_in, env_buffer, arg_val
+  CHARACTER(LEN=100) :: csp_dummy_file
+  CHARACTER(LEN=20) :: isoc_arg, spec_arg, dust_arg
+  LOGICAL :: isoc_set, spec_set, dust_set
+  
+  ! Dimensions read from file
+  INTEGER :: file_nspec, file_ntfull, file_nbands
+
+  ! Comparison stats
+  REAL(SP) :: rtol
+  LOGICAL :: test_passed
+
+  ! Configuration
+  unit_in = 40
+  csp_dummy_file = 'dummy_csp.out'
+  test_passed = .TRUE.
+  isoc_set = .FALSE.
+  spec_set = .FALSE.
+  dust_set = .FALSE.
+  filename_in = ''
+
+  WRITE(*,*) '========================================='
+  WRITE(*,*) 'FSPS TEST RUNNER'
+  WRITE(*,*) '========================================='
+
+  ! Parse command line arguments
+  arg_count = COMMAND_ARGUMENT_COUNT()
+  i = 1
+  DO WHILE (i <= arg_count)
+     CALL GET_COMMAND_ARGUMENT(i, arg_val, STATUS=status)
+     IF (status /= 0) EXIT
+     
+     IF (TRIM(arg_val) == '--isoc') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, isoc_arg, STATUS=status)
+           isoc_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --isoc requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     ELSE IF (TRIM(arg_val) == '--spec') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, spec_arg, STATUS=status)
+           spec_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --spec requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     ELSE IF (TRIM(arg_val) == '--dust') THEN
+        i = i + 1
+        IF (i <= arg_count) THEN
+           CALL GET_COMMAND_ARGUMENT(i, dust_arg, STATUS=status)
+           dust_set = .TRUE.
+        ELSE
+           WRITE(*,*) 'ERROR: --dust requires an argument'
+           STOP EXIT_FAILURE
+        END IF
+     ELSE
+        ! Assume it's the filename if it doesn't start with --
+        ! Or if we haven't found one yet.
+        IF (LEN_TRIM(filename_in) == 0) THEN
+           filename_in = TRIM(arg_val)
+        END IF
+     END IF
+     i = i + 1
+  END DO
+
+  IF (LEN_TRIM(filename_in) == 0) THEN
+     WRITE(*,*) 'ERROR: Must provide reference filename as argument.'
+     WRITE(*,*) 'Usage: ./test_runner [--isoc type] [--spec type] [--dust type] tests/data/sps_ref_XXX.bin'
+     STOP EXIT_FAILURE
+  END IF
+
+  ! Read tolerance from environment variable
+  CALL GET_ENVIRONMENT_VARIABLE("FSPS_TEST_RTOL", VALUE=env_buffer, STATUS=status)
+  IF (status == 0) THEN
+     READ(env_buffer, *) rtol
+     WRITE(*,*) 'Using RTOL from environment: ', rtol
+  ELSE
+     rtol = DEFAULT_RTOL
+     WRITE(*,*) 'Using default RTOL: ', rtol
+  END IF
+
+  ! Initialize FSPS and check dimensions
+  ! Note: We must initialize FSPS before allocating, but we must read the 
+  ! file header before we know if dimensions match.
+  
+  imf_type = 1          
+  pset%zmet = 10
+  
+  WRITE(*,*) 'Initializing FSPS...'
+  IF (.NOT. isoc_set) isoc_arg = ''
+  IF (.NOT. spec_set) spec_arg = ''
+  IF (.NOT. dust_set) dust_arg = ''
+  
+  ! Use optional arguments if set, otherwise rely on defaults handled by SPS_SETUP logic
+  ! Since Fortran optional arguments must be passed if present in call, and we want 
+  ! to mix and match, it's cleaner to build call permutations or just pass blanks if allowed/handled.
+  ! However, sps_setup uses PRESENT(). 
+  
+  IF (isoc_set .AND. spec_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), spec_type_in=TRIM(spec_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (isoc_set .AND. spec_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), spec_type_in=TRIM(spec_arg))
+  ELSE IF (isoc_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (spec_set .AND. dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, spec_type_in=TRIM(spec_arg), dust_type_in=TRIM(dust_arg))
+  ELSE IF (isoc_set) THEN
+      CALL SPS_SETUP(pset%zmet, isoc_type_in=TRIM(isoc_arg))
+  ELSE IF (spec_set) THEN
+      CALL SPS_SETUP(pset%zmet, spec_type_in=TRIM(spec_arg))
+  ELSE IF (dust_set) THEN
+      CALL SPS_SETUP(pset%zmet, dust_type_in=TRIM(dust_arg))
+  ELSE
+      CALL SPS_SETUP(pset%zmet)
+  END IF
+
+  ! Allocate pset allocatable components
+  IF (ALLOCATED(pset%mag_compute)) DEALLOCATE(pset%mag_compute)
+  ALLOCATE(pset%mag_compute(nbands))
+  pset%mag_compute = 1
+  
+  IF (ALLOCATED(pset%ssp_gen_age)) DEALLOCATE(pset%ssp_gen_age)
+  ALLOCATE(pset%ssp_gen_age(nt))
+  pset%ssp_gen_age = 1
+
+  ! Open the reference file
+  OPEN(UNIT=unit_in, FILE=TRIM(filename_in), STATUS='OLD', &
+       FORM='UNFORMATTED', ACCESS='STREAM', IOSTAT=status)
+  IF (status /= 0) THEN
+     WRITE(*,*) 'ERROR: Could not open file: ', TRIM(filename_in)
+     STOP EXIT_FAILURE
+  END IF
+
+  ! Read Header
+  READ(unit_in) file_nspec
+  READ(unit_in) file_ntfull
+  READ(unit_in) file_nbands
+
+  WRITE(*,*) 'Reference Dimensions: nspec=', file_nspec, ' nt=', file_ntfull
+  WRITE(*,*) 'Compiled Dimensions:  nspec=', nspec,      ' nt=', ntfull
+
+  ! Strict dimension check
+  IF (file_nspec /= nspec .OR. file_ntfull /= ntfull) THEN
+     WRITE(*,*) 'FATAL: Binary dimensions do not match compiled FSPS dimensions.'
+     WRITE(*,*) 'Ensure you are running the test with the same flags/arguments used to generate the data.'
+     STOP EXIT_FAILURE
+  END IF
+
+  ! Allocate and read reference data
+  ALLOCATE(ref_spec_ssp(ntfull,nspec))
+  ALLOCATE(ref_mass_ssp(ntfull))
+  ALLOCATE(ref_lbol_ssp(ntfull))
+  ALLOCATE(ref_ocompsp(ntfull))
+
+  ALLOCATE(new_spec_ssp(ntfull,nspec))
+  ALLOCATE(new_mass_ssp(ntfull))
+  ALLOCATE(new_lbol_ssp(ntfull))
+  ALLOCATE(new_ocompsp(ntfull))
+
+  WRITE(*,*) 'Reading SSP reference data...'
+  READ(unit_in) ref_mass_ssp
+  READ(unit_in) ref_lbol_ssp
+  READ(unit_in) ref_spec_ssp
+
+  WRITE(*,*) 'Reading CSP reference data...'
+  DO i = 1, ntfull
+     ! Allocate components of derived type before reading
+     ALLOCATE(ref_ocompsp(i)%mags(nbands))
+     ALLOCATE(ref_ocompsp(i)%spec(nspec))
+     ALLOCATE(ref_ocompsp(i)%indx(nindx))
+     ALLOCATE(ref_ocompsp(i)%emlines(nemline))
+     
+     READ(unit_in) ref_ocompsp(i)%age
+     READ(unit_in) ref_ocompsp(i)%mass_csp
+     READ(unit_in) ref_ocompsp(i)%lbol_csp
+     READ(unit_in) ref_ocompsp(i)%sfr
+     READ(unit_in) ref_ocompsp(i)%mdust
+     READ(unit_in) ref_ocompsp(i)%mformed
+     READ(unit_in) ref_ocompsp(i)%mags      
+     READ(unit_in) ref_ocompsp(i)%spec      
+     READ(unit_in) ref_ocompsp(i)%indx      
+     READ(unit_in) ref_ocompsp(i)%emlines   
+  END DO
+  CLOSE(unit_in)
+
+  ! Generate new data
+  WRITE(*,*) 'Generating new SSP data...'
+  pset%sfh   = 0     
+  pset%const = 0.0
+  pset%zred  = 0.0
+  pset%dust1 = 0.0
+  pset%dust2 = 0.0
+  add_neb_emission = 1 
+  CALL SSP_GEN(pset, new_mass_ssp, new_lbol_ssp, new_spec_ssp)
+
+  WRITE(*,*) 'Generating new CSP data...'
+  pset%sfh   = 1    
+  pset%tau   = 2.0   
+  pset%dust1 = 1.0   
+  pset%dust2 = 0.3
+  CALL SSP_GEN(pset, new_mass_ssp, new_lbol_ssp, new_spec_ssp)
+  
+  DO i = 1, ntfull
+     IF (.NOT. ALLOCATED(new_ocompsp(i)%mags)) ALLOCATE(new_ocompsp(i)%mags(nbands))
+     IF (.NOT. ALLOCATED(new_ocompsp(i)%spec)) ALLOCATE(new_ocompsp(i)%spec(nspec))
+     IF (.NOT. ALLOCATED(new_ocompsp(i)%indx)) ALLOCATE(new_ocompsp(i)%indx(nindx))
+     IF (.NOT. ALLOCATED(new_ocompsp(i)%emlines)) ALLOCATE(new_ocompsp(i)%emlines(nemline))
+  END DO
+
+  CALL COMPSP(3, 1, csp_dummy_file, new_mass_ssp, new_lbol_ssp, new_spec_ssp, pset, new_ocompsp)
+
+  ! Compare results
+  WRITE(*,*) 'Verifying results (RTOL = ', rtol, ')...'
+
+  ! Helper internal subroutine to check arrays
+  CALL CHECK_ARRAY_2D("SSP Spectra", ref_spec_ssp, new_spec_ssp, ntfull, nspec)
+  CALL CHECK_ARRAY_1D("SSP Mass", ref_mass_ssp, new_mass_ssp, ntfull)
+  CALL CHECK_ARRAY_1D("SSP Lbol", ref_lbol_ssp, new_lbol_ssp, ntfull)
+
+  ! Check CSP structure components manually
+  DO i = 1, ntfull
+     ! Check Scalars
+     CALL CHECK_VAL("CSP Lbol", i, ref_ocompsp(i)%lbol_csp, new_ocompsp(i)%lbol_csp)
+     CALL CHECK_VAL("CSP Mass", i, ref_ocompsp(i)%mass_csp, new_ocompsp(i)%mass_csp)
+     CALL CHECK_VAL("CSP SFR", i, ref_ocompsp(i)%sfr, new_ocompsp(i)%sfr)
+     CALL CHECK_VAL("CSP Dust Mass", i, ref_ocompsp(i)%mdust, new_ocompsp(i)%mdust)
+     CALL CHECK_VAL("CSP Mass Formed", i, ref_ocompsp(i)%mformed, new_ocompsp(i)%mformed)
+
+     ! Check Arrays for EVERY time step
+     CALL CHECK_ARRAY_1D("CSP Mags", ref_ocompsp(i)%mags, new_ocompsp(i)%mags, nbands)
+     CALL CHECK_ARRAY_1D("CSP Indx", ref_ocompsp(i)%indx, new_ocompsp(i)%indx, nindx)
+     CALL CHECK_ARRAY_1D("CSP Emlines", ref_ocompsp(i)%emlines, new_ocompsp(i)%emlines, nemline)
+     
+     ! Check Spectrum
+     CALL CHECK_ARRAY_1D("CSP Spec", ref_ocompsp(i)%spec, new_ocompsp(i)%spec, nspec)
+  END DO
+
+  ! Report results
+  WRITE(*,*) '--------------------------------------------------'  
+
+  CALL SPS_TAKEDOWN()
+
+  IF (test_passed) THEN
+     WRITE(*,*) 'TEST RESULT: PASS'
+     STOP EXIT_SUCCESS
+  ELSE
+     WRITE(*,*) 'TEST RESULT: FAIL'
+     STOP EXIT_FAILURE
+  END IF
+
+
+CONTAINS
+
+  SUBROUTINE CHECK_ARRAY_2D(label, ref, new, d1, d2)
+    CHARACTER(*), INTENT(IN) :: label
+    INTEGER, INTENT(IN) :: d1, d2
+    REAL(SP), DIMENSION(d1,d2), INTENT(IN) :: ref, new
+    REAL(SP) :: delta, threshold
+    INTEGER :: j, k
+
+    DO k = 1, d2
+       DO j = 1, d1
+          IF (IEEE_IS_NAN(ref(j,k)) .OR. IEEE_IS_NAN(new(j,k))) THEN
+             WRITE(*,*) 'FAIL: ', label, ' contains NaN at index (',j,',',k,')'
+             test_passed = .FALSE.
+             RETURN 
+          END IF
+
+          delta = ABS(ref(j,k) - new(j,k))
+          ! If ref is close to zero, use absolute tolerance, else relative
+          threshold = MAX(ABS(ref(j,k)) * rtol, 1.0E-30) 
+          
+          IF (delta > threshold) THEN
+             WRITE(*,*) 'FAIL: ', label, ' mismatch at index (',j,',',k,')'
+             WRITE(*,*) '  Ref:', ref(j,k), ' New:', new(j,k), ' Diff:', delta
+             test_passed = .FALSE.
+             RETURN ! Return early on failure to avoid log spam
+          END IF
+       END DO
+    END DO
+  END SUBROUTINE CHECK_ARRAY_2D
+
+  SUBROUTINE CHECK_ARRAY_1D(label, ref, new, d1)
+    CHARACTER(*), INTENT(IN) :: label
+    INTEGER, INTENT(IN) :: d1
+    REAL(SP), DIMENSION(d1), INTENT(IN) :: ref, new
+    REAL(SP) :: delta, threshold
+    INTEGER :: j
+
+    DO j = 1, d1
+       IF (IEEE_IS_NAN(ref(j)) .OR. IEEE_IS_NAN(new(j))) THEN
+          WRITE(*,*) 'FAIL: ', label, ' contains NaN at index (',j,')'
+          test_passed = .FALSE.
+          RETURN
+       END IF
+       
+       delta = ABS(ref(j) - new(j))
+       threshold = MAX(ABS(ref(j)) * rtol, 1.0E-30)
+       
+       IF (delta > threshold) THEN
+          WRITE(*,*) 'FAIL: ', label, ' mismatch at index (',j,')'
+          WRITE(*,*) '  Ref:', ref(j), ' New:', new(j), ' Diff:', delta
+          test_passed = .FALSE.
+          RETURN
+       END IF
+    END DO
+  END SUBROUTINE CHECK_ARRAY_1D
+
+  SUBROUTINE CHECK_VAL(label, idx, r, n)
+    CHARACTER(*), INTENT(IN) :: label
+    INTEGER, INTENT(IN) :: idx
+    REAL(SP), INTENT(IN) :: r, n
+    REAL(SP) :: delta, threshold
+
+    IF (IEEE_IS_NAN(r) .OR. IEEE_IS_NAN(n)) THEN
+       WRITE(*,*) 'FAIL: ', label, ' contains NaN at step ', idx
+       test_passed = .FALSE.
+       RETURN
+    END IF
+
+    delta = ABS(r - n)
+    threshold = MAX(ABS(r) * rtol, 1.0E-30)
+
+    IF (delta > threshold) THEN
+       WRITE(*,*) 'FAIL: ', label, ' mismatch at step ', idx
+       WRITE(*,*) '  Ref:', r, ' New:', n, ' Diff:', delta
+       test_passed = .FALSE.
+    END IF
+  END SUBROUTINE CHECK_VAL
+
+END PROGRAM TEST_RUNNER