diff --git a/tests/bflux/run.sh b/tests/bflux/run.sh
deleted file mode 100755
index 929912d3..00000000
--- a/tests/bflux/run.sh
+++ /dev/null
@@ -1,247 +0,0 @@
-#!/bin/bash 
-# Hyerin (02/17/23) copied from Ben's code
-
-# Bash script testing nonzero b flux
-
-# User specified values here
-KERR=false
-JITTER=true #false #
-ONEZONE=false #true # 
-SHORT_T_OUT=false #true # 
-ROT=false #true #
-ROT_UPHI=0.1
-bz=1e-4 #2e-8 #1e-8 #0 #
-DIM=3 #2 # 
-NZONES=3 #4 # 8 #
-BASE=8
-NRUNS=3000
-START_RUN=494 # if !=0, change start_time, out_to_in, iteration, r_out, r_in
-res=128 #64 #96 #32 #
-#DRTAG="bondi_multizone_050423_bflux0_${bz}_${res}^3_n8_test_faster_rst"
-DRTAG="bondi_multizone_050523_bflux0_${bz}_${res}^3_n3_noshort"
-#DRTAG="bondi_multizone_050423_onezone_bflux0_${bz}_2d_n4"
-
-# Set paths
-KHARMADIR=../..
-PDR="/n/holylfs05/LABS/bhi/Users/hyerincho/grmhd/" ## parent directory
-DR="${PDR}data/${DRTAG}"
-parfilename="${PDR}/kharma/pars/bondi_multizone/bondi_multizone_00000.par" # parameter file
-
-# other values determined automatically
-turn_around=$(($NZONES-1))
-start_time=1424043 #0 #`echo "scale=20; 85184.643106985808" | bc -l` # #
-out_to_in=1 #-1 #
-iteration=247 #
-if [[ $ONEZONE == "true" ]]; then
-  r_out=$((${BASE}**($turn_around+2))) # 4096 
-  r_in=1 
-  NRUNS=4 #1
-  #NZONES=1
-else
-  r_out=64 #$((${BASE}**($turn_around+2))) # 
-  r_in=1 #$((${BASE}**$turn_around)) #
-fi
-#nlim=0
-
-# if the directories are not present, make them.
-if [ ! -d "${DR}" ]; then
-  mkdir "${DR}"
-fi
-if [ ! -d "${PDR}logs/${DRTAG}" ]; then
-  mkdir "${PDR}logs/${DRTAG}"
-fi
-
-outermost_zone=$((2*($NZONES-1)))
-
-### Start running zone by zone
-for (( VAR=$START_RUN; VAR<$NRUNS; VAR++ ))
-do
-  args=()
-  echo "${DRTAG}: iter $iteration, $VAR : t = $start_time, r_out = $r_out, r_in = $r_in"
-  if [[ $NZONES -eq 3 ]]; then
-    r_b=256
-  elif [[ $NZONES -eq 4 ]]; then
-    r_b=256
-  else
-    r_b=100000
-  fi
-  logruntime=`echo "scale=20; l($r_out)*3./2-l(1.+$r_out/$r_b)/2." | bc -l` # round to an integer for the free-fall time (cs^2=0.01 should be updated from the desired rs value) # GIZMO
-  runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-  runtime=$(($runtime/2)) # test
-  if [[ $SHORT_T_OUT == "true" ]]; then
-    # use the next largest annulus' runtime
-    if [ $(($VAR % $outermost_zone)) -eq 0 ]; then
-      r_out_sm=$((${r_out}/${BASE}))
-      echo "r_out=${r_out}, but I will use next largest annulus' r_out=${r_out_sm} for the runtime"
-      logruntime=`echo "scale=20; l($r_out_sm)*3./2-l(1.+$r_out_sm/$r_b)/2." | bc -l`
-      runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-      runtime=$(($runtime/2)) # test
-    fi
-  fi
-
-  #runtime=$(($runtime/4)) # test
-  logrho=-8.2014518
-  log_u_over_rho=-5.2915149 # test same vacuum conditions as r_shell when (rs=1e2.5)
-  #if [ $VAR -ne 0 ]; then
-  #  tag=($( tail -n 5 ${PDR}/logs/${DRTAG}/log_multizone$(printf %05d $((${VAR}-1)))_out ))
-  #  echo ${tag[0]:5}
-  #  start_time=$(printf "%.18g" "${tag[0]:5}") # last cycle number
-  #  start_time=`echo "scale=0; start_time" | bc -l`
-  #  echo "start time $start_time"
-  #fi
-  start_time=$(($start_time+$runtime))  
-
-  #parfilename="../../kharma/pars/bondi_multizone/bondi_multizone_$(printf %05d ${VAR}).par" # parameter file
-  
-  # set problem type and cleanup
-  if [ $VAR -eq 0 ]; then
-    prob="bondi"
-  else
-    prob="resize_restart_kharma"
-  fi
-  
-  # set BH spin
-  if [[ $KERR == "true" ]]; then
-    spin=0.99
-  else
-    spin=0.0
-  fi
-
-  # set flow rotation
-  if [[ $ROT == "true" ]]; then
-    args+=(" bondi/uphi=$ROT_UPHI ")
-  else
-    args+=(" bondi/uphi=0 ")
-  fi
-  
-  # output time steps
-  #output0_dt=$((${runtime}/100*10))
-  output0_dt=$((${runtime}/200*10))
-  if (( $(echo "$output0_dt < 0.00001" |bc -l) )); then
-    output0_dt=1
-  fi
-  #output1_dt=$((${runtime}/100*10))
-  output1_dt=$((${runtime}/50*10))
-  output2_dt=$((${runtime}/1000*10))
-  
-  # dt, fname, fname_fill
-  if [ $VAR -ne 0 ]; then
-    # update dt from the previous run
-    tag=($( tail -n 10 ${PDR}/logs/${DRTAG}/log_multizone$(printf %05d $((${VAR}-1)))_out ))
-    dt=$(printf "%.18g" "${tag[2]:3}") # previous dt
-    #tag=($( tail -n 5 ${PDR}/logs/${DRTAG}/log_multizone$(printf %05d $((${VAR}-1)))_out ))
-    #nlim=$(printf "%d" "${tag[1]:6}") # last cycle number
-
-    dt_new=$(echo "scale=14; $dt*sqrt($BASE^(-3*$out_to_in))/4" | bc -l) # new dt ## TODO: r^3/2
-    if (( $(echo "$dt_new > 0.00001" |bc -l) )); then
-      dt_new=$dt_new
-    else
-      dt_new=0.00001
-    fi
-    fname_dir="${DR}/bondi_multizone_$(printf %05d $((${VAR}-1)))"
-    fname=$(find ${fname_dir} -type f -iname "*final.rhdf")
-    if [ $VAR -ge $NZONES ]; then
-      fname_fill_num=$((2*($iteration-1)*(${NZONES}-1)-${VAR}))
-      fname_fill_dir="${DR}/bondi_multizone_$(printf %05d $fname_fill_num)"
-      fname_fill=$(find ${fname_fill_dir} -type f -iname "*final.rhdf")
-    else
-      fname_fill="none"
-    fi
-    args+=(" resize_restart/fname=$fname parthenon/time/dt_min=$dt_new")
-    args+=(" resize_restart/fname_fill=$fname_fill ")
-  else
-    r_shell=$((${r_out}/2))
-    args+=(" bondi/r_shell=$r_shell ")
-    
-    if [[ $JITTER == "true" ]]; then
-        args+=(" perturbation/u_jitter=0.1 ")
-    else
-        args+=(" perturbation/u_jitter=0.0 ")
-    fi
-  fi
-
-  # data_dir, logfiles
-  data_dir="${DR}/bondi_multizone_$(printf %05d ${VAR})"
-  out_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_out"
-  err_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_err"
-  
-  # configuration
-  if [[ $DIM -gt 2 ]]; then
-    nx3_m=$res #32 #
-    nx3_mb=$(($res/2)) #$res #
-    nx2_mb=$(($res/2))
-  else
-    nx3_m=1
-    nx3_mb=1
-    nx2_mb=$res  #
-  fi
-  if [[ $NZONES -eq 3 ]]; then
-    logrho=-4.13354231 #-5.686638255139154
-    log_u_over_rho=-2.57960521 #-3.6150030239527497
-    args+=(" bondi/rs=16 ")
-  fi
-  if [[ $NZONES -eq 4 ]]; then
-    logrho=-4.200592800419657
-    log_u_over_rho=-2.62430556
-    args+=(" bondi/rs=16 ")
-  fi
-  if [[ $ONEZONE == "true" ]]; then
-      #args+=(" bondi/r_shell=256 ")
-      nx1_m=$((32*($NZONES+1)))
-      nx1_mb=$nx1_m #128
-      start_time=5100000 #10000000
-      output1_dt=1000 #660 #6420 #2310
-      nlim=-1
-  else
-      nx1_m=$res #64
-      nx1_mb=$res #64
-      nlim=$((100000*$res*$res/64/64)) #$(($nlim+50000)) ##$((50000*($VAR+1)))
-  fi
-  if [[ $bz == 0 ]]; then
-    nlim=-1
-  fi
-
-  #srun --mpi=pmix ${PDR}/kharma_nlim.cuda -i ${parfilename} \
-  srun --mpi=pmix ${PDR}/kharma_faster_rst_fixed.cuda -i ${parfilename} \
-                                    parthenon/mesh/nx1=$nx1_m parthenon/mesh/nx2=$res parthenon/mesh/nx3=$nx3_m \
-                                    parthenon/meshblock/nx1=$nx1_mb parthenon/meshblock/nx2=$nx2_mb parthenon/meshblock/nx3=$nx3_mb \
-                                    parthenon/job/problem_id=$prob \
-                                    parthenon/time/tlim=${start_time} \
-                                    parthenon/time/nlim=$nlim \
-                                    coordinates/r_in=${r_in} coordinates/r_out=${r_out} coordinates/a=$spin coordinates/ext_g=false\
-                                    coordinates/transform=mks coordinates/hslope=0.3 \
-                                    bondi/vacuum_logrho=${logrho} bondi/vacuum_log_u_over_rho=${log_u_over_rho} \
-                                    floors/disable_floors=false floors/rho_min_geom=1e-6 floors/u_min_geom=1e-8 floors/bsq_over_u_max=1e+20 floors/adjust_k=0 \
-                                    floors/bsq_over_rho_max=100 floors/u_over_rho_max=100 floors/gamma_max=10 floors/frame=drift \
-                                    GRMHD/reconstruction=linear_vl GRMHD/cfl=0.5 GRMHD/add_jcon=0 \
-                                    b_field/type=vertical b_field/solver=flux_ct b_field/bz=${bz} \
-                                    b_field/fix_flux_x1=1 b_field/initial_cleanup=0 \
-                                    resize_restart/base=$BASE resize_restart/nzone=$NZONES  resize_restart/iteration=$iteration\
-                                    parthenon/output0/dt=$output0_dt \
-                                    parthenon/output1/dt=$output1_dt \
-                                    parthenon/output2/dt=$output2_dt \
-                                    debug/flag_verbose=2 debug/extra_checks=1 \
-                                    ${args[@]} \
-                                    -d ${data_dir} 1> ${out_fn} 2>${err_fn}
-                                    #parthenon/time/tlim=10000000 \
-                                    #parthenon/output1/dt=2310 \
-                                    #  b_field/type=vertical b_field/bz=${bz}flux_ct
-                                    #floors/bsq_over_rho_max=100 floors/u_over_rho_max=2 \  
-
-  if [ $VAR -ne 0 ]; then
-    if [ $(($VAR % ($NZONES-1))) -eq 0 ]; then
-      out_to_in=$(($out_to_in*(-1)))
-      iteration=$(($iteration+1))
-    fi
-  fi
-
-  if [ $out_to_in -gt 0 ]; then
-    # half the radii
-    r_out=$((${r_out}/$BASE))
-    r_in=$((${r_in}/$BASE))
-  else
-    # double the radii
-    r_out=$((${r_out}*$BASE))
-    r_in=$((${r_in}*$BASE))
-  fi
-done
diff --git a/tests/bondi_multizone/bondi_multizone.par b/tests/bondi_multizone/bondi_multizone.par
deleted file mode 100755
index ee6c2e45..00000000
--- a/tests/bondi_multizone/bondi_multizone.par
+++ /dev/null
@@ -1,103 +0,0 @@
-# Multizone Bondi flow template
-# Many options replaced by run script
-
-<parthenon/job>
-problem_id = bondi
-
-<parthenon/mesh>
-refinement = none
-numlevel = 1
-nx1 = 128
-nx2 = 128
-nx3 = 128
-
-<parthenon/meshblock>
-nx1 = 128
-nx2 = 64
-nx3 = 64
-
-<coordinates>
-base = ks
-transform = exponential
-a = 0.0 # spin updated from run_kharma.sh
-
-<parthenon/time>
-tlim = 5289680481 # updated from run_kharma.sh
-nlim = -1 # updated from run_kharma.sh
-dt_min = 0.00001
-
-<GRMHD>
-cfl = 0.9
-gamma = 1.666667
-reconstruction = linear_vl
-
-<bondi>
-mdot = 1.0
-rs = 316.22776601683796  # (1e2.5)#1000.0 #8.0 #300.0 #
-vacuum_logrho= -8.2014518
-vacuum_log_u_over_rho = -5.2915149 # updated from run_kharma.sh
-r_shell = 8388608 # updated from run_kharma.sh
-use_gizmo = false
-
-<gizmo_shell>
-datfn = none
-
-<resize_restart>
-fname = none # updated from run_kharma.sh
-fname_fill = none #, updated from run_kharma.sh
-use_dt = false
-base = 8
-nzone = 7
-
-<floors>
-disable_floors = true
-rho_min_geom = 1.0e-12
-u_min_geom = 1.0e-15
-
-<bounds>
-check_inflow_outer = false
-check_inflow_inner = false # Hyerin test (12/22/22)
-fix_flux_pole      = true
-
-<perturbation>
-u_jitter=0
-
-<b_field>
-type = vertical # b_field_type updated from run_kharma.sh
-solver = none # b_field_solver updated from run_kharma.sh
-#norm = true # Hyerin (12/29/22) this increases divB in the boundaries, so won't use here
-#beta_min=1000 # Hyerin (12/29/22)
-bz = 1e-4
-#fix_polar_flux = 0 # Hyerin test (12/22/22)
-fix_flux_x1      = 0 #1 # Hyerin test (02/16/23)
-initial_cleanup = true # updated from run_kharma.sh
-
-<b_cleanup>
-rel_tolerance=1.e-8
-
-<debug>
-verbose = 1
-
-<driver>
-type = harm
-two_sync = true
-
-<implicit>
-max_nonlinear_iter = 3
-
-<parthenon/output0>
-file_type = hdf5
-dt = 528968040 # output0_dt updated from run_kharma.sh
-single_precision_output = true
-variables = prims.rho, prims.u, prims.uvec, prims.B, fflag, pflag
-ghost_zones = true
-
-<parthenon/output1>
-file_type = rst
-dt = 2644840240 # output1_dt updated from run_kharma.sh
-ghost_zones = true
-
-<parthenon/output2>
-file_type = hst
-dt = 52896800 # output2_dt updated from run_kharma.sh
-
diff --git a/tests/bondi_multizone/kharma_multizone.sb b/tests/bondi_multizone/kharma_multizone.sb
deleted file mode 100755
index 9b9f3655..00000000
--- a/tests/bondi_multizone/kharma_multizone.sb
+++ /dev/null
@@ -1,40 +0,0 @@
-#!/bin/bash
-# Admin stuff
-#SBATCH -J KHARMA
-#SBATCH -t 2:00:00
-#SBATCH -N 1
-#SBATCH -o "out-%j.txt"
-
-# Partition
-#SBATCH -p gpu_debug --reservation gpu_debug --qos debug
-##SBATCH -p gpu
-
-# Nodes we want
-#SBATCH --tasks-per-node=4
-# TODO 80GB specific?
-
-# Node options
-# Note we could do 32 if HT is faster
-#SBATCH --cpus-per-task=16
-# ALWAYS reserve full nodes to mitigate memory leaks
-#SBATCH --exclusive
-#SBATCH --mem=0
-
-# LANL Chicoma run script
-
-# KHARMA directory
-KHARMA_DIR="$HOME/Code/kharma-multizone"
-
-# OpenMP directives: use all available threads
-# TODO Chicoma doesn't allocate correctly
-#export OMP_PROC_BIND=spread
-#export OMP_PLACES=threads
-
-# Load any defaults/modules from the machine file
-HOST=ch-fe1.lanl.gov
-ARGS=$(cat $KHARMA_DIR/make_args)
-source $KHARMA_DIR/machines/lanl.sh
-
-# Run with srun
-# TODO auto-switch to mpirun in interactive?
-./run.sh -t 11:55:00 "$@"
diff --git a/tests/bondi_multizone/run.sh b/tests/bondi_multizone/run.sh
deleted file mode 100755
index 0aef4b18..00000000
--- a/tests/bondi_multizone/run.sh
+++ /dev/null
@@ -1,195 +0,0 @@
-#!/bin/bash 
-# Hyerin (02/17/23) copied from Ben's code
-
-# Bash script testing HD bondi
-
-# User specified values here
-KERR=false
-JITTER=false #true #
-ROT=false #true #
-SUPEREXP=false #true #
-SHORT_T_OUT=true
-ROT_UPHI=0.1 # 1e-10
-NZONES=8
-BASE=8
-NRUNS=300
-
-# Set paths
-KHARMA_DIR=$HOME/Code/kharma-multizone
-PDR=$(pwd)/..
-DR=$(pwd)
-parfilename="bondi_multizone.par" # parameter file
-
-turn_around=$(($NZONES-1))
-
-# Read some values if we're restarting
-if [ -f start_run ]; then
-  START_RUN=$(cat start_run.var)
-  start_time=$(cat start_time.var)
-  out_to_in=$(cat out_to_in.var)
-  iteration=$(cat iteration.var)
-  r_in=$(cat r_in.var)
-  r_out=$(cat r_out.var)
-
-  last_run=$(($START_RUN - 1))
-  data_dir="${DR}/$(printf %05d $last_run)"
-  out_fn="${DR}/$(printf %05d $last_run).log"
-else
-  START_RUN=0
-  start_time=0
-  out_to_in=1
-  iteration=1 # eq : (iteration-1)*(NZONES-1)<VAR<=iteration*(NZONES-1)
-  r_in=$((${BASE}**$turn_around))
-  r_out=$((${BASE}**($turn_around+2)))
-  echo "STARTING MULTIZONE RUN: $START_RUN $start_time $out_to_in $iteration $r_in $r_out"
-fi
-
-# other values determined automatically
-ZONE0=$((2*($NZONES-1)))
-r_b=100000
-log_u_over_rho=-5.2915149 # test same vacuum conditions as r_shell when (rs=1e2.5)
-outermost_zone=$((2*($NZONES-1)))
-
-# if the directories are not present, make them.
-mkdir -p "${DR}"
-
-outermost_zone=$((2*($NZONES-1)))
-r_b=100000
-
-### Start running zone by zone
-for (( VAR=$START_RUN; VAR<$NRUNS; VAR++ ))
-do
-  # Save variables as they are when entering the loop, for restarting
-  echo $VAR > start_run.var
-  echo $start_time > start_time.var
-  echo $out_to_in > out_to_in.var
-  echo $iteration > iteration.var
-  echo $r_in > r_in.var
-  echo $r_out > r_out.var
-
-  args=()
-  echo "iter $iteration, $VAR : t = $start_time, r_out = $r_out, r_in = $r_in"
-  logruntime=`echo "scale=20; l($r_out)*3./2-l(1.+$r_out/$r_b)/2." | bc -l` # round to an integer for the free-fall time (cs^2=0.01 should be updated from the desired rs value) # GIZMO
-  runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-  if [[ $SHORT_T_OUT == "true" ]]; then
-    # use the next largest annulus' runtime
-    if [ $(($VAR % $outermost_zone)) -eq 0 ]; then
-      r_out_sm=$((${r_out}/${BASE}))
-      echo "r_out=${r_out}, but I will use next largest annulus' r_out=${r_out_sm} for the runtime"
-      logruntime=`echo "scale=20; l($r_out_sm)*3./2-l(1.+$r_out_sm/$r_b)/2." | bc -l`
-      runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-    fi
-  fi
-  start_time=$(($start_time+$runtime))  
-
-  # set problem type and cleanup
-  if [ $VAR -eq 0 ]; then
-    prob="bondi"
-  else
-    prob="resize_restart_kharma"
-  fi
-  
-  # set BH spin
-  if [[ $KERR == "true" ]]; then
-    spin=0.99
-  else
-    spin=0.0
-  fi
-  
-  # set flow rotation
-  if [[ $ROT == "true" ]]; then
-    args+=(" bondi/uphi=$ROT_UPHI ")
-  else
-    args+=(" bondi/uphi=0 ")
-  fi
-
-  
-  # output time steps
-  output0_dt=$((${runtime}/100*10))
-  output1_dt=$((${runtime}/20*10))
-  output2_dt=$((${runtime}/1000*10))
-  
-  # dt, fname, fname_fill
-  if [ $VAR -ne 0 ]; then
-    # update dt from the previous run
-    tag=($( tail -n 10 $out_fn))
-    dt=$(printf "%.18g" "${tag[2]:3}") # previous dt
-    dt_new=$(echo "scale=14; $dt*sqrt($BASE^(-3*$out_to_in))/4" | bc -l) # new dt ## TODO: r^3/2
-    if (( $(echo "$dt_new > 0.00001" |bc -l) )); then
-      dt_new=$dt_new
-    else
-      dt_new=0.00001
-    fi
-    fname=$(find $data_dir -type f -iname "*final.rhdf")
-    if [ $VAR -ge $NZONES ]; then
-      fname_fill_num=$((2*($iteration-1)*(${NZONES}-1)-${VAR}))
-      fname_fill_dir="${DR}/$(printf %05d $fname_fill_num)"
-      fname_fill=$(find ${fname_fill_dir} -type f -iname "*final.rhdf")
-    else
-      fname_fill="none"
-    fi
-    args+=(" resize_restart/fname=$fname resize_restart/use_dt=false parthenon/time/dt_min=$dt_new")
-    args+=(" resize_restart/fname_fill=$fname_fill ")
-  else
-    r_shell=$((${r_out}/2))
-    args+=(" bondi/r_shell=$r_shell ")
-    if [[ $JITTER == "true" ]]; then
-        args+=(" perturbation/u_jitter=0.1 ")
-    else
-        args+=(" perturbation/u_jitter=0.0 ")
-    fi
-  fi
-
-  if [[ $SUPEREXP == "true" ]]; then
-    r_br=4e6 # such that it does not lie exactly on the zone centers or faces #$((2*${BASE}**($NZONES))) # break radius (double, in order to be safe from boundaries)
-    r_superexp=1e8 # outermost radius of super exponential grid
-
-    if [ $(($VAR % $ZONE0)) -eq 0 ]; then
-      args+=(" coordinates/transform=superexp ")
-      args+=(" coordinates/r_br=${r_br} coordinates/npow=2.0 coordinates/cpow=1.0 ")
-      #args+=(" coordinates/r_out=${r_superexp} ")
-      #args+=(" coordinates/r_out=${r_out} ")
-    else
-      args+=(" coordinates/transform=eks ")
-
-    fi
-    args+=(" coordinates/r_out=${r_out} ")
-  else
-    args+=(" coordinates/transform=mks coordinates/hslope=1 coordinates/r_out=${r_out} ")
-  fi
-  
-
-  # data_dir, logfiles
-  data_dir="${DR}/$(printf %05d ${VAR})"
-  out_fn="${DR}/$(printf %05d ${VAR}).log"
-  
-  ${KHARMA_DIR}/run.sh -i ${parfilename} \
-                       parthenon/job/problem_id=$prob \
-                       parthenon/time/tlim=${start_time} \
-                       coordinates/r_in=${r_in} coordinates/a=$spin coordinates/ext_g=false \
-                       bondi/vacuum_log_u_over_rho=${log_u_over_rho} bondi/use_gizmo=false \
-                       b_field/type=none b_field/solver=none b_field/bz=1e-3 \
-                       b_field/fix_flux_x1=false b_field/initial_cleanup=false \
-                       resize_restart/base=$BASE resize_restart/nzone=$NZONES resize_restart/iteration=$iteration\
-                       parthenon/output0/dt=$output0_dt \
-                       parthenon/output1/dt=$output1_dt \
-                       parthenon/output2/dt=$output2_dt \
-                       ${args[@]} "$@" \
-                       -d ${data_dir} 1>${out_fn} 2>&1
-  if [ $VAR -ne 0 ]; then
-    if [ $(($VAR % ($NZONES-1))) -eq 0 ]; then
-      out_to_in=$(($out_to_in*(-1)))
-      iteration=$(($iteration+1))
-    fi
-  fi
-
-  if [ $out_to_in -gt 0 ]; then
-    # half the radii
-    r_out=$((${r_out}/$BASE))
-    r_in=$((${r_in}/$BASE))
-  else
-    # double the radii
-    r_out=$((${r_out}*$BASE))
-    r_in=$((${r_in}*$BASE))
-  fi
-done
diff --git a/tests/bondi_multizone/run_prev.sh b/tests/bondi_multizone/run_prev.sh
deleted file mode 100755
index 581eb7e4..00000000
--- a/tests/bondi_multizone/run_prev.sh
+++ /dev/null
@@ -1,189 +0,0 @@
-#!/bin/bash 
-# Hyerin (02/17/23) copied from Ben's code
-
-# Bash script testing HD bondi
-
-# User specified values here
-KERR=false
-JITTER=false #true #
-ROT=false #true #
-SUPEREXP=false #true #
-SHORT_T_OUT=true #false
-ROT_UPHI=0.1 # 1e-10
-DIM=3 #2 #
-NZONES=8 #7
-BASE=8
-NRUNS=300
-START_RUN=0 # if this is not 0, then update start_time, out_to_in, iteration, r_out, r_in to values that you are re-starting from
-#DRTAG="bondi_multizone_041723_bondi_64^3_rot${ROT_UPHI}keprst_jit"
-DRTAG="bondi_multizone_052423_bondi_64^3_n8_noshock"
-
-# Set paths
-KHARMA_DIR=../..
-PDR="/n/holylfs05/LABS/bhi/Users/hyerincho/grmhd/" ## parent directory
-DR="${PDR}data/${DRTAG}"
-parfilename="${PDR}/kharma/pars/bondi_multizone/bondi_multizone_00000.par" # parameter file
-
-# other values determined automatically
-turn_around=$(($NZONES-1))
-start_time=0 #100337021259 #
-out_to_in=1 # -1 #
-iteration=1 # eq : (iteration-1)*(NZONES-1)<VAR<=iteration*(NZONES-1)
-r_out=$((${BASE}**($turn_around+2))) #512 #
-r_in=$((${BASE}**$turn_around)) #8 #
-
-# if the directories are not present, make them.
-mkdir -p "${DR}"
-mkdir -p "${PDR}logs/${DRTAG}"
-
-outermost_zone=$((2*($NZONES-1)))
-r_b=100000
-
-### Start running zone by zone
-for (( VAR=$START_RUN; VAR<$NRUNS; VAR++ ))
-do
-
-  args=()
-  echo "${DRTAG} iter $iteration, $VAR : t = $start_time, r_out = $r_out, r_in = $r_in"
-  logruntime=`echo "scale=20; l($r_out)*3./2-l(1.+$r_out/$r_b)/2." | bc -l` # round to an integer for the free-fall time (cs^2=0.01 should be updated from the desired rs value) # GIZMO
-  runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-  if [[ $SHORT_T_OUT == "true" ]]; then
-    # use the next largest annulus' runtime
-    if [ $(($VAR % $outermost_zone)) -eq 0 ]; then
-      r_out_sm=$((${r_out}/${BASE}))
-      echo "r_out=${r_out}, but I will use next largest annulus' r_out=${r_out_sm} for the runtime"
-      logruntime=`echo "scale=20; l($r_out_sm)*3./2-l(1.+$r_out_sm/$r_b)/2." | bc -l`
-      runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-    fi
-  fi
-  log_u_over_rho=0 # test case when less shock #-5.2915149 # test same vacuum conditions as r_shell when (rs=1e2.5)
-  start_time=$(($start_time+$runtime))  
-
-  # set problem type and cleanup
-  if [ $VAR -eq 0 ]; then
-    prob="bondi"
-  else
-    prob="resize_restart_kharma"
-  fi
-  
-  # set BH spin
-  if [[ $KERR == "true" ]]; then
-    spin=0.99
-  else
-    spin=0.0
-  fi
-  
-  # set flow rotation
-  if [[ $ROT == "true" ]]; then
-    args+=(" bondi/uphi=$ROT_UPHI ")
-  else
-    args+=(" bondi/uphi=0 ")
-  fi
-
-  
-  # output time steps
-  output0_dt=$((${runtime}/100*10))
-  output1_dt=$((${runtime}/20*10))
-  output2_dt=$((${runtime}/1000*10))
-  
-  # dt, fname, fname_fill
-  if [ $VAR -ne 0 ]; then
-    # update dt from the previous run
-    tag=($( tail -n 10 $out_fn))
-    dt=$(printf "%.18g" "${tag[2]:3}") # previous dt
-    dt_new=$(echo "scale=14; $dt*sqrt($BASE^(-3*$out_to_in))/4" | bc -l) # new dt ## TODO: r^3/2
-    if (( $(echo "$dt_new > 0.00001" |bc -l) )); then
-      dt_new=$dt_new
-    else
-      dt_new=0.00001
-    fi
-    fname=$(find $data_dir -type f -iname "*final.rhdf")
-    if [ $VAR -ge $NZONES ]; then
-      fname_fill_num=$((2*($iteration-1)*(${NZONES}-1)-${VAR}))
-      fname_fill_dir="${DR}/bondi_multizone_$(printf %05d $fname_fill_num)"
-      fname_fill=$(find ${fname_fill_dir} -type f -iname "*final.rhdf")
-    else
-      fname_fill="none"
-    fi
-    args+=(" resize_restart/fname=$fname resize_restart/use_dt=false parthenon/time/dt_min=$dt_new")
-    args+=(" resize_restart/fname_fill=$fname_fill ")
-  else
-    r_shell=$((${r_out}/2))
-    args+=(" bondi/r_shell=$r_shell ")
-    if [[ $JITTER == "true" ]]; then
-        args+=(" perturbation/u_jitter=0.1 ")
-    else
-        args+=(" perturbation/u_jitter=0.0 ")
-    fi
-  fi
-
-  if [[ $SUPEREXP == "true" ]]; then
-    r_br=4e6 # such that it does not lie exactly on the zone centers or faces #$((2*${BASE}**($NZONES))) # break radius (double, in order to be safe from boundaries)
-    r_superexp=1e8 # outermost radius of super exponential grid
-
-    if [ $(($VAR % $ZONE0)) -eq 0 ]; then
-      args+=(" coordinates/transform=superexp ")
-      args+=(" coordinates/r_br=${r_br} coordinates/npow=2.0 coordinates/cpow=1.0 ")
-      #args+=(" coordinates/r_out=${r_superexp} ")
-      #args+=(" coordinates/r_out=${r_out} ")
-    else
-      args+=(" coordinates/transform=eks ")
-
-    fi
-    args+=(" coordinates/r_out=${r_out} ")
-  else
-    args+=(" coordinates/transform=mks coordinates/hslope=1 coordinates/r_out=${r_out} ")
-    #args+=(" coordinates/transform=eks coordinates/r_out=${r_out} ")
-  fi
-  
-
-  # data_dir, logfiles
-  data_dir="${DR}/bondi_multizone_$(printf %05d ${VAR})"
-  out_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_out"
-  err_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_err"
-  
-  # configuration
-  if [[ $DIM -gt 2 ]]; then
-    nx3_m=64 #
-    nx3_mb=64 #
-  else
-    nx3_m=1
-    nx3_mb=1
-  fi
-
-  #srun --mpi=pmix ${PDR}/kharma_rocky_2.cuda -i ${parfilename} \
-  srun --mpi=pmix ${PDR}/kharma_faster_rst_fixed.cuda -i ${parfilename} \
-                                    parthenon/job/problem_id=$prob \
-                                    parthenon/time/tlim=${start_time} parthenon/time/nlim=-1 \
-                                    parthenon/mesh/nx1=64 parthenon/mesh/nx2=64 parthenon/mesh/nx3=$nx3_m \
-                                    parthenon/meshblock/nx1=32 parthenon/meshblock/nx2=32 parthenon/meshblock/nx3=$nx3_mb \
-                                    coordinates/r_in=${r_in} coordinates/a=$spin coordinates/ext_g=false \
-                                    GRMHD/reconstruction=linear_vl \
-                                    bounds/fix_flux_pole=1 \
-                                    bondi/vacuum_logrho=-8.2014518 bondi/vacuum_log_u_over_rho=${log_u_over_rho} bondi/use_gizmo=false \
-                                    b_field/type=none b_field/solver=none b_field/bz=1e-3 \
-                                    b_field/fix_flux_x1=0 b_field/initial_cleanup=0 \
-                                    resize_restart/base=$BASE resize_restart/nzone=$NZONES resize_restart/iteration=$iteration\
-                                    parthenon/output0/dt=$output0_dt \
-                                    parthenon/output1/dt=$output1_dt \
-                                    parthenon/output2/dt=$output2_dt \
-                                    ${args[@]} \
-                                    -d ${data_dir} 1> ${out_fn} 2>${err_fn}
-
-  if [ $VAR -ne 0 ]; then
-    if [ $(($VAR % ($NZONES-1))) -eq 0 ]; then
-      out_to_in=$(($out_to_in*(-1)))
-      iteration=$(($iteration+1))
-    fi
-  fi
-
-  if [ $out_to_in -gt 0 ]; then
-    # half the radii
-    r_out=$((${r_out}/$BASE))
-    r_in=$((${r_in}/$BASE))
-  else
-    # double the radii
-    r_out=$((${r_out}*$BASE))
-    r_in=$((${r_in}*$BASE))
-  fi
-done
diff --git a/tests/gizmo_shell/run.sh b/tests/gizmo_shell/run.sh
deleted file mode 100755
index 592c98e9..00000000
--- a/tests/gizmo_shell/run.sh
+++ /dev/null
@@ -1,155 +0,0 @@
-#!/bin/bash 
-# Hyerin (02/18/23) copied from Ben's code
-
-# Bash script testing gizmo shell run (no b field)
-
-# User specified values here
-KERR=false
-EXT_G=false #true #
-SHORT_T_OUT=true #false
-DIM=3
-NZONES=8 #7
-BASE=8
-NRUNS=300
-START_RUN=0 # if this is not 0, then update start_time, out_to_in, iteration, r_out, r_in to values that you are re-starting from
-DRTAG="bondi_multizone_052523_gizmo_n8_64^3_noshock"
-
-# Set paths
-KHARMADIR=../..
-PDR="/n/holylfs05/LABS/bhi/Users/hyerincho/grmhd/" ## parent directory
-DR="${PDR}data/${DRTAG}"
-parfilename="${PDR}/kharma/pars/bondi_multizone/bondi_multizone_00000.par" # parameter file
-
-# other values determined automatically
-turn_around=$(($NZONES-1))
-start_time=0 #19493738388 #
-out_to_in=1 #-1 # 
-iteration=1 #6 # eq : (iteration-1)*(NZONES-1)<VAR<=iteration*(NZONES-1)
-r_out=$((${BASE}**($turn_around+2))) #512 #
-r_in=$((${BASE}**$turn_around)) #8 #
-
-# if the directories are not present, make them.
-if [ ! -d "${DR}" ]; then
-  mkdir "${DR}"
-fi
-if [ ! -d "${PDR}logs/${DRTAG}" ]; then
-  mkdir "${PDR}logs/${DRTAG}"
-fi
-
-outermost_zone=$((2*($NZONES-1)))
-r_b=100000
-
-# GIZMO data
-#datfn="${PDR}/data/gizmo/042623_3d/new_snap_deposit_010.hdf5"
-datfn="${PDR}/data/gizmo/031623_100Myr/dat.txt"
-logrho=-7.80243572
-log_u_over_rho=0 #-5.34068635
-
-### Start running zone by zone
-for (( VAR=$START_RUN; VAR<$NRUNS; VAR++ ))
-do
-  args=()
-  echo "${DRTAG} iter $iteration, $VAR : t = $start_time, r_out = $r_out, r_in = $r_in"
-  logruntime=`echo "scale=20; l($r_out)*3./2-l(1.+$r_out/$r_b)/2." | bc -l` # round to an integer for the free-fall time (cs^2=0.01 should be updated from the desired rs value) # GIZMO
-  if [[ $SHORT_T_OUT == "true" ]]; then
-    # use the next largest annulus' runtime
-    if [ $(($VAR % $outermost_zone)) -eq 0 ]; then
-      r_out_sm=$((${r_out}/${BASE}))
-      echo "r_out=${r_out}, but I will use next largest annulus' r_out=${r_out_sm} for the runtime"
-      logruntime=`echo "scale=20; l($r_out_sm)*3./2-l(1.+$r_out_sm/$r_b)/2." | bc -l`
-    fi
-  fi
-  runtime=`echo "scale=0; e($logruntime)+1" | bc -l`
-  #log_u_over_rho=-5.2915149 # test same vacuum conditions as r_shell when (rs=1e2.5)
-  start_time=$(($start_time+$runtime))  
-  
-  # set problem type and cleanup
-  if [ $VAR -eq 0 ]; then
-    prob="bondi" #"gizmo_shell"
-  else
-    prob="resize_restart_kharma"
-  fi
-  
-  # set BH spin
-  if [[ $KERR == "true" ]]; then
-    spin=0.99
-  else
-    spin=0.0
-  fi
-  
-  # output time steps
-  output0_dt=$((${runtime}/100*10))
-  output1_dt=$((${runtime}/20*10))
-  output2_dt=$((${runtime}/1000*10))
-  
-  # dt, fname, fname_fill
-  if [ $VAR -ne 0 ]; then
-    # update dt from the previous run
-    tag=($( tail -n 10 ${PDR}/logs/${DRTAG}/log_multizone$(printf %05d $((${VAR}-1)))_out ))
-    dt=$(printf "%.18g" "${tag[2]:3}") # previous dt
-    dt_new=$(echo "scale=14; $dt*sqrt($BASE^(-3*$out_to_in))/4" | bc -l) # new dt ## TODO: r^3/2
-    if (( $(echo "$dt_new > 0.00001" |bc -l) )); then
-      dt_new=$dt_new
-    else
-      dt_new=0.00001
-    fi
-    fname_dir="${DR}/bondi_multizone_$(printf %05d $((${VAR}-1)))"
-    fname=$(find ${fname_dir} -type f -iname "*final.rhdf")
-    if [ $VAR -ge $NZONES ]; then
-      fname_fill_num=$((2*($iteration-1)*(${NZONES}-1)-${VAR}))
-      fname_fill_dir="${DR}/bondi_multizone_$(printf %05d $fname_fill_num)"
-      fname_fill=$(find ${fname_fill_dir} -type f -iname "*final.rhdf")
-    else
-      fname_fill="none"
-    fi
-    args+=(" resize_restart/fname=$fname resize_restart/use_dt=false parthenon/time/dt_min=$dt_new")
-    args+=(" resize_restart/fname_fill=$fname_fill ")
-  else
-    r_shell=$r_in
-    args+=(" bondi/r_shell=$r_shell ")
-  fi
-
-  # data_dir, logfiles
-  data_dir="${DR}/bondi_multizone_$(printf %05d ${VAR})"
-  out_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_out"
-  err_fn="${PDR}/logs/${DRTAG}/log_multizone$(printf %05d ${VAR})_err"
-
-  #srun --mpi=pmix ${PDR}/kharma_fork/kharma.cuda -i ${parfilename}  \
-  srun --mpi=pmix ${PDR}/kharma_faster_rst_fixed.cuda -i ${parfilename}  \
-                                    parthenon/job/problem_id=$prob \
-                                    parthenon/mesh/nx1=64 parthenon/mesh/nx2=64 parthenon/mesh/nx3=64 \
-                                    parthenon/meshblock/nx1=64 parthenon/meshblock/nx2=32 parthenon/meshblock/nx3=32 \
-                                    parthenon/time/tlim=${start_time} \
-                                    coordinates/r_in=${r_in} coordinates/r_out=${r_out} coordinates/a=$spin coordinates/ext_g=$EXT_G \
-                                    coordinates/transform=mks coordinates/hslope=1 \
-                                    GRMHD/reconstruction=linear_vl \
-                                    bondi/vacuum_logrho=${logrho} bondi/vacuum_log_u_over_rho=${log_u_over_rho} \
-                                    bondi/use_gizmo=true gizmo_shell/datfn=${datfn} \
-                                    b_field/type=none b_field/solver=none \
-                                    b_field/fix_flux_x1=0 b_field/initial_cleanup=0 \
-                                    resize_restart/base=$BASE resize_restart/nzone=$NZONES resize_restart/iteration=$iteration \
-                                    parthenon/output0/dt=$output0_dt \
-                                    parthenon/output1/dt=$output1_dt \
-                                    parthenon/output2/dt=$output2_dt \
-                                    ${args[@]} \
-                                    -d ${data_dir} 1> ${out_fn} 2>${err_fn}
-                                    # logrho=-5.34068635
-                                    #gizmo_shell/datfn="/n/holylfs05/LABS/bhi/Users/hyerincho/grmhd/data/gizmo/041423_3d/new_snap_deposit_020.hdf5" 
-
-  if [ $VAR -ne 0 ]; then
-    if [ $(($VAR % ($NZONES-1))) -eq 0 ]; then
-      out_to_in=$(($out_to_in*(-1)))
-      iteration=$(($iteration+1))
-    fi
-  fi
-
-  if [ $out_to_in -gt 0 ]; then
-    # half the radii
-    r_out=$((${r_out}/$BASE))
-    r_in=$((${r_in}/$BASE))
-  else
-    # double the radii
-    r_out=$((${r_out}*$BASE))
-    r_in=$((${r_in}*$BASE))
-  fi
-done