diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90
index a9eed78a..059dbcfd 100644
--- a/src/physics/grad_shaf_prof_phys.F90
+++ b/src/physics/grad_shaf_prof_phys.F90
@@ -22,7 +22,7 @@ module grad_shaf_prof_phys
   oft_increase_indent, oft_decrease_indent, gsinv_interp, gs_prof_interp, &
   gs_get_qprof, gs_ani_press, gs_epsilon
 use tracing_2d, only: set_tracer, active_tracer, tracinginv_fs
-use oft_gs_profiles, only: spline_flux_func, linterp_flux_func
+use oft_gs_profiles, only: spline_flux_func, linterp_flux_func, get_max_threads
 use spline_mod
 USE mhd_utils, ONLY: mu0
 implicit none
@@ -137,7 +137,7 @@ SUBROUTINE create_mercier_ff(func,npsi)
     self%npsi=npsi
     CALL spline_alloc(self%func,self%npsi,1)
     CALL spline_alloc(self%funcp,self%npsi,1)
-    DO i=1,omp_get_max_threads()
+    DO i=1,get_max_threads()
       CALL spline_alloc(self%fun_loc(i),self%npsi,1)
     END DO
 end select
@@ -252,7 +252,7 @@ subroutine mercier_update(self,gseq)
 self%func%xs=self%funcp%xs
 self%func%fs(:,1)=self%funcp%fsi(:,1)
 CALL spline_fit(self%func,"extrap")
-DO k=1,omp_get_max_threads()
+DO k=1,get_max_threads()
   CALL spline_copy(self%func,self%fun_loc(k))
 END DO
 !
@@ -451,7 +451,7 @@ SUBROUTINE create_dipole_b0_prof(func,npsi)
     !---
     self%npsi=npsi
     CALL spline_alloc(self%func,self%npsi,1)
-    DO i=1,omp_get_max_threads()
+    DO i=1,get_max_threads()
       CALL spline_alloc(self%fun_loc(i),self%npsi,1)
     END DO
 end select
@@ -542,7 +542,7 @@ subroutine dipole_b0_update(self,gseq)
 self%ypn=0.d0
 !---Setup Spline
 CALL spline_fit(self%func,"extrap")
-DO k=1,omp_get_max_threads()
+DO k=1,get_max_threads()
   CALL spline_copy(self%func,self%fun_loc(k))
 END DO
 IF(oft_debug_print(2))CALL oft_decrease_indent
@@ -616,7 +616,7 @@ subroutine dipole_ani_delete(self)
   SELECT TYPE(this=>self%B0_prof)
   TYPE IS(dipole_b0_flux_func)
     CALL spline_dealloc(this%func)
-    DO i=1,omp_get_max_threads()
+    DO i=1,get_max_threads()
       CALL spline_dealloc(this%fun_loc(i))
     END DO
   END SELECT
@@ -661,7 +661,7 @@ SUBROUTINE create_mirror_b0_prof(func,npsi)
     !---
     self%npsi=npsi
     CALL spline_alloc(self%func,self%npsi,1)
-    DO i=1,omp_get_max_threads()
+    DO i=1,get_max_threads()
       CALL spline_alloc(self%fun_loc(i),self%npsi,1)
     END DO
 end select
@@ -723,7 +723,7 @@ subroutine mirror_b0_update(self,gseq)
 self%xmax=self%func%xs(self%npsi)+(self%func%xs(1)-self%func%xs(0))*4.d0
 !---Setup Spline
 CALL spline_fit(self%func,"extrap")
-DO k=1,omp_get_max_threads()
+DO k=1,get_max_threads()
   CALL spline_copy(self%func,self%fun_loc(k))
 END DO
 IF(oft_debug_print(2))CALL oft_decrease_indent
@@ -759,7 +759,7 @@ subroutine mirror_slosh_delete(self)
   SELECT TYPE(this=>self%B0_prof)
   TYPE IS(mirror_b0_flux_func)
     CALL spline_dealloc(this%func)
-    DO i=1,omp_get_max_threads()
+    DO i=1,get_max_threads()
       CALL spline_dealloc(this%fun_loc(i))
     END DO
   END SELECT
diff --git a/src/physics/grad_shaf_profiles.F90 b/src/physics/grad_shaf_profiles.F90
index 83b3d674..99efb212 100644
--- a/src/physics/grad_shaf_profiles.F90
+++ b/src/physics/grad_shaf_profiles.F90
@@ -89,6 +89,7 @@ module oft_gs_profiles
 !------------------------------------------------------------------------------
 !> Needs docs
 !------------------------------------------------------------------------------
+integer(4), parameter :: max_threads = 24
 type, extends(flux_func) :: spline_flux_func
   INTEGER(4) :: npsi = 0 !< Needs docs
   REAL(8) :: xmin = 0.d0 !< Needs docs
@@ -98,7 +99,7 @@ module oft_gs_profiles
   REAL(8) :: yp1 = 0.d0 !< Needs docs
   REAL(8) :: ypn = 0.d0 !< Needs docs
   TYPE(spline_type) :: func !< Needs docs
-  TYPE(spline_type) :: fun_loc(24) !< Needs docs
+  TYPE(spline_type) :: fun_loc(max_threads) !< Needs docs
 contains
   !> Needs docs
   procedure :: f => spline_f
@@ -551,7 +552,7 @@ SUBROUTINE create_spline_ff(func,npsi,psimin,psimax,psivals)
   self%npsi=npsi
   self%ncofs=self%npsi
   CALL spline_alloc(self%func,self%npsi-1,1)
-  DO i=1,omp_get_max_threads()
+  DO i=1,get_max_threads()
     CALL spline_alloc(self%fun_loc(i),self%npsi-1,1)
   END DO
   IF(PRESENT(psivals))THEN
@@ -578,7 +579,7 @@ SUBROUTINE create_spline_ff(func,npsi,psimin,psimax,psivals)
   WRITE(*,*)'Fitting',self%func%xs
   !CALL spline_fit(self%func,"extrap")
   CALL spline_fit(self%func,"not-a-knot")
-  DO i=1,omp_get_max_threads()
+  DO i=1,get_max_threads()
     CALL spline_copy(self%func,self%fun_loc(i))
   END DO
   ALLOCATE(c(self%ncofs))
@@ -645,7 +646,7 @@ subroutine spline_update(self,gseq)
 CALL spline_eval(self%func,self%xmax,1)
 self%ypn=self%func%f1(1)
 !---
-DO i=1,omp_get_max_threads()
+DO i=1,get_max_threads()
   CALL spline_copy(self%func,self%fun_loc(i))
 END DO
 end subroutine spline_update
@@ -679,7 +680,7 @@ function spline_cofs_update(self,c) result(ierr)
 self%fn=self%func%f(1)-self%f1
 self%ypn=self%func%f1(1)
 !---
-DO i=1,omp_get_max_threads()
+DO i=1,get_max_threads()
   CALL spline_copy(self%func,self%fun_loc(i))
 END DO
 ierr=0
@@ -696,6 +697,14 @@ subroutine spline_cofs_get(self,c)
 END DO
 end subroutine spline_cofs_get
 !------------------------------------------------------------------------------
+!> Return the maximum number of threads that can be created, limited by the
+!> allocated size of the fun_loc array.
+!------------------------------------------------------------------------------
+function get_max_threads() result(threads)
+integer(4) :: threads
+threads = min(omp_get_max_threads(), max_threads)
+end function get_max_threads
+!------------------------------------------------------------------------------
 !> Needs docs
 !------------------------------------------------------------------------------
 SUBROUTINE create_linterp_ff(func,npsi,psivals,yvals,y0)