diff --git a/moment_kinetics/src/analysis.jl b/moment_kinetics/src/analysis.jl
index b6e59edc2..fc7722975 100644
--- a/moment_kinetics/src/analysis.jl
+++ b/moment_kinetics/src/analysis.jl
@@ -15,7 +15,7 @@ using ..interpolation: interpolate_to_grid_1d
 using ..load_data: open_readonly_output_file, get_nranks, load_pdf_data, load_rank_data
 using ..load_data: load_distributed_ion_pdf_slice
 using ..looping
-using ..type_definitions: mk_int, mk_float
+using ..type_definitions
 using ..velocity_moments: integrate_over_vspace
 
 using FFTW
@@ -690,7 +690,7 @@ function steady_state_square_residuals(variable, variable_at_previous_time, dt,
             end
             variable_max = MPI.Allreduce(local_max, max, comm_world)
         end
-        if isa(dt, MKOrAbtractVector)
+        if isa(dt, MKOrAbstractVector)
             reshaped_dt = reshape(dt, tuple((1 for _ ∈ 1:t_dim-1)..., size(dt)...))
         else
             reshaped_dt = dt
@@ -699,9 +699,9 @@ function steady_state_square_residuals(variable, variable_at_previous_time, dt,
             if size(dt) == ()
                 # local_max_absolute should always be at least a 1d array of size 1, not
                 # a 0d array, so that the MPI.Gather() below works correctly.
-                local_max_absolute = zeros(1)
+                local_max_absolute = mk_zeros(1)
             else
-                local_max_absolute = zeros(size(dt))
+                local_max_absolute = mk_zeros(size(dt))
             end
         else
             if size(dt) == ()
@@ -732,7 +732,7 @@ function steady_state_square_residuals(variable, variable_at_previous_time, dt,
                 # have size 1.
                 this_dims = tuple((1:t_dim-3)...)
                 if this_dims === ()
-                    local_max_absolute = max.(local_max_absolute, [absolute_residual])
+                    local_max_absolute[1] = max(local_max_absolute[1], absolute_residual[])
                 else
                     local_max_absolute = max.(local_max_absolute,
                                               vec(maximum(absolute_residual,
@@ -996,7 +996,7 @@ phi_fit_result struct whose fields are:
 """
 function fit_delta_phi_mode(t, z, delta_phi)
     # First fit a cosine to each time slice
-    results = allocate_float(3, size(delta_phi)[2])
+    results = zeros(mk_float, 3, size(delta_phi)[2])
     amplitude_guess = 1.0
     offset_guess = 0.0
     for (i, phi_z) in enumerate(eachcol(delta_phi))
diff --git a/moment_kinetics/src/chebyshev.jl b/moment_kinetics/src/chebyshev.jl
index 9fc5ab9f0..ea159c751 100644
--- a/moment_kinetics/src/chebyshev.jl
+++ b/moment_kinetics/src/chebyshev.jl
@@ -238,7 +238,7 @@ function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
     if irank == 0 # use a Chebyshev-Gauss-Radau element for the lowest element on rank 0
         scale_factor = element_scale[1]
         shift = element_shift[1]
-        grid[imin[1]:imax[1]] .= (chebyshev_radau_grid[1:ngrid] * scale_factor) .+ shift
+        grid[imin[1]:imax[1]] .= (chebyshev_radau_grid[1:ngrid] .* scale_factor) .+ shift
         # account for the fact that the minimum index needed for the chebyshev_grid
         # within each element changes from 1 to 2 in going from the first element
         # to the remaining elements
@@ -248,7 +248,7 @@ function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
             shift = element_shift[j]
             # reverse the order of the original chebyshev_grid (ran from [1,-1])
             # and apply the scale factor and shift
-            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
+            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] .* scale_factor) .+ shift
         end
         wgts = clenshaw_curtis_radau_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     else
@@ -261,7 +261,7 @@ function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
             shift = element_shift[j]
             # reverse the order of the original chebyshev_grid (ran from [1,-1])
             # and apply the scale factor and shift
-            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
+            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] .* scale_factor) .+ shift
             # after first element, increase minimum index for chebyshev_grid to 2
             # to avoid double-counting boundary element
             k = 2
@@ -578,10 +578,10 @@ function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_s
     @inbounds begin
         # calculate the weights within a single element and
         # scale to account for modified domain (not [-1,1])
-        wgts[1:ngrid] = w*element_scale[1]
+        wgts[1:ngrid] .= w .* element_scale[1]
         if nelement_local > 1
             for j ∈ 2:nelement_local
-                wgts[imin[j]-1:imax[j]] .+= w*element_scale[j]
+                wgts[imin[j]-1:imax[j]] .+= w .* element_scale[j]
             end
         end
     end
@@ -598,13 +598,13 @@ function clenshaw_curtis_radau_weights(ngrid, nelement_local, n, imin, imax, ele
     @inbounds begin
         # calculate the weights within a single element and
         # scale to account for modified domain (not [-1,1])
-        wgts[1:ngrid] .= wgts_radau[1:ngrid]*element_scale[1]
+        wgts[1:ngrid] .= wgts_radau[1:ngrid] .* element_scale[1]
         if nelement_local > 1
             for j ∈ 2:nelement_local
                 # account for double-counting of points at inner element boundaries
-                wgts[imin[j]-1] += wgts_lobatto[1]*element_scale[j]
+                wgts[imin[j]-1] += wgts_lobatto[1] .* element_scale[j]
                 # assign weights for interior of elements and one boundary point
-                wgts[imin[j]:imax[j]] .= wgts_lobatto[2:ngrid]*element_scale[j]
+                wgts[imin[j]:imax[j]] .= wgts_lobatto[2:ngrid] .* element_scale[j]
             end
         end
     end
@@ -616,7 +616,7 @@ compute and return modified Chebyshev moments of the first kind:
 ∫dx Tᵢ(x) over range [-1,1]
 """
 function chebyshevmoments(N)
-    μ = zeros(N)
+    μ = mk_zeros(N)
     @inbounds for i = 0:2:N-1
         μ[i+1] = 2/(1-i^2)
     end
diff --git a/moment_kinetics/src/clenshaw_curtis.jl b/moment_kinetics/src/clenshaw_curtis.jl
index df20dbb7d..8c97dd415 100644
--- a/moment_kinetics/src/clenshaw_curtis.jl
+++ b/moment_kinetics/src/clenshaw_curtis.jl
@@ -24,7 +24,7 @@ clenshawcurtisweights(μ) = clenshawcurtisweights!(copy(μ))
 clenshawcurtisweights!(μ) = clenshawcurtisweights!(μ, plan_clenshawcurtis(μ))
 function clenshawcurtisweights!(μ, plan)
     N = length(μ)
-    rmul!(μ, inv(N-1))
+    μ .*= inv(N-1)
     plan*μ
     μ[1] *= 0.5; μ[N] *= 0.5
     return μ
diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl
index 7b3544f08..08db47b41 100644
--- a/moment_kinetics/src/em_fields.jl
+++ b/moment_kinetics/src/em_fields.jl
@@ -129,10 +129,13 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
     ## calculate the electric fields after obtaining phi
     #Er = - d phi / dr 
     if r.n > 1
-        derivative_r!(fields.Er,-fields.phi,
+        derivative_r!(fields.Er,fields.phi,
                 scratch_dummy.buffer_z_1, scratch_dummy.buffer_z_2,
                 scratch_dummy.buffer_z_3, scratch_dummy.buffer_z_4,
                 r_spectral,r)
+        @loop_r_z ir iz begin
+            fields.Er[iz,ir] *= -1.0
+        end
         if z.irank == 0 && fields.force_Er_zero_at_wall
             fields.Er[1,:] .= 0.0
         end
@@ -150,10 +153,13 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
                                        kinetic_electrons_with_temperature_equation)
         if z.n > 1
             # Ez = - d phi / dz
-            @views derivative_z!(fields.Ez,-fields.phi,
+            @views derivative_z!(fields.Ez, fields.phi,
                     scratch_dummy.buffer_rs_1[:,1], scratch_dummy.buffer_rs_2[:,1],
                     scratch_dummy.buffer_rs_3[:,1], scratch_dummy.buffer_rs_4[:,1],
                     z_spectral,z)
+            @loop_r_z ir iz begin
+                fields.Ez[iz,ir] *= -1.0
+            end
         end
     end
 
diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl
index 6d92d8014..e9e6641cf 100644
--- a/moment_kinetics/src/file_io.jl
+++ b/moment_kinetics/src/file_io.jl
@@ -2710,7 +2710,7 @@ function write_timing_data(io_moments, t_idx, dfns=false)
 
         # Gather the strings
         gathered_char_vector = MKVector{Char}(undef, sum(string_sizes))
-        local_start_index = sum(string_sizes[1:comm_rank]) + 1
+        local_start_index = sum(@view string_sizes[1:comm_rank]) + 1
         local_end_index = local_start_index - 1 + string_sizes[comm_rank+1]
         gathered_char_vector[local_start_index:local_end_index] .= [new_names_string...]
         gathered_buffer = MPI.VBuffer(gathered_char_vector, string_sizes)
diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl
index 52eb838d5..a6d023fbc 100644
--- a/moment_kinetics/src/gauss_legendre.jl
+++ b/moment_kinetics/src/gauss_legendre.jl
@@ -1146,7 +1146,7 @@ function setup_global_strong_form_matrix!(QQ_global::MKOrAbstractArray{mk_float,
         # upper boundary assembly on element
         if j == coord.nelement_local
             if periodic_bc && coord.nrank == 1
-                QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] / 2.0
+                QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] ./ 2.0
             else
                 QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:]
             end
diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl
index 0a54dbe41..b6ce7c862 100644
--- a/moment_kinetics/src/load_data.jl
+++ b/moment_kinetics/src/load_data.jl
@@ -3819,9 +3819,11 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
             !isa(it, mk_int) && push!(dims, nt)
             vartype = typeof(variable[1][1])
             if vartype == mk_int
-                result = allocate_int(dims...)
+                #result = allocate_int(dims...)
+                result = zeros(mk_int, dims...)
             elseif vartype == mk_float
-                result = allocate_float(dims...)
+                #result = allocate_float(dims...)
+                result = zeros(mk_float, dims...)
             else
                 error("Unsupported dtype for 1D variable $(variable.dtype)")
             end
@@ -3835,9 +3837,11 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
             !isa(it, mk_int) && push!(dims, nt)
             vartype = typeof(variable[1][1,1])
             if vartype == mk_int
-                result = allocate_int(dims...)
+                #result = allocate_int(dims...)
+                result = zeros(mk_int, dims...)
             elseif vartype == mk_float
-                result = allocate_float(dims...)
+                #result = allocate_float(dims...)
+                result = zeros(mk_float, dims...)
             else
                 error("Unsupported dtype for 1D variable $(variable.dtype)")
             end
@@ -3847,7 +3851,8 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
             !isa(iz, mk_int) && push!(dims, nz)
             !isa(ir, mk_int) && push!(dims, nr)
             !isa(it, mk_int) && push!(dims, nt)
-            result = allocate_float(dims...)
+            #result = allocate_float(dims...)
+            result = zeros(mk_float, dims...)
         elseif nd == 4
             # moment variable with dimensions (z,r,s,t)
             # Get nspecies from the variable, not from run_info, because it might be
@@ -3862,7 +3867,8 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
                 push!(dims, nspecies)
             end
             !isa(it, mk_int) && push!(dims, nt)
-            result = allocate_float(dims...)
+            #result = allocate_float(dims...)
+            result = zeros(mk_float, dims...)
         elseif nd == 5
             # electron distribution function variable with dimensions (vpa,vperp,z,r,t)
             dims = Vector{mk_int}()
@@ -3871,7 +3877,8 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
             !isa(iz, mk_int) && push!(dims, nz)
             !isa(ir, mk_int) && push!(dims, nr)
             !isa(it, mk_int) && push!(dims, nt)
-            result = allocate_float(dims...)
+            #result = allocate_float(dims...)
+            result = zeros(mk_float, dims...)
         elseif nd == 6
             # ion distribution function variable with dimensions (vpa,vperp,z,r,s,t)
             nspecies = size(variable[1], 5)
@@ -3886,7 +3893,8 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
                 push!(dims, nspecies)
             end
             !isa(it, mk_int) && push!(dims, nt)
-            result = allocate_float(dims...)
+            #result = allocate_float(dims...)
+            result = zeros(mk_float, dims...)
         elseif nd == 7
             # neutral distribution function variable with dimensions (vz,vr,vzeta,z,r,s,t)
             nspecies = size(variable[1], 6)
@@ -3902,7 +3910,8 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing,
                 push!(dims, nspecies)
             end
             !isa(it, mk_int) && push!(dims, nt)
-            result = allocate_float(dims...)
+            #result = allocate_float(dims...)
+            result = zeros(mk_float, dims...)
         else
             error("Unsupported number of dimensions ($nd) for '$variable_name'.")
         end
diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl
index d0398a1d4..d2f13178e 100644
--- a/moment_kinetics/src/moment_kinetics.jl
+++ b/moment_kinetics/src/moment_kinetics.jl
@@ -4,8 +4,13 @@ module moment_kinetics
 
 export run_moment_kinetics
 
+using HDF5
 using MPI
 using StatsBase
+try
+    using NCDatasets
+catch
+end
 
 # Include submodules from other source files
 # Note that order of includes matters - things used in one module must already
diff --git a/moment_kinetics/src/nonlinear_solvers.jl b/moment_kinetics/src/nonlinear_solvers.jl
index 63c5c1690..6338d4300 100644
--- a/moment_kinetics/src/nonlinear_solvers.jl
+++ b/moment_kinetics/src/nonlinear_solvers.jl
@@ -564,7 +564,7 @@ old_precon_iterations = nl_solver_params.precon_iterations[]
 end
 
 @timeit_debug global_timer distributed_norm(
-                               ::Val{:z}, residual::MKArray{mk_float, 1}, coords,
+                               ::Val{:z}, residual::AbstractMKArray{mk_float, 1}, coords,
                                rtol, atol, x) = begin
     z = coords.z
 
@@ -600,7 +600,7 @@ end
 end
 
 @timeit_debug global_timer distributed_norm(
-                               ::Val{:vpa}, residual::MKArray{mk_float, 1}, coords,
+                               ::Val{:vpa}, residual::AbstractMKArray{mk_float, 1}, coords,
                                rtol, atol, x) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
@@ -616,7 +616,7 @@ end
 
 @timeit_debug global_timer distributed_norm(
                                ::Val{:zvperpvpa},
-                               residual::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}},
+                               residual::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}},
                                coords, rtol, atol, x) = begin
     ppar_residual, pdf_residual = residual
     x_ppar, x_pdf = x
@@ -679,7 +679,7 @@ end
 end
 
 @timeit_debug global_timer distributed_norm(
-                  ::Val{:srzvperpvpa}, residual::MKArray{mk_float, 5}, coords, rtol,
+                  ::Val{:srzvperpvpa}, residual::AbstractMKArray{mk_float, 5}, coords, rtol,
                   atol, x) = begin
     n_ion_species = coords.s
     r = coords.r
@@ -724,7 +724,7 @@ end
 end
 
 @timeit_debug global_timer distributed_dot(
-                  ::Val{:z}, v::MKArray{mk_float, 1}, w::MKArray{mk_float, 1},
+                  ::Val{:z}, v::AbstractMKArray{mk_float, 1}, w::AbstractMKArray{mk_float, 1},
                   coords, rtol, atol, x) = begin
 
     z = coords.z
@@ -761,7 +761,7 @@ end
 end
 
 @timeit_debug global_timer distributed_dot(
-                  ::Val{:vpa}, v::MKArray{mk_float, 1}, w::MKArray{mk_float, 1}, coords,
+                  ::Val{:vpa}, v::AbstractMKArray{mk_float, 1}, w::AbstractMKArray{mk_float, 1}, coords,
                   rtol, atol, x) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
@@ -774,8 +774,8 @@ end
 end
 
 @timeit_debug global_timer distributed_dot(
-                  ::Val{:zvperpvpa}, v::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}},
-                  w::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}}, coords,
+                  ::Val{:zvperpvpa}, v::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}},
+                  w::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}}, coords,
                   rtol, atol, x) = begin
     v_ppar, v_pdf = v
     w_ppar, w_pdf = w
@@ -835,8 +835,8 @@ end
 end
 
 @timeit_debug global_timer distributed_dot(
-                  ::Val{:srzvperpvpa}, v::MKArray{mk_float, 5},
-                  w::MKArray{mk_float, 5}, coords, rtol, atol, x) = begin
+                  ::Val{:srzvperpvpa}, v::AbstractMKArray{mk_float, 5},
+                  w::AbstractMKArray{mk_float, 5}, coords, rtol, atol, x) = begin
     n_ion_species = coords.s
     r = coords.r
     z = coords.z
@@ -880,7 +880,7 @@ end
 # slow code
 
 @timeit_debug global_timer parallel_map(
-                  ::Val{:z}, func, result::MKArray{mk_float, 1}) = begin
+                  ::Val{:z}, func, result::AbstractMKArray{mk_float, 1}) = begin
 
     begin_z_region()
 
@@ -891,7 +891,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:z}, func, result::MKArray{mk_float, 1}, x1) = begin
+                  ::Val{:z}, func, result::AbstractMKArray{mk_float, 1}, x1) = begin
 
     begin_z_region()
 
@@ -902,11 +902,11 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:z}, func, result::MKArray{mk_float, 1}, x1, x2) = begin
+                  ::Val{:z}, func, result::AbstractMKArray{mk_float, 1}, x1, x2) = begin
 
     begin_z_region()
 
-    if isa(x2, MKArray)
+    if isa(x2, AbstractMKArray)
         @loop_z iz begin
             result[iz] = func(x1[iz], x2[iz])
         end
@@ -919,11 +919,11 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:z}, func, result::MKArray{mk_float, 1}, x1, x2, x3) = begin
+                  ::Val{:z}, func, result::AbstractMKArray{mk_float, 1}, x1, x2, x3) = begin
 
     begin_z_region()
 
-    if isa(x3, MKArray)
+    if isa(x3, AbstractMKArray)
         @loop_z iz begin
             result[iz] = func(x1[iz], x2[iz], x3[iz])
         end
@@ -937,7 +937,7 @@ end
 end
 
 @timeit_debug global_timer parallel_map(
-                  ::Val{:vpa}, func, result::MKArray{mk_float, 1}) = begin
+                  ::Val{:vpa}, func, result::AbstractMKArray{mk_float, 1}) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
     for i ∈ eachindex(result)
@@ -946,7 +946,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:vpa}, func, result::MKArray{mk_float, 1}, x1) = begin
+                  ::Val{:vpa}, func, result::AbstractMKArray{mk_float, 1}, x1) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
     for i ∈ eachindex(result)
@@ -955,10 +955,10 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:vpa}, func, result::MKArray{mk_float, 1}, x1, x2) = begin
+                  ::Val{:vpa}, func, result::AbstractMKArray{mk_float, 1}, x1, x2) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
-    if isa(x2, MKArray)
+    if isa(x2, AbstractMKArray)
         for i ∈ eachindex(result)
             result[i] = func(x1[i], x2[i])
         end
@@ -970,10 +970,10 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:vpa}, func, result::MKArray{mk_float, 1}, x1, x2, x3) = begin
+                  ::Val{:vpa}, func, result::AbstractMKArray{mk_float, 1}, x1, x2, x3) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
-    if isa(x3, MKArray)
+    if isa(x3, AbstractMKArray)
         for i ∈ eachindex(result)
             result[i] = func(x1[i], x2[i], x3[i])
         end
@@ -986,7 +986,7 @@ end
 end
 
 @timeit_debug global_timer parallel_map(
-                  ::Val{:zvperpvpa}, func, result::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}}) = begin
+                  ::Val{:zvperpvpa}, func, result::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}}) = begin
 
     result_ppar, result_pdf = result
 
@@ -1005,7 +1005,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:zvperpvpa}, func, result::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}},
+                  ::Val{:zvperpvpa}, func, result::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}},
                   x1) = begin
 
     result_ppar, result_pdf = result
@@ -1026,7 +1026,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:zvperpvpa}, func, result::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}},
+                  ::Val{:zvperpvpa}, func, result::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}},
                   x1, x2) = begin
 
     result_ppar, result_pdf = result
@@ -1062,7 +1062,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:zvperpvpa}, func, result::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}},
+                  ::Val{:zvperpvpa}, func, result::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}},
                   x1, x2, x3) = begin
 
     result_ppar, result_pdf = result
@@ -1100,7 +1100,7 @@ end
 end
 
 @timeit_debug global_timer parallel_map(
-                  ::Val{:srzvperpvpa}, func, result::MKArray{mk_float, 5}) = begin
+                  ::Val{:srzvperpvpa}, func, result::AbstractMKArray{mk_float, 5}) = begin
 
     begin_s_r_z_vperp_vpa_region()
 
@@ -1111,7 +1111,7 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:srzvperpvpa}, func, result::MKArray{mk_float, 5}, x1) = begin
+                  ::Val{:srzvperpvpa}, func, result::AbstractMKArray{mk_float, 5}, x1) = begin
 
     begin_s_r_z_vperp_vpa_region()
 
@@ -1122,11 +1122,11 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:srzvperpvpa}, func, result::MKArray{mk_float, 5}, x1, x2) = begin
+                  ::Val{:srzvperpvpa}, func, result::AbstractMKArray{mk_float, 5}, x1, x2) = begin
 
     begin_s_r_z_vperp_vpa_region()
 
-    if isa(x2, MKArray)
+    if isa(x2, AbstractMKArray)
         @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
             result[ivpa,ivperp,iz,ir,is] = func(x1[ivpa,ivperp,iz,ir,is], x2[ivpa,ivperp,iz,ir,is])
         end
@@ -1139,12 +1139,12 @@ end
     return nothing
 end
 @timeit_debug global_timer parallel_map(
-                  ::Val{:srzvperpvpa}, func, result::MKArray{mk_float, 5}, x1, x2,
+                  ::Val{:srzvperpvpa}, func, result::AbstractMKArray{mk_float, 5}, x1, x2,
                   x3) = begin
 
     begin_s_r_z_vperp_vpa_region()
 
-    if isa(x3, MKArray)
+    if isa(x3, AbstractMKArray)
         @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
             result[ivpa,ivperp,iz,ir,is] = func(x1[ivpa,ivperp,iz,ir,is], x2[ivpa,ivperp,iz,ir,is], x3[ivpa,ivperp,iz,ir,is])
         end
@@ -1158,7 +1158,7 @@ end
 end
 
 @timeit_debug global_timer parallel_delta_x_calc(
-                  ::Val{:z}, delta_x::MKArray{mk_float, 1}, V, y) = begin
+                  ::Val{:z}, delta_x::AbstractMKArray{mk_float, 1}, V, y) = begin
 
     begin_z_region()
 
@@ -1173,7 +1173,7 @@ end
 end
 
 @timeit_debug global_timer parallel_delta_x_calc(
-                  ::Val{:vpa}, delta_x::MKArray{mk_float, 1}, V, y) = begin
+                  ::Val{:vpa}, delta_x::AbstractMKArray{mk_float, 1}, V, y) = begin
     # No parallelism needed when the implicit solve is over vpa - assume that this will be
     # called inside a parallelised s_r_z_vperp loop.
     ny = length(y)
@@ -1186,7 +1186,7 @@ end
 end
 
 @timeit_debug global_timer parallel_delta_x_calc(
-                  ::Val{:zvperpvpa}, delta_x::Tuple{MKArray{mk_float, 1},MKArray{mk_float, 3}}, V,
+                  ::Val{:zvperpvpa}, delta_x::Tuple{AbstractMKArray{mk_float, 1},AbstractMKArray{mk_float, 3}}, V,
                   y) = begin
 
     delta_x_ppar, delta_x_pdf = delta_x
@@ -1214,7 +1214,7 @@ end
 end
 
 @timeit_debug global_timer parallel_delta_x_calc(
-                  ::Val{:srzvperpvpa}, delta_x::MKArray{mk_float, 5}, V, y) = begin
+                  ::Val{:srzvperpvpa}, delta_x::AbstractMKArray{mk_float, 5}, V, y) = begin
 
     begin_s_r_z_vperp_vpa_region()
 
diff --git a/moment_kinetics/src/type_definitions.jl b/moment_kinetics/src/type_definitions.jl
index bc248404c..2ca0f5acf 100644
--- a/moment_kinetics/src/type_definitions.jl
+++ b/moment_kinetics/src/type_definitions.jl
@@ -8,6 +8,9 @@ export mk_zeros
 export MKArray
 export MKVector
 export MKMatrix
+export AbstractMKArray
+export AbstractMKVector
+export AbstractMKMatrix
 export MKOrAbstractArray
 export MKOrAbstractVector
 export MKOrAbstractMatrix
@@ -36,6 +39,18 @@ const MKVector{T} = MKArray{T,1}
 """
 const MKMatrix{T} = MKArray{T,2}
 
+"""
+"""
+const AbstractMKArray{T,N} = AbstractInboundsArray{T,N}
+
+"""
+"""
+const AbstractMKVector{T} = AbstractInboundsArray{T,1}
+
+"""
+"""
+const AbstractMKMatrix{T} = AbstractInboundsArray{T,2}
+
 """
 """
 const MKOrAbstractArray{T,N} = Union{AbstractInboundsArray{T,N}, AbstractArray{T,N}} where {T,N}
diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl
index 23a77e71a..4ea470b2f 100644
--- a/moment_kinetics/test/calculus_tests.jl
+++ b/moment_kinetics/test/calculus_tests.jl
@@ -2,6 +2,7 @@ module CalculusTests
 
 include("setup.jl")
 
+using moment_kinetics.array_allocation: allocate_float
 using moment_kinetics.coordinates: define_test_coordinate
 using moment_kinetics.calculus: derivative!, second_derivative!, integral
 using moment_kinetics.calculus: laplacian_derivative!
@@ -47,9 +48,9 @@ function runtests()
                                                      element_spacing_option=element_spacing_option,
                                                      collision_operator_dim=false)
                 # create array for the function f(x) to be differentiated/integrated
-                f = Array{Float64,1}(undef, x.n)
+                f = allocate_float(x.n)
                 # create array for the derivative df/dx
-                df = Array{Float64,1}(undef, x.n)
+                df = allocate_float(x.n)
                 # initialize f
                 for ix ∈ 1:x.n
                     f[ix] = ( (cospi(2.0*x.grid[ix]/x.L)+sinpi(2.0*x.grid[ix]/x.L))
@@ -94,7 +95,7 @@ function runtests()
                                                      collision_operator_dim=false)
 
                 # create array for the derivative df/dx and the expected result
-                df = Array{Float64,1}(undef, x.n)
+                df = allocate_float(x.n)
 
                 # initialize f and expected df
                 offset = randn(rng)
@@ -142,7 +143,7 @@ function runtests()
                                                      collision_operator_dim=false)
 
                 # create array for the derivative df/dx and the expected result
-                df = Array{Float64,1}(undef, x.n)
+                df = allocate_float(x.n)
 
                 # initialize f and expected df
                 offset = randn(rng)
@@ -186,7 +187,7 @@ function runtests()
                                                      collision_operator_dim=false)
 
                 # create array for the derivative df/dx and the expected result
-                df = Array{Float64,1}(undef, x.n)
+                df = allocate_float(x.n)
 
                 # initialize f and expected df
                 offset = undef
@@ -236,7 +237,7 @@ function runtests()
                                                      collision_operator_dim=false)
 
                 # create array for the derivative df/dx and the expected result
-                df = Array{Float64,1}(undef, x.n)
+                df = allocate_float(x.n)
 
                 # initialize f and expected df
                 offset = undef
@@ -682,7 +683,7 @@ function runtests()
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
                     # create array for the function f(x) to be differentiated/integrated
-                    f = Array{Float64,1}(undef, x.n)
+                    f = allocate_float(x.n)
                     # create array for the derivative df/dx and the expected result
                     df = similar(f)
                     expected_df = similar(f)
@@ -728,7 +729,7 @@ function runtests()
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
                     # create array for the function f(x) to be differentiated/integrated
-                    f = Array{Float64,1}(undef, x.n)
+                    f = allocate_float(x.n)
                     # create array for the derivative df/dx and the expected result
                     df = similar(f)
                     expected_df = similar(f)
@@ -1009,7 +1010,7 @@ function runtests()
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
                     # create array for the function f(x) to be differentiated/integrated
-                    f = Array{Float64,1}(undef, x.n)
+                    f = allocate_float(x.n)
                     # create array for the derivative df/dx and the expected result
                     df = similar(f)
                     expected_df = similar(f)
@@ -1054,7 +1055,7 @@ function runtests()
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
                     # create array for the function f(x) to be differentiated/integrated
-                    f = Array{Float64,1}(undef, x.n)
+                    f = allocate_float(x.n)
                     # create array for the derivative df/dx and the expected result
                     df = similar(f)
                     expected_df = similar(f)
diff --git a/moment_kinetics/test/jacobian_matrix_tests.jl b/moment_kinetics/test/jacobian_matrix_tests.jl
index 6721f7b30..370681afc 100644
--- a/moment_kinetics/test/jacobian_matrix_tests.jl
+++ b/moment_kinetics/test/jacobian_matrix_tests.jl
@@ -231,7 +231,7 @@ function cleanup_mk_state!(args...)
     return nothing
 end
 
-function generate_norm_factor(perturbed_residual::AbstractArray{mk_float,3})
+function generate_norm_factor(perturbed_residual::MKArray{mk_float,3})
     # half-width of the window for moving average
     w = 3
     norm_factor_unsmoothed = mean(abs.(perturbed_residual); dims=3)
@@ -249,7 +249,7 @@ function generate_norm_factor(perturbed_residual::AbstractArray{mk_float,3})
     end
     return norm_factor
 end
-function generate_norm_factor(perturbed_residual::AbstractArray{mk_float,1})
+function generate_norm_factor(perturbed_residual::MKArray{mk_float,1})
     norm_factor_unsmoothed = mean(abs.(perturbed_residual); dims=1)
 end
 
diff --git a/moment_kinetics/test/nonlinear_solver_tests.jl b/moment_kinetics/test/nonlinear_solver_tests.jl
index 81d4c87c7..6f5d123f7 100644
--- a/moment_kinetics/test/nonlinear_solver_tests.jl
+++ b/moment_kinetics/test/nonlinear_solver_tests.jl
@@ -70,11 +70,11 @@ function linear_test()
 
         function rhs_func!(residual, x; krylov=false)
             if serial_solve
-                residual .= A * x - b
+                residual .= A * x .- b
             else
                 begin_serial_region()
                 @serial_region begin
-                    residual .= A * x - b
+                    residual .= A * x .- b
                 end
             end
             return nothing