diff --git a/Makefile b/Makefile
index 46d620dae..714a72d4b 100644
--- a/Makefile
+++ b/Makefile
@@ -60,9 +60,9 @@ SWIFTEST_MODULES =   swiftest_globals.f90 \
 
 include Makefile.Defines
 
-MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
-IMKL = -I$(MKLROOT)/include
-LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
+#MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
+#IMKL = -I$(MKLROOT)/include
+#LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
 
 MODULES         = $(SWIFTEST_MODULES) $(USER_MODULES) 
 
diff --git a/Makefile.Defines b/Makefile.Defines
index 3789bb8a8..c85166569 100644
--- a/Makefile.Defines
+++ b/Makefile.Defines
@@ -64,18 +64,18 @@ GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none
 GPAR = -fopenmp #-ftree-parallelize-loops=4
 GMEM =  -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
-GPRODUCTION = -O2 -ffree-line-length-none $(GPAR)
+GPRODUCTION = -O3 -ffree-line-length-none $(GPAR)
 
 #FFLAGS  = $(IDEBUG) $(SIMDVEC) $(PAR)
 #FFASTFLAGS  = $(IDEBUG) $(SIMDVEC) $(PAR)
-FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT)  #$(ADVIXE_FLAGS)
-FFLAGS = $(IPRODUCTION) -fp-model=fast $(OPTREPORT)  #$(ADVIXE_FLAGS)
-FORTRAN = ifort
+#FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT)  #$(ADVIXE_FLAGS)
+#FFLAGS = $(IPRODUCTION) -fp-model=fast $(OPTREPORT)  #$(ADVIXE_FLAGS)
+#FORTRAN = ifort
 AR     = xiar
 
-#FORTRAN = gfortran
-#FFLAGS =  $(GDEBUG) # $(GMEM) $(GPAR)
-#FFLAGS =  $(GPRODUCTION) -g -fbacktrace #-fcheck=all #-Wall AR     = ar # DO NOT include in CFLAGS the "-c" option to compile object only
+FORTRAN = gfortran
+FFLAGS =  $(GDEBUG) # $(GMEM) $(GPAR)
+#FFLAGS =  $(GPRODUCTION) -g -fbacktrace -fcheck=all #-Wall AR     = ar # DO NOT include in CFLAGS the "-c" option to compile object only
 # this is done explicitly as needed in the Makefile
 
 CC              = icc
diff --git a/src/drift/drift.f90 b/src/drift/drift.f90
index 0437a18d1..dd53a40b3 100644
--- a/src/drift/drift.f90
+++ b/src/drift/drift.f90
@@ -88,7 +88,7 @@ end subroutine drift_all
 
 
    module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
-      !$omp declare simd(drift_one)
+      !!$omp declare simd(drift_one)
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
       !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient
@@ -119,7 +119,7 @@ end subroutine drift_one
 
 
    pure subroutine drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
-      !$omp declare simd(drift_dan)
+      !!$omp declare simd(drift_dan)
       !! author: David A. Minton
       !!
       !! Perform Kepler drift, solving Kepler's equation in appropriate variables
@@ -197,7 +197,7 @@ end subroutine drift_dan
 
 
    pure subroutine drift_kepmd(dm, es, ec, x, s, c)
-      !$omp declare simd(drift_kepmd)
+      !!$omp declare simd(drift_kepmd)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in difference form for an ellipse for small input dm and eccentricity
@@ -243,7 +243,7 @@ end subroutine drift_kepmd
 
 
    pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag)
-      !$omp declare simd(drift_kepu)
+      !!$omp declare simd(drift_kepu)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables
@@ -272,7 +272,7 @@ end subroutine drift_kepu
 
 
    pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f)
-      !$omp declare simd(drift_kepu_fchk)
+      !!$omp declare simd(drift_kepu_fchk)
       !! author: David A. Minton
       !!
       !! Computes the value of f, the function whose root we are trying to find in universal variables
@@ -303,7 +303,7 @@ end subroutine drift_kepu_fchk
 
 
    pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s)
-      !$omp declare simd(drift_kepu_guess)
+      !!$omp declare simd(drift_kepu_guess)
       !! author: David A. Minton
       !!
       !! Compute initial guess for solving Kepler's equation using universal variables
@@ -348,7 +348,7 @@ end subroutine drift_kepu_guess
 
 
    pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag)
-      !$omp declare simd(drift_kepu_lag)
+      !!$omp declare simd(drift_kepu_lag)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables using Laguerre's method
@@ -403,7 +403,7 @@ end subroutine drift_kepu_lag
 
 
    pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag)
-      !$omp declare simd(drift_kepu_new)
+      !!$omp declare simd(drift_kepu_new)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables using Newton's method
@@ -455,7 +455,7 @@ end subroutine drift_kepu_new
 
 
    pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag)
-      !$omp declare simd(drift_kepu_p3solve)
+      !!$omp declare simd(drift_kepu_p3solve)
       !! author: David A. Minton
       !!
       !! Computes real root of cubic involved in setting initial guess for solving Kepler's equation in universal variables
@@ -506,7 +506,7 @@ end subroutine drift_kepu_p3solve
    
 
    pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3)
-      !$omp declare simd(drift_kepu_stumpff)
+      !!$omp declare simd(drift_kepu_stumpff)
       !! author: David A. Minton
       !!
       !! Compute Stumpff functions needed for Kepler drift in universal variables
diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90
index 122473eaa..1cf71296f 100644
--- a/src/encounter/encounter_check.f90
+++ b/src/encounter/encounter_check.f90
@@ -770,7 +770,7 @@ end subroutine encounter_check_all_triangular_pltp
 
 
    module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr)
-      !$omp declare simd(encounter_check_one)
+      !!$omp declare simd(encounter_check_one)
       !! author: David A. Minton
       !!
       !! Determine whether a test particle and planet are having or will have an encounter within the next time step
@@ -1068,7 +1068,7 @@ end subroutine encounter_check_sweep_aabb_all_single_list
 
 
    pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
-      !$omp declare simd(encounter_check_sweep_aabb_one_double_list)
+      !!$omp declare simd(encounter_check_sweep_aabb_one_double_list)
       !! author: David A. Minton
       !!
       !! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only)
@@ -1101,7 +1101,7 @@ end subroutine encounter_check_sweep_aabb_one_double_list
  
 
    pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
-      !$omp declare simd(encounter_check_sweep_aabb_one_single_list)
+      !!$omp declare simd(encounter_check_sweep_aabb_one_single_list)
       !! author: David A. Minton
       !!
       !! Performs a sweep operation on a single body. Mutual encounters allowed (e.g. pl-pl)
diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
index f8d0d398d..0b0126702 100644
--- a/src/fraggle/fraggle_generate.f90
+++ b/src/fraggle/fraggle_generate.f90
@@ -547,21 +547,19 @@ function radial_objective_function(v_r_mag_input) result(fval)
             real(DP), dimension(frag%nbody)       :: kearr
             real(DP)                              :: keo, ke_radial, rotmag2, vmag2
       
-            associate(nfrag => frag%nbody)
-               allocate(v_shift, mold=frag%vb)
-               v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) 
-               !$omp do simd firstprivate(nfrag) lastprivate(rotmag2, vmag2) 
-               do i = 1,nfrag
-                  rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2
-                  vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
-                  kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) 
-               end do
-               !$omp end do simd
-               keo = 2 * frag%ke_budget - sum(kearr(:))
-               ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin
-               ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for
-               fval = (keo / (2 * ke_radial))**2
-            end associate
+            allocate(v_shift, mold=frag%vb)
+            v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) 
+            !$omp do simd firstprivate(frag)
+            do i = 1,frag%nbody
+               rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2
+               vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
+               kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) 
+            end do
+            !$omp end do simd
+            keo = 2 * frag%ke_budget - sum(kearr(:))
+            ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin
+            ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for
+            fval = (keo / (2 * ke_radial))**2
       
             return
          end function radial_objective_function
diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90
index 03ebfbd30..2b63e1408 100644
--- a/src/helio/helio_drift.f90
+++ b/src/helio/helio_drift.f90
@@ -76,7 +76,7 @@ end subroutine helio_drift_tp
 
 
    pure elemental subroutine helio_drift_linear_one(xhx, xhy, xhz, ptx, pty, ptz, dt)
-      !$omp declare simd(helio_drift_linear_one)
+      !!$omp declare simd(helio_drift_linear_one)
       implicit none
       real(DP), intent(inout) :: xhx, xhy, xhz
       real(DP), intent(in) :: ptx, pty, ptz, dt
diff --git a/src/kick/kick.f90 b/src/kick/kick.f90
index 717a96753..7208de656 100644
--- a/src/kick/kick.f90
+++ b/src/kick/kick.f90
@@ -108,7 +108,8 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
          zr = x(3, j) - x(3, i) 
          rji2 = xr**2 + yr**2 + zr**2
          rlim2 = (radius(i) + radius(j))**2
-         if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
+         if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
+                                 ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
       end do
       !$omp end parallel do 
      
@@ -156,7 +157,8 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius
             zr = x(3, j) - x(3, i) 
             rji2 = xr**2 + yr**2 + zr**2
             rlim2 = (radius(i) + radius(j))**2
-            if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
+            if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
+                                    ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
          end do
       end do
       !$omp end parallel do
@@ -210,7 +212,7 @@ end subroutine kick_getacch_int_all_tp
 
 
    module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
-      !$omp declare simd(kick_getacch_int_one_pl)
+      !!$omp declare simd(kick_getacch_int_one_pl)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies
@@ -242,7 +244,7 @@ end subroutine kick_getacch_int_one_pl
 
 
    module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
-      !$omp declare simd(kick_getacch_int_one_tp)
+      !!$omp declare simd(kick_getacch_int_one_tp)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair.
diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90
index ea9bb3c4a..09e2c3324 100644
--- a/src/modules/encounter_classes.f90
+++ b/src/modules/encounter_classes.f90
@@ -137,7 +137,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in
          logical,              dimension(:), allocatable, intent(out), optional :: lvdotr      !! Array indicating which bodies are approaching
       end subroutine encounter_check_collapse_ragged_list
 
-      module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
+      module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
          implicit none
          class(encounter_bounding_box_1D), intent(inout) :: self       !! Bounding box structure along a single dimension
          integer(I4B),                     intent(in)    :: n          !! Number of bodies with extents
diff --git a/src/operators/operator_cross.f90 b/src/operators/operator_cross.f90
index 99de730fe..17b9ff929 100644
--- a/src/operators/operator_cross.f90
+++ b/src/operators/operator_cross.f90
@@ -8,7 +8,7 @@
 contains
 
    module pure function operator_cross_sp(A, B) result(C)
-      !$omp declare simd(operator_cross_sp)
+      !!$omp declare simd(operator_cross_sp)
       implicit none
       real(SP), dimension(:), intent(in) :: A, B
       real(SP), dimension(3) :: C
@@ -19,7 +19,7 @@ module pure function operator_cross_sp(A, B) result(C)
    end function operator_cross_sp
 
    module pure function operator_cross_dp(A, B) result(C)
-      !$omp declare simd(operator_cross_dp)
+      !!$omp declare simd(operator_cross_dp)
       implicit none
       real(DP), dimension(:), intent(in) :: A, B
       real(DP), dimension(3) :: C
@@ -30,7 +30,7 @@ module pure function operator_cross_dp(A, B) result(C)
    end function operator_cross_dp
 
    module pure function operator_cross_qp(A, B) result(C)
-      !$omp declare simd(operator_cross_qp)
+      !!$omp declare simd(operator_cross_qp)
       implicit none
       real(QP), dimension(:), intent(in) :: A, B
       real(QP), dimension(3) :: C
@@ -41,7 +41,7 @@ module pure function operator_cross_qp(A, B) result(C)
    end function operator_cross_qp
 
    module pure function operator_cross_i1b(A, B) result(C)
-      !$omp declare simd(operator_cross_i1b)
+      !!$omp declare simd(operator_cross_i1b)
       implicit none
       integer(I1B), dimension(:), intent(in) :: A, B
       integer(I1B), dimension(3) :: C
@@ -52,7 +52,7 @@ module pure function operator_cross_i1b(A, B) result(C)
    end function operator_cross_i1b
 
    module pure function operator_cross_i2b(A, B) result(C)
-      !$omp declare simd(operator_cross_i2b)
+      !!$omp declare simd(operator_cross_i2b)
       implicit none
       integer(I2B), dimension(:), intent(in) :: A, B
       integer(I2B), dimension(3) :: C
@@ -63,7 +63,7 @@ module pure function operator_cross_i2b(A, B) result(C)
    end function operator_cross_i2b
 
    module pure function operator_cross_i4b(A, B) result(C)
-      !$omp declare simd(operator_cross_i4b)
+      !!$omp declare simd(operator_cross_i4b)
       implicit none
       integer(I4B), dimension(:), intent(in) :: A, B
       integer(I4B), dimension(3) :: C
@@ -74,7 +74,7 @@ module pure function operator_cross_i4b(A, B) result(C)
    end function operator_cross_i4b
 
    module pure function operator_cross_i8b(A, B) result(C)     
-      !$omp declare simd(operator_cross_i8b)
+      !!$omp declare simd(operator_cross_i8b)
       implicit none
       integer(I8B), dimension(:), intent(in) :: A, B
       integer(I8B), dimension(3) :: C
diff --git a/src/operators/operator_mag.f90 b/src/operators/operator_mag.f90
index f74555775..4056e696e 100644
--- a/src/operators/operator_mag.f90
+++ b/src/operators/operator_mag.f90
@@ -7,7 +7,7 @@
    contains
 
    module pure function operator_mag_sp(A) result(B)
-      !$omp declare simd(operator_mag_sp)
+      !!$omp declare simd(operator_mag_sp)
       implicit none
       real(SP), dimension(:), intent(in) :: A
       real(SP)                           :: B
@@ -16,7 +16,7 @@ module pure function operator_mag_sp(A) result(B)
    end function operator_mag_sp
 
    module pure function operator_mag_dp(A) result(B)
-      !$omp declare simd(operator_mag_dp)
+      !!$omp declare simd(operator_mag_dp)
       implicit none
       real(DP), dimension(:), intent(in) :: A
       real(DP)                           :: B
diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90
index 5de45e296..6bb9ece5e 100644
--- a/src/orbel/orbel.f90
+++ b/src/orbel/orbel.f90
@@ -129,7 +129,7 @@ end subroutine orbel_el2xv
 
 
    module pure subroutine orbel_scget(angle, sx, cx)
-      !$omp declare simd(orbel_scget)
+      !!$omp declare simd(orbel_scget)
       !! author: David A. Minton
       !!
       !! Efficiently compute the sine and cosine of an input angle
@@ -182,7 +182,7 @@ end subroutine orbel_scget
    !     REVISIONS:
    !**********************************************************************
    pure subroutine orbel_schget(angle,shx,chx)
-      !$omp declare simd(orbel_schget)
+      !!$omp declare simd(orbel_schget)
       real(DP), intent(in)  ::  angle
       real(DP), intent(out) :: shx,chx
       
@@ -213,7 +213,7 @@ end subroutine orbel_schget
    !     REVISIONS:
    !**********************************************************************
    real(DP) pure function orbel_flon(e,icapn)
-      !$omp declare simd(orbel_flon)
+      !!$omp declare simd(orbel_flon)
       implicit none
       real(DP), intent(in) ::  e, icapn
       integer(I4B) :: iflag,i
@@ -301,7 +301,7 @@ end function orbel_flon
    !     REVISIONS: 2/26/93 hfl
    !**********************************************************************
    real(DP) pure function orbel_fget(e,capn)
-      !$omp declare simd(orbel_fget)
+      !!$omp declare simd(orbel_fget)
       implicit none
 
       real(DP), intent(in) ::  e,capn
@@ -372,7 +372,7 @@ end function orbel_fget
    !       series for small Q.
    !**********************************************************************
    real(DP) pure function orbel_zget(iq)
-      !$omp declare simd(orbel_zget)
+      !!$omp declare simd(orbel_zget)
       implicit none
 
       real(DP), intent(in)  :: iq
@@ -428,7 +428,7 @@ end function orbel_zget
    !     REVISIONS: 2/26/93 hfl
    !**********************************************************************
    real(DP) pure function orbel_esolmd(e,m)
-      !$omp declare simd(orbel_esolmd)
+      !!$omp declare simd(orbel_esolmd)
       implicit none
 
       real(DP), intent(in)  :: e
@@ -484,7 +484,7 @@ end function orbel_esolmd
    !     REVISIONS:
    !**********************************************************************
    real(DP) pure function orbel_ehie(e,im)
-      !$omp declare simd(orbel_ehie)
+      !!$omp declare simd(orbel_ehie)
       implicit none
 
       real(DP), intent(in) :: e,im
@@ -560,7 +560,7 @@ end function orbel_ehie
    !          we have an ellipse with e between 0.15 and 0.8
    !**********************************************************************
    real(DP) pure function orbel_eget(e,m)
-      !$omp declare simd(orbel_eget)
+      !!$omp declare simd(orbel_eget)
       implicit none
       
       real(DP), intent(in) ::  e,m
@@ -634,7 +634,7 @@ end function orbel_eget
    !     REVISIONS: 2/26/93 hfl
    !**********************************************************************
    real(DP) pure function orbel_ehybrid(e,m)
-      !$omp declare simd(orbel_ehybrid)
+      !!$omp declare simd(orbel_ehybrid)
       implicit none
 
       real(DP), intent(in) :: e,m
@@ -674,7 +674,7 @@ end function orbel_ehybrid
    !     REVISIONS: 2/26/93 hfl
    !**********************************************************************
    real(DP) pure function orbel_fhybrid(e,n)
-      !$omp declare simd(orbel_fhybrid)
+      !!$omp declare simd(orbel_fhybrid)
       implicit none
       real(DP), intent(in) :: e,n
 
@@ -694,7 +694,7 @@ end function orbel_fhybrid
    
 
    module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
-      !$omp declare simd(orbel_xv2aeq)
+      !!$omp declare simd(orbel_xv2aeq)
       !! author: David A. Minton
       !!
       !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity
@@ -760,7 +760,7 @@ end subroutine orbel_xv2aeq
 
 
    module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
-      !$omp declare simd(orbel_xv2aqt)
+      !!$omp declare simd(orbel_xv2aqt)
       !! author: David A. Minton
       !!
       !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from
@@ -899,7 +899,7 @@ end subroutine orbel_xv2el_vec
 
 
    pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
-      !$omp declare simd(orbel_xv2el)
+      !!$omp declare simd(orbel_xv2el)
       !! author: David A. Minton
       !!
       !! Compute osculating orbital elements from relative Cartesian position and velocity
diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90
index b637a9df4..fefad67e7 100644
--- a/src/symba/symba_kick.f90
+++ b/src/symba/symba_kick.f90
@@ -113,7 +113,8 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg)
       if (self%nbody == 0) return
       select type(system)
       class is (symba_nbody_system)
-         associate(tp => self, cb => system%cb, pl => system%pl, pltpenc_list => system%pltpenc_list, npltpenc => system%pltpenc_list%nenc)
+         associate(tp => self, cb => system%cb, pl => system%pl, &
+            pltpenc_list => system%pltpenc_list, npltpenc => system%pltpenc_list%nenc)
             call helio_kick_getacch_tp(tp, system, param, t, lbeg)
             ! Remove accelerations from encountering pairs
             do k = 1, npltpenc
diff --git a/src/util/util_flatten.f90 b/src/util/util_flatten.f90
index e3e97c1e7..0dddd95d5 100644
--- a/src/util/util_flatten.f90
+++ b/src/util/util_flatten.f90
@@ -3,7 +3,7 @@
 contains
 
    module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
-      !$omp declare simd(util_flatten_eucl_ij_to_k)
+      !!$omp declare simd(util_flatten_eucl_ij_to_k)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions.
diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90
index 3278b6de4..cac951a99 100644
--- a/src/util/util_minimize_bfgs.f90
+++ b/src/util/util_minimize_bfgs.f90
@@ -38,7 +38,7 @@ module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1)
       real(DP) :: astar                         !! 1D minimized value 
       real(DP), dimension(N) :: y, P
       real(DP), dimension(N,N) :: PP, PyH, HyP
-      real(DP) :: yHy, Py
+      real(DP), save :: yHy, Py
       type(ieee_status_type) :: original_fpe_status
       logical, dimension(:), allocatable :: fpe_flag