diff --git a/examples/gk-ions/2D-periodic-gk.toml b/examples/gk-ions/2D-periodic-gk.toml
index 8e47f20b7..99fcbed29 100644
--- a/examples/gk-ions/2D-periodic-gk.toml
+++ b/examples/gk-ions/2D-periodic-gk.toml
@@ -45,7 +45,7 @@ discretization = "chebyshev_pseudospectral"
 n_ion_species = 1
 n_neutral_species = 0
 electron_physics = "boltzmann_electron_response"
-gyrokinetic_ions = true
+ion_physics = "gyrokinetic_ions"
 T_e = 1.0
 T_wall = 1.0
 
diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl
index 1034df1a0..42a39ef06 100644
--- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl
+++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl
@@ -261,6 +261,12 @@ function makie_post_process(run_dir::Union{String,Tuple},
             mkpath(comparison_plot_dir)
             plot_prefix = joinpath(comparison_plot_dir, "compare_")
         end
+    else
+        if length(run_info) != 1
+            comparison_plot_dir = "comparison_plots_$plot_prefix"
+            mkpath(comparison_plot_dir)
+            plot_prefix = joinpath(comparison_plot_dir, "compare_")
+        end
     end
 
     timestep_diagnostics(run_info, run_info_dfns; plot_prefix=plot_prefix)
@@ -349,6 +355,8 @@ function makie_post_process(run_dir::Union{String,Tuple},
 
     sound_wave_plots(run_info; plot_prefix=plot_prefix)
 
+    collisionality_plots(run_info, plot_prefix)
+
     if all(ri === nothing for ri ∈ run_info_dfns)
         nvperp = nothing
     else
@@ -818,7 +826,29 @@ function _setup_single_input!(this_input_dict::OrderedDict{String,Any},
         plot_steady_state_residual=false,
         animate_steady_state_residual=false,
        )
-
+    
+    set_defaults_and_check_section!(
+        this_input_dict, "collisionality_plots";
+        plot=true,
+        plot_dT_dz_vs_z=false,
+        animate_dT_dz_vs_z=false,
+        plot_mfp_vs_z=false,
+        animate_mfp_vs_z=false,
+        plot_nu_ii_vth_mfp_vs_z = false,
+        plot_LT_mfp_vs_z = false,
+        animate_LT_mfp_vs_z = false,
+        plot_LT_dT_dz_temp_vs_z = false,
+        plot_Ln_mfp_vs_z = false,
+        animate_Ln_mfp_vs_z = false,
+        plot_Lupar_mfp_vs_z = false,
+        animate_Lupar_mfp_vs_z = false,
+        plot_Lupar_Ln_LT_mfp_vs_z = false,
+        animate_Lupar_Ln_LT_mfp_vs_z = false,
+        plot_overlay_coll_krook_heat_flux = false,
+        animate_overlay_coll_krook_heat_flux = false,
+        animation_ext = "gif"
+       )
+       
     set_defaults_and_check_section!(
         this_input_dict, "timing_data";
         plot=false,
@@ -8426,6 +8456,357 @@ function timestep_diagnostics(run_info, run_info_dfns; plot_prefix=nothing, it=n
     return steps_fig, dt_fig, CFL_fig
 end
 
+"""
+A function to plot collisionalities. The mean free path is plotted (or animated) 
+along with the lengthscales of the gradients of density, parallel flow and temperature.
+
+There are also functions to check the calculations of the mean free path and the 
+comparison of temperature, L_T and dT_dz. They would only be for making sure
+lengthscales and mean free path calculations are sensible.
+"""
+function collisionality_plots(run_info, plot_prefix=nothing)
+    if !isa(run_info, Tuple)
+        run_info = (run_info,)
+    end
+    input = Dict_to_NamedTuple(input_dict["collisionality_plots"])
+
+    if input.plot
+        println("Making plots for collisionality")
+
+        temperature_input = Dict_to_NamedTuple(input_dict["temperature"])
+        mfp = get_variable(run_info, "mfp")
+        L_T = get_variable(run_info, "L_T")
+        L_n = get_variable(run_info, "L_n")
+        L_upar = get_variable(run_info, "L_upar")
+        nt = minimum(length(mfp[ri][1,1,1,:]) for ri in eachindex(run_info))
+        # print warning if the lengths of all the mfp[ri][1,1,1,:] are not the same
+        if any(length(mfp[ri][1,1,1,:]) != nt for ri in eachindex(run_info))
+            println("Warning: The number of timesteps of some simulations are different, " *
+                    "so only the first common timesteps will be animated.")
+        end
+
+        # write function to check that mfp[ri][1, 1, 1, :] is the same length (i.e. nt) for all ri 
+        if input.plot_mfp_vs_z
+            variable_name = "mean_free_path"
+            variable = mfp
+            variable_prefix = plot_prefix * variable_name
+            plot_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                    outfile=variable_prefix * "_vs_z.pdf")
+        end
+
+        if input.animate_mfp_vs_z
+            variable_name = "mean_free_path"
+            variable = mfp
+            variable_prefix = plot_prefix * variable_name
+            animate_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                            outfile=variable_prefix * "_vs_z." * input.animation_ext)
+        end
+
+        if input.plot_dT_dz_vs_z
+            variable_name = "dT_dz"
+            variable = nothing
+            try
+                variable = get_variable(run_info, variable_name)
+            catch e
+                return makie_post_processing_error_handler(
+                        e,
+                        "collisionality_plots () failed for $variable_name - could not load data.")
+            end
+
+            variable_prefix = plot_prefix * variable_name
+            plot_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                    outfile=variable_prefix * "_vs_z.pdf")
+        end
+
+        if input.animate_dT_dz_vs_z
+            variable_name = "dT_dz"
+            variable = nothing
+            try
+                variable = get_variable(run_info, variable_name)
+            catch e
+                return makie_post_processing_error_handler(
+                        e,
+                        "collisionality_plots () failed for $variable_name - could not load data.")
+            end
+
+            variable_prefix = plot_prefix * variable_name
+            animate_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                            outfile=variable_prefix * "_vs_z." * input.animation_ext)
+        end
+
+        if input.plot_nu_ii_vth_mfp_vs_z
+            vth = get_variable(run_info, "thermal_speed")
+            nu_ii = get_variable(run_info, "collision_frequency_ii")
+            variable_prefix = plot_prefix * "checking_mfp_vth_and_nu_ii"
+
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, vth[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"vth", ax=ax[1], title = "checking_mfp_vth")
+                plot_1d(run_info[ri].z.grid, nu_ii[ri][:,1,1,end], label=run_label*"nu_ii", ax=ax[1])
+                plot_1d(run_info[ri].z.grid, mfp[ri][:,1,1,end], label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.plot_LT_dT_dz_temp_vs_z
+            dT_dz = get_variable(run_info, "dT_dz")
+            temp = get_variable(run_info, "temperature")
+            variable_prefix = plot_prefix * "LT_dT_dz_temp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, L_T[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"L_T", ax=ax[1], title = "LT_dT_dz_temp")
+                plot_1d(run_info[ri].z.grid, dT_dz[ri][:,1,1,end], label=run_label*"dT_dz", ax=ax[1])
+                plot_1d(run_info[ri].z.grid, temp[ri][:,1,1,end], label=run_label*"temp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.plot_LT_mfp_vs_z
+            variable_prefix = plot_prefix * "LT_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, L_T[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"L_T", ax=ax[1], title = "L_T and mean free path comparison")
+                plot_1d(run_info[ri].z.grid, mfp[ri][:,1,1,end], label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.animate_LT_mfp_vs_z
+            variable_prefix = plot_prefix * "LT_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            frame_index = Observable(1)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                animate_1d(run_info[ri].z.grid, L_T[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_T", ax=ax[1], title = "L_T and mean free path comparison")
+                animate_1d(run_info[ri].z.grid, mfp[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z." * input.animation_ext
+            save_animation(fig, frame_index, nt, outfile)
+        end
+
+        if input.plot_Ln_mfp_vs_z
+            variable_prefix = plot_prefix * "Ln_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, L_n[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"L_n", ax=ax[1], title = "Ln and mean free path comparison")
+                plot_1d(run_info[ri].z.grid, mfp[ri][:,1,1,end], label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.animate_Ln_mfp_vs_z
+            variable_prefix = plot_prefix * "Ln_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            frame_index = Observable(1)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                animate_1d(run_info[ri].z.grid, L_n[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_n", ax=ax[1], title = "L_n and mean free path comparison")
+                animate_1d(run_info[ri].z.grid, mfp[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z." * input.animation_ext
+            save_animation(fig, frame_index, nt, outfile)
+        end
+
+        if input.plot_Lupar_mfp_vs_z
+            variable_prefix = plot_prefix * "Lupar_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, L_upar[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"L_upar", ax=ax[1], title = "Lupar and mean free path comparison")
+                plot_1d(run_info[ri].z.grid, mfp[ri][:,1,1,end], label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.animate_Lupar_mfp_vs_z
+            variable_prefix = plot_prefix * "Lupar_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            frame_index = Observable(1)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                animate_1d(run_info[ri].z.grid, L_upar[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_upar", ax=ax[1], title = "L_upar and mean free path comparison")
+                animate_1d(run_info[ri].z.grid, mfp[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z." * input.animation_ext
+            save_animation(fig, frame_index, nt, outfile)
+        end
+
+        if input.plot_Lupar_Ln_LT_mfp_vs_z
+            variable_prefix = plot_prefix * "Lupar_Ln_LT_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                plot_1d(run_info[ri].z.grid, L_upar[ri][:,1,1,end], xlabel="z",
+                        ylabel="values", label=run_label*"L_upar", ax=ax[1], title = "Lupar Ln LT and mean free path comparison")
+                plot_1d(run_info[ri].z.grid, L_n[ri][:,1,1,end], label=run_label*"L_n", ax=ax[1])
+                plot_1d(run_info[ri].z.grid, L_T[ri][:,1,1,end], label=run_label*"L_T", ax=ax[1])
+                plot_1d(run_info[ri].z.grid, mfp[ri][:,1,1,end], label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:horizontal)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.animate_Lupar_Ln_LT_mfp_vs_z
+            variable_prefix = plot_prefix * "Lupar_Ln_LT_mfp"
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            frame_index = Observable(1)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                animate_1d(run_info[ri].z.grid, L_upar[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_upar", ax=ax[1], title = "Lupar Ln LT and mean free path comparison")
+                animate_1d(run_info[ri].z.grid, L_n[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_n", ax=ax[1])
+                animate_1d(run_info[ri].z.grid, L_T[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"L_T", ax=ax[1])
+                animate_1d(run_info[ri].z.grid, mfp[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"mfp", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z." * input.animation_ext
+            save_animation(fig, frame_index, nt, outfile)
+        end
+
+        if input.plot_overlay_coll_krook_heat_flux
+            variable_prefix = plot_prefix * "coll_krook_vs_original_heat_flux"
+            coll_krook_q = get_variable(run_info, "coll_krook_heat_flux")
+            original_q = get_variable(run_info, "parallel_heat_flux")
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                if get(run_info[ri].input["composition"], "ion_physics", "") !== coll_krook_ions
+                    plot_1d(run_info[ri].z.grid, coll_krook_q[ri][:,1,1,end], xlabel="z",
+                            ylabel="values", label=run_label*"coll_krook_q_overlay", ax=ax[1], title = "coll_krook heat flux overlay")
+                end
+                plot_1d(run_info[ri].z.grid, original_q[ri][:,1,1,end], label=run_label*"original_q", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z.pdf"
+            save(outfile, fig)
+        end
+
+        if input.animate_overlay_coll_krook_heat_flux
+            variable_prefix = plot_prefix * "coll_krook_vs_original_heat_flux"
+            coll_krook_q = get_variable(run_info, "coll_krook_heat_flux")
+            original_q = get_variable(run_info, "parallel_heat_flux")
+            fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+            frame_index = Observable(1)
+            for ri ∈ eachindex(run_info)
+                if length(run_info) > 1
+                    run_label = run_info[ri].run_name * " "
+                else
+                    run_label = " "
+                end
+                if get(run_info[ri].input["composition"], "ion_physics", "") !== coll_krook_ions
+                    animate_1d(run_info[ri].z.grid, coll_krook_q[ri][:,1,1,:],
+                            frame_index=frame_index, xlabel="z", ylabel="values",
+                            label=run_label*"coll_krook_q_overlay", ax=ax[1], title = "coll_krook heat flux overlay")
+                end
+                animate_1d(run_info[ri].z.grid, original_q[ri][:,1,1,:],
+                        frame_index=frame_index, xlabel="z", ylabel="values",
+                        label=run_label*"original_q", ax=ax[1])
+            end
+            Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+                orientation=:vertical)
+            outfile = variable_prefix * "_vs_z." * input.animation_ext
+            save_animation(fig, frame_index, nt, outfile)
+        end
+    end
+end
+
 """
     timing_data(run_info; plot_prefix=nothing, threshold=nothing,
                 include_patterns=nothing, exclude_patterns=nothing, ranks=nothing,
diff --git a/moment_kinetics/debug_test/gyroaverage_inputs.jl b/moment_kinetics/debug_test/gyroaverage_inputs.jl
index 9855050c3..39c544d50 100644
--- a/moment_kinetics/debug_test/gyroaverage_inputs.jl
+++ b/moment_kinetics/debug_test/gyroaverage_inputs.jl
@@ -6,7 +6,7 @@ test_input = OptionsDict(
     "output" => OptionsDict("run_name" => "gyroaverage"),
     "composition" => OptionsDict("n_ion_species" => 1,
                           "n_neutral_species" => 0,
-                          "gyrokinetic_ions" => true,
+                          "ion_physics" => "gyrokinetic_ions",
                           "T_e" => 1.0,
                           "T_wall" => 1.0),
     "evolve_moments" => OptionsDict("density" => false,
diff --git a/moment_kinetics/src/collision_frequencies.jl b/moment_kinetics/src/collision_frequencies.jl
new file mode 100644
index 000000000..3da8aa5b4
--- /dev/null
+++ b/moment_kinetics/src/collision_frequencies.jl
@@ -0,0 +1,118 @@
+module collision_frequencies
+
+export get_collision_frequency_ii, get_collision_frequency_ee, get_collision_frequency_ei
+
+using ..reference_parameters: get_reference_collision_frequency_ii,
+                              get_reference_collision_frequency_ee,
+                              get_reference_collision_frequency_ei
+using ..reference_parameters: setup_reference_parameters
+
+"""
+    get_collision_frequency_ii(collisions, n, vth)
+
+Calculate the ion-ion collision frequency, depending on the settings/parameters in
+`collisions`, for the given density `n` and thermal speed `vth`.
+
+`n` and `vth` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ii(collisions, n, vth)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuii0 = colk.nuii0
+    frequency_option = colk.frequency_option
+    if frequency_option ∈ ("reference_parameters", "collisionality_scan")
+        return @. nuii0 * n * vth^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuii0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+"""
+    get_collision_frequency_ee(collisions, n, vthe)
+
+Calculate the electron-electron collision frequency, depending on the settings/parameters
+in `collisions`, for the given density `n` and electron thermal speed `vthe`.
+
+`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ee(collisions, n, vthe)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuee0 = colk.nuee0
+    frequency_option = colk.frequency_option
+    if frequency_option == "reference_parameters"
+        return @. nuee0 * n * vthe^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuee0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+"""
+    get_collision_frequency_ei(collisions, n, vthe)
+
+Calculate the electron-electron collision frequency, depending on the settings/parameters
+in `collisions`, for the given density `n` and electron thermal speed `vthe`.
+
+`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ei(collisions, n, vthe)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuei0 = colk.nuei0
+    frequency_option = colk.frequency_option
+    if frequency_option == "reference_parameters"
+        return @. nuei0 * n * vthe^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuei0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+end
\ No newline at end of file
diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl
index 61380097b..629695ef7 100644
--- a/moment_kinetics/src/electron_kinetic_equation.jl
+++ b/moment_kinetics/src/electron_kinetic_equation.jl
@@ -38,8 +38,9 @@ using ..em_fields: update_phi!
 using ..external_sources: total_external_electron_sources!,
                           add_total_external_electron_source_to_Jacobian!
 using ..file_io: get_electron_io_info, write_electron_state, finish_electron_io
-using ..krook_collisions: electron_krook_collisions!, get_collision_frequency_ee,
-                          get_collision_frequency_ei,
+using ..collision_frequencies: get_collision_frequency_ee,
+                                 get_collision_frequency_ei
+using ..krook_collisions: electron_krook_collisions!,
                           add_electron_krook_collisions_to_Jacobian!
 using ..timer_utils
 using ..moment_constraints: hard_force_moment_constraints!,
diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl
index 462d23959..7b3544f08 100644
--- a/moment_kinetics/src/em_fields.jl
+++ b/moment_kinetics/src/em_fields.jl
@@ -158,7 +158,7 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
     end
 
     # get gyroaveraged field arrays for distribution function advance
-    gkions = composition.gyrokinetic_ions
+    gkions = composition.ion_physics == gyrokinetic_ions
     if gkions
         gyroaverage_field!(fields.gphi,fields.phi,gyroavs,vperp,z,r,composition)
         gyroaverage_field!(fields.gEz,fields.Ez,gyroavs,vperp,z,r,composition)
diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl
index 6fbd35a96..de20e7a44 100644
--- a/moment_kinetics/src/external_sources.jl
+++ b/moment_kinetics/src/external_sources.jl
@@ -74,6 +74,9 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
                      PI_density_target_z_profile="constant",
                      PI_density_target_z_width=1.0,
                      PI_density_target_z_relative_minimum=0.0,
+                     PI_temperature_controller_P=0.0,
+                     PI_temperature_controller_I=0.0,
+                     PI_temperature_target_amplitude=1.0,
                      recycling_controller_fraction=0.0,
                     )
 
@@ -104,6 +107,10 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
             PI_density_target_ir = nothing
             PI_density_target_iz = nothing
             PI_density_target_rank = nothing
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
         elseif input["source_type"] == "density_midpoint_control"
             PI_density_target = input["PI_density_target_amplitude"]
 
@@ -143,6 +150,53 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
             else
                 PI_density_target_rank = nothing
             end
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
+        elseif input["source_type"] == "temperature_midpoint_control"
+            PI_temperature_target = input["PI_temperature_target_amplitude"]
+            PI_density_target = nothing
+            PI_density_target_ir = nothing
+            PI_density_target_iz = nothing
+            PI_density_target_rank = nothing
+
+            if comm_block[] != MPI.COMM_NULL
+                PI_controller_amplitude = allocate_shared_float(1)
+                controller_source_profile = allocate_shared_float(z.n, r.n)
+            else
+                PI_controller_amplitude = allocate_float(1)
+                controller_source_profile = allocate_float(z.n, r.n)
+            end
+            for ir ∈ 1:r.n, iz ∈ 1:z.n
+                controller_source_profile[iz,ir] = r_amplitude[ir] * z_amplitude[iz]
+            end
+
+            # Find the indices, and process rank of the point at r=0, z=0.
+            # The result of findfirst() will be `nothing` if the point was not found.
+            PI_temperature_target_ir = findfirst(x->abs(x)<1.e-14, r.grid)
+            PI_temperature_target_iz = findfirst(x->abs(x)<1.e-14, z.grid)
+            if block_rank[] == 0
+                # Only need to do communications from the root process of each
+                # shared-memory block
+                if PI_temperature_target_ir !== nothing && PI_temperature_target_iz !== nothing
+                    PI_temperature_target_rank = iblock_index[]
+                else
+                    PI_temperature_target_rank = 0
+                end
+                if comm_inter_block[] != MPI.COMM_NULL
+                    PI_temperature_target_rank = MPI.Allreduce(PI_temperature_target_rank, +,
+                                                           comm_inter_block[])
+                end
+                if PI_temperature_target_rank == 0 && iblock_index[] == 0 &&
+                        (PI_temperature_target_ir === nothing ||
+                         PI_temperature_target_iz === nothing)
+                    error("No grid point with r=0 and z=0 was found for the "
+                          * "'temperature_midpoint' controller.")
+                end
+            else
+                PI_temperature_target_rank = nothing
+            end
         elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses")
             PI_density_target = nothing
             PI_controller_amplitude = nothing
@@ -150,16 +204,22 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
             PI_density_target_ir = nothing
             PI_density_target_iz = nothing
             PI_density_target_rank = nothing
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
         else
             error("Unrecognised ion source_type=$(input["source_type"])."
                   * "Possible values are: Maxwellian, density_profile_control, "
-                  * "density_midpoint_control, energy, alphas, alphas-with-losses, "
-                  * "beam, beam-with-losses")
+                  * "density_midpoint_control, temperature_midpoint_control, energy, "
+                  * "alphas, alphas-with-losses, beam, beam-with-losses")
         end
         return ion_source_data(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude,
                 z_amplitude=z_amplitude, PI_density_target=PI_density_target,
                 PI_controller_amplitude, controller_source_profile,
-                PI_density_target_ir, PI_density_target_iz, PI_density_target_rank)
+                PI_density_target_ir, PI_density_target_iz, PI_density_target_rank,
+                PI_temperature_target, PI_temperature_target_ir, PI_temperature_target_iz,
+                PI_temperature_target_rank)
     end
 
     function get_settings_neutrals(source_index, active_flag)
@@ -328,7 +388,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
                      source_strength=ion_settings.source_strength,
                      source_T=ion_settings.source_T,
                     )
-        if ion_settings.source_type != "energy"
+        if ion_settings.source_type != "energy" && ion_settings.source_type != "temperature_midpoint_control"
             # Need to keep same amplitude for ions and electrons so there is no charge
             # source.
             if input["source_strength"] != ion_settings.source_strength
@@ -413,6 +473,9 @@ function get_source_profile(profile_type, width, relative_minimum, coord)
         L = coord.L
         return @. (1.0 - relative_minimum) * exp(-(x+0.5*L) / width) + relative_minimum +
                   (1.0 - relative_minimum) * exp(-(0.5*L-x) / width) + relative_minimum
+    elseif profile_type == "super_gaussian_4"
+        x = coord.grid
+        return @. (1.0 - relative_minimum) * exp(-(x / width)^4) + relative_minimum
     else
         error("Unrecognised source profile type '$profile_type'.")
     end
@@ -447,7 +510,7 @@ function initialize_external_source_amplitude!(moments, external_source_settings
 
     for index ∈ eachindex(ion_source_settings)
         if ion_source_settings[index].active
-            if ion_source_settings[index].source_type == "energy"
+            if ion_source_settings[index].source_type ∈ ("energy", "temperature_midpoint_control")
                 @loop_r_z ir iz begin
                     moments.ion.external_source_amplitude[iz,ir,index] =
                         ion_source_settings[index].source_strength *
@@ -471,7 +534,7 @@ function initialize_external_source_amplitude!(moments, external_source_settings
                 if moments.evolve_ppar
                     @loop_r_z ir iz begin
                         moments.ion.external_source_pressure_amplitude[iz,ir,index] =
-                            (0.5 * ion_source.source_T +
+                            (0.5 * ion_source_settings[index].source_T +
                             moments.ion.upar[iz,ir]^2 - moments.ion.ppar[iz,ir]) *
                             ion_source_settings[index].source_strength *
                             ion_source_settings[index].r_amplitude[ir] *
@@ -665,7 +728,7 @@ function initialize_external_source_controller_integral!(
         for index ∈ eachindex(ion_source_settings)
             if ion_source_settings[index].active && 
                ion_source_settings[index].PI_density_controller_I != 0.0 && 
-               ion_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control")
+               ion_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control", "temperature_midpoint_control")
                 moments.ion.external_source_controller_integral[:, :, index] .= 0.0
             end
         end
@@ -721,7 +784,7 @@ Add external source term to the ion kinetic equation.
     end
     vpa_grid = vpa.grid
     vperp_grid = vperp.grid
-    if source_type in ("Maxwellian","energy")
+    if source_type in ("Maxwellian","energy","density_midpoint_control","density_profile_control","temperature_midpoint_control")
         begin_s_r_z_vperp_region()
         if moments.evolve_ppar && moments.evolve_upar && moments.evolve_density
             vth = moments.ion.vth
@@ -785,7 +848,7 @@ Add external source term to the ion kinetic equation.
                   * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)")
         end
 
-        if source_type == "energy"
+        if source_type ∈ ("energy", "temperature_midpoint_control")
             if moments.evolve_density
                 # Take particles out of pdf so source does not change density
                 @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
@@ -1254,7 +1317,7 @@ source amplitude.
                     dt * ion_source_settings.PI_density_controller_I * n_error
 
                 # Only want a source, so never allow amplitude to be negative
-                amplitude = max(
+                amplitude = max(ion_source_settings.source_strength +
                     ion_source_settings.PI_density_controller_P * n_error +
                     ion_moments.external_source_controller_integral[1,1,index],
                     0)
@@ -1286,6 +1349,65 @@ source amplitude.
                     amplitude * ion_source_settings.controller_source_profile[iz,ir]
             end
         end
+    elseif ion_source_settings.source_type == "temperature_midpoint_control"
+        begin_serial_region()
+        ion_moments.temp .= 2 .* ppar ./ density
+        # controller_amplitude error is a shared memory Vector of length 1
+        controller_amplitude = ion_source_settings.PI_controller_amplitude
+        @serial_region begin
+            if ion_source_settings.PI_temperature_target_ir !== nothing &&
+                    ion_source_settings.PI_temperature_target_iz !== nothing
+                # This process has the target point
+
+                T_mid = ion_moments.temp[ion_source_settings.PI_temperature_target_iz,
+                                ion_source_settings.PI_temperature_target_ir, is]
+                T_error = ion_source_settings.PI_temperature_target - T_mid
+
+                ion_moments.external_source_controller_integral[1,1,index] +=
+                    dt * ion_source_settings.PI_temperature_controller_I * T_error
+
+                # Only want a source, so never allow amplitude to be negative
+                amplitude = max(ion_source_settings.source_strength +
+                    ion_source_settings.PI_temperature_controller_P * T_error +
+                    ion_moments.external_source_controller_integral[1,1,index],
+                    0)
+            else
+                amplitude = nothing
+            end
+            controller_amplitude[1] =
+                MPI.Bcast(amplitude, ion_source_settings.PI_temperature_target_rank,
+                          comm_inter_block[])
+        end
+
+        begin_r_z_region()
+
+        amplitude = controller_amplitude[1]
+        @loop_r_z ir iz begin
+            ion_moments.external_source_amplitude[iz,ir,index] =
+                amplitude * ion_source_settings.controller_source_profile[iz,ir]
+        end
+
+        if moments.evolve_upar
+            @loop_r_z ir iz begin
+                ion_moments.external_source_momentum_amplitude[iz,ir,index] =
+                    - density[iz,ir] * upar[iz,ir] * amplitude *
+                    ion_source_settings.controller_source_profile[iz,ir]
+            end
+        end
+        if moments.evolve_ppar
+            @loop_r_z ir iz begin
+                ion_moments.external_source_pressure_amplitude[iz,ir,index] =
+                    ((0.5 * ion_source_settings.source_T + 2 * upar[iz,ir]^2) *
+                    amplitude) * ion_source_settings.controller_source_profile[iz,ir]
+            end
+        end
+        #if moments.evolve_ppar
+        #    @loop_r_z ir iz begin
+        #        ion_moments.external_source_pressure_amplitude[iz,ir,index] =
+        #            (0.5 * ion_source_settings.source_T + upar[iz,ir]^2 - ppar[iz,ir]) *
+        #            amplitude * ion_source_settings.controller_source_profile[iz,ir]
+        #    end
+        #end
     elseif ion_source_settings.source_type == "density_profile_control"
         begin_r_z_region()
 
diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl
index 227d6eb66..81853dbdc 100644
--- a/moment_kinetics/src/file_io.jl
+++ b/moment_kinetics/src/file_io.jl
@@ -1482,8 +1482,10 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate,
                     parallel_io=parallel_io,
                     description="Integral term for the PID controller of the external source for ions")
             else
+                r_midpoint = (name="midpoint_controller_r", n=1)
+                z_midpoint = (name="midpoint_controller_z", n=1)
                 external_source_controller_integral = create_dynamic_variable!(
-                    dynamic, "external_source_controller_integral", mk_float;
+                    dynamic, "external_source_controller_integral", mk_float, r_midpoint, z_midpoint, n_sources;
                     parallel_io=parallel_io,
                     description="Integral term for the PID controller of the external source for ions")
             end
@@ -3163,10 +3165,11 @@ function write_ion_moments_data_to_binary(scratch, moments, n_ion_species, t_par
             end
         end
         if io_moments.external_source_controller_integral !== nothing
+            n_sources = size(moments.ion.external_source_amplitude)[3]
             if size(moments.ion.external_source_controller_integral) == (1,1, n_sources)
                 append_to_dynamic_var(io_moments.external_source_controller_integral,
                                       moments.ion.external_source_controller_integral,
-                                      t_idx, parallel_io)
+                                      t_idx, parallel_io, 1, 1, n_sources)
             else
                 append_to_dynamic_var(io_moments.external_source_controller_integral,
                                       moments.ion.external_source_controller_integral,
@@ -3375,14 +3378,15 @@ function write_neutral_moments_data_to_binary(scratch, moments, n_neutral_specie
             end
         end
         if io_moments.external_source_neutral_controller_integral !== nothing
-            if size(moments.neutral.external_source_neutral_controller_integral) == (1,1)
+            n_sources = size(moments.neutral.external_source_amplitude)[3]
+            if size(moments.neutral.external_source_neutral_controller_integral) == (1,1, n_sources)
                 append_to_dynamic_var(io_moments.external_source_neutral_controller_integral,
                                       moments.neutral.external_source_controller_integral,
-                                      t_idx, parallel_io)
+                                      t_idx, parallel_io, 1, 1, n_sources)
             else
                 append_to_dynamic_var(io_moments.external_source_neutral_controller_integral,
                                       moments.neutral.external_source_controller_integral,
-                                      t_idx, parallel_io, z, r)
+                                      t_idx, parallel_io, z, r, n_sources)
             end
         end
         if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar
diff --git a/moment_kinetics/src/gyroaverages.jl b/moment_kinetics/src/gyroaverages.jl
index df649fa31..19a6ddad9 100644
--- a/moment_kinetics/src/gyroaverages.jl
+++ b/moment_kinetics/src/gyroaverages.jl
@@ -13,6 +13,7 @@ using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_shared_float
 using ..array_allocation: allocate_int, allocate_shared_int
 using ..lagrange_polynomials: lagrange_poly
+using ..input_structs: gyrokinetic_ions
 using ..looping
 using ..timer_utils
 using ..communication: MPISharedArray, comm_block, _block_synchronize
@@ -40,7 +41,7 @@ other nonzero boundary conditions for the z and r domains.
 """
 
 function init_gyro_operators(vperp,z,r,gyrophase,geometry,composition;print_info=false)
-    gkions = composition.gyrokinetic_ions
+    gkions = composition.ion_physics == gyrokinetic_ions
     if !gkions
         gyromatrix =  allocate_shared_float(1,1,1,1,1,1)
         gyroloopsizes = allocate_shared_int(1,1,1,1)
diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl
index 2ba8ccd4f..a99ccb53e 100644
--- a/moment_kinetics/src/initial_conditions.jl
+++ b/moment_kinetics/src/initial_conditions.jl
@@ -33,7 +33,7 @@ using ..nonlinear_solvers: nl_solver_info
 using ..velocity_moments: integrate_over_vspace, integrate_over_neutral_vspace
 using ..velocity_moments: integrate_over_positive_vz, integrate_over_negative_vz
 using ..velocity_moments: create_moments_ion, create_moments_electron, create_moments_neutral
-using ..velocity_moments: update_qpar!
+using ..velocity_moments: update_ion_qpar!
 using ..velocity_moments: update_neutral_density!, update_neutral_pz!, update_neutral_pr!, update_neutral_pzeta!
 using ..velocity_moments: update_neutral_uz!, update_neutral_ur!, update_neutral_uzeta!, update_neutral_qz!
 using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, update_vth!, reset_moments_status!
@@ -142,7 +142,7 @@ function init_pdf_and_moments!(pdf, moments, fields, boundary_distributions, geo
                                              r, composition.n_ion_species,
                                              composition.n_neutral_species,
                                              geometry.input, composition, species,
-                                             manufactured_solns_input)
+                                             manufactured_solns_input, collisions)
     else
         n_ion_species = composition.n_ion_species
         n_neutral_species = composition.n_neutral_species
@@ -201,10 +201,11 @@ function init_pdf_and_moments!(pdf, moments, fields, boundary_distributions, geo
                         vpa, vzeta, vr, vz, vpa_spectral, vz_spectral, species)
 
         begin_s_r_z_region()
-        # calculate the initial parallel heat flux from the initial un-normalised pdf
-        update_qpar!(moments.ion.qpar, moments.ion.qpar_updated,
-                     moments.ion.dens, moments.ion.upar, moments.ion.vth,
-                     pdf.ion.norm, vpa, vperp, z, r, composition,
+        # calculate the initial parallel heat flux from the initial un-normalised pdf. Even if coll_krook fluid is being
+        # advanced, initialised ion_qpar uses the pdf 
+        update_ion_qpar!(moments.ion.qpar, moments.ion.qpar_updated,
+                     moments.ion.dens, moments.ion.upar, moments.ion.vth, moments.ion.dT_dz,
+                     pdf.ion.norm, vpa, vperp, z, r, composition, drift_kinetic_ions, collisions,
                      moments.evolve_density, moments.evolve_upar, moments.evolve_ppar)
 
         begin_serial_region()
@@ -1652,7 +1653,8 @@ function init_electron_pdf_over_density_and_boundary_phi!(pdf, phi, density, upa
     return nothing
 end
 
-function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, vperp, z, r, n_ion_species, n_neutral_species, geometry,composition)
+function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, vperp, z, r, n_ion_species, n_neutral_species, 
+                                              geometry, composition, species, manufactured_solns_input, collisions)
     manufactured_solns_list = manufactured_solutions(r.L,z.L,r.bc,z.bc,geometry,composition,r.n)
     dfni_func = manufactured_solns_list.dfni_func
     densi_func = manufactured_solns_list.densi_func
@@ -1679,10 +1681,10 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
                  vpa, vperp, z, r, composition, moments.evolve_density,
                  moments.evolve_upar)
     update_pperp!(moments.ion.pperp, pdf.ion.norm, vpa, vperp, z, r, composition)
-    update_qpar!(moments.ion.qpar, moments.ion.qpar_updated,
+    update_ion_qpar!(moments.ion.qpar, moments.ion.qpar_updated,
                  moments.ion.dens, moments.ion.upar,
-                 moments.ion.vth, pdf.ion.norm, vpa, vperp, z, r,
-                 composition, moments.evolve_density, moments.evolve_upar,
+                 moments.ion.vth, moments.ion.dT_dz, pdf.ion.norm, vpa, vperp, z, r,
+                 composition, drift_kinetic_ions, collisions, moments.evolve_density, moments.evolve_upar,
                  moments.evolve_ppar)
     update_vth!(moments.ion.vth, moments.ion.ppar, moments.ion.pperp, moments.ion.dens, vperp, z, r, composition)
 
diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl
index f427c3a83..13d6de0cf 100644
--- a/moment_kinetics/src/input_structs.jl
+++ b/moment_kinetics/src/input_structs.jl
@@ -165,6 +165,15 @@ export braginskii_fluid
 export kinetic_electrons
 export kinetic_electrons_with_temperature_equation
 
+@enum ion_physics_type begin
+    gyrokinetic_ions
+    drift_kinetic_ions
+    coll_krook_ions
+end
+export ion_physics_type
+export gyrokinetic_ions
+export drift_kinetic_ions
+export coll_krook_ions
 """
 """
 Base.@kwdef struct spatial_initial_condition_input
@@ -264,6 +273,11 @@ Base.@kwdef struct species_composition
     #   density is fixed to be Nₑ*(eϕ/T_e) and N_e is calculated using a current
     #   condition at the wall
     electron_physics::electron_physics_type
+    # ion physics can be drift_kinetic_ions, gyrokinetic_ions and coll_krook_ions
+    # gyrokinetic_ions (originally gyrokinetic_ions = true) -> use gyroaveraged fields at fixed guiding 
+    # centre and moments of the pdf computed at fixed r
+    # drift_kinetic_ions (originally gyrokinetic_ions = false) -> use drift kinetic approximation
+    ion_physics::ion_physics_type
     # if false -- wall bc uses true Knudsen cosine to specify neutral pdf leaving the wall
     # if true -- use a simpler pdf that is easier to integrate
     use_test_neutral_wall_pdf::Bool
@@ -280,10 +294,6 @@ Base.@kwdef struct species_composition
     # The ion flux reaching the wall that is recycled as neutrals is reduced by
     # `recycling_fraction` to account for ions absorbed by the wall.
     recycling_fraction::mk_float
-    # gyrokinetic_ions is a flag determining if the ion species is gyrokinetic
-    # gyrokinetic_ions = true -> use gyroaveraged fields at fixed guiding centre and moments of the pdf computed at fixed r
-    # gyrokinetic_ions = false -> use drift kinetic approximation
-    gyrokinetic_ions::Bool
     # array of structs of parameters for each ion species
     ion::Vector{ion_species_parameters}
     # array of structs of parameters for each neutral species
@@ -352,6 +362,9 @@ Base.@kwdef struct ion_source_data
     PI_density_target_z_profile::String
     PI_density_target_z_width::mk_float
     PI_density_target_z_relative_minimum::mk_float
+    PI_temperature_controller_P::mk_float
+    PI_temperature_controller_I::mk_float
+    PI_temperature_target_amplitude::mk_float
     recycling_controller_fraction::mk_float
     # r_amplitude through the r coordinate (in 1D this can just be set to 1.0)
     r_amplitude::Vector{mk_float}
@@ -359,12 +372,16 @@ Base.@kwdef struct ion_source_data
     # constant profile, parabolic, etc..
     z_amplitude::Vector{mk_float}
     PI_density_target::Union{mk_float, Nothing, MPISharedArray{mk_float,2}}
+    PI_temperature_target::Union{mk_float, Nothing, MPISharedArray{mk_float,2}}
     PI_controller_amplitude::Union{Nothing, MPISharedArray{mk_float,1}}
     controller_source_profile::Union{Nothing, MPISharedArray{mk_float,2}, Array{mk_float, 2}}
     PI_density_target_ir::Union{mk_int, Nothing}
     PI_density_target_iz::Union{mk_int, Nothing}
     PI_density_target_rank::Union{mk_int, Nothing} #possibly this should have Int64 as well, 
     # in the event that the code is running with mk_int = Int32 but the rank is set to 0::Int64
+    PI_temperature_target_ir::Union{mk_int, Nothing}
+    PI_temperature_target_iz::Union{mk_int, Nothing}
+    PI_temperature_target_rank::Union{mk_int, Nothing}
 end
 
 Base.@kwdef struct electron_source_data
diff --git a/moment_kinetics/src/krook_collisions.jl b/moment_kinetics/src/krook_collisions.jl
index 6e37d6dfc..b0339dce9 100644
--- a/moment_kinetics/src/krook_collisions.jl
+++ b/moment_kinetics/src/krook_collisions.jl
@@ -2,19 +2,15 @@
 """
 module krook_collisions
 
-export setup_krook_collisions_input, get_collision_frequency_ii, get_collision_frequency_ee,
-       get_collision_frequency_ei, krook_collisions!, electron_krook_collisions!,
+export setup_krook_collisions_input, krook_collisions!, electron_krook_collisions!,
        add_electron_krook_collisions_to_Jacobian!
 
 using ..looping
 using ..boundary_conditions: skip_f_electron_bc_points_in_Jacobian
 using ..input_structs: krook_collisions_input, set_defaults_and_check_section!
-using ..timer_utils
-using ..reference_parameters: get_reference_collision_frequency_ii,
-                              get_reference_collision_frequency_ee,
-                              get_reference_collision_frequency_ei
-using ..reference_parameters: setup_reference_parameters
-
+using ..timer_utils 
+using ..collision_frequencies
+using ..reference_parameters
 """
 Function for reading Krook collision operator input parameters. 
 Structure the namelist as follows.
@@ -68,113 +64,6 @@ function setup_krook_collisions_input(toml_input::Dict)
     return krook_collisions_input(; input...)
 end
 
-"""
-    get_collision_frequency_ii(collisions, n, vth)
-
-Calculate the ion-ion collision frequency, depending on the settings/parameters in
-`collisions`, for the given density `n` and thermal speed `vth`.
-
-`n` and `vth` may be scalars or arrays, but should have shapes that can be broadcasted
-together.
-"""
-function get_collision_frequency_ii(collisions, n, vth)
-    # extract krook options from collisions struct
-    colk = collisions.krook
-    nuii0 = colk.nuii0
-    frequency_option = colk.frequency_option
-    if frequency_option ∈ ("reference_parameters", "collisionality_scan")
-        return @. nuii0 * n * vth^(-3)
-    elseif frequency_option == "manual"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. nuii0 + 0.0 * n
-    elseif frequency_option == "none"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. 0.0 * n
-    else
-        error("Unrecognised option [krook_collisions] "
-              * "frequency_option=$(frequency_option)")
-    end
-end
-
-"""
-    get_collision_frequency_ee(collisions, n, vthe)
-
-Calculate the electron-electron collision frequency, depending on the settings/parameters
-in `collisions`, for the given density `n` and electron thermal speed `vthe`.
-
-`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
-together.
-"""
-function get_collision_frequency_ee(collisions, n, vthe)
-    # extract krook options from collisions struct
-    colk = collisions.krook
-    nuee0 = colk.nuee0
-    frequency_option = colk.frequency_option
-    if frequency_option == "reference_parameters"
-        return @. nuee0 * n * vthe^(-3)
-    elseif frequency_option == "manual"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. nuee0 + 0.0 * n
-    elseif frequency_option == "none"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. 0.0 * n
-    else
-        error("Unrecognised option [krook_collisions] "
-              * "frequency_option=$(frequency_option)")
-    end
-end
-
-"""
-    get_collision_frequency_ei(collisions, n, vthe)
-
-Calculate the electron-electron collision frequency, depending on the settings/parameters
-in `collisions`, for the given density `n` and electron thermal speed `vthe`.
-
-`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
-together.
-"""
-function get_collision_frequency_ei(collisions, n, vthe)
-    # extract krook options from collisions struct
-    colk = collisions.krook
-    nuei0 = colk.nuei0
-    frequency_option = colk.frequency_option
-    if frequency_option == "reference_parameters"
-        return @. nuei0 * n * vthe^(-3)
-    elseif frequency_option == "manual"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. nuei0 + 0.0 * n
-    elseif frequency_option == "none"
-        # Include 0.0*n so that the result gets promoted to an array if n is an array,
-        # which hopefully means this function will have a fixed return type given the
-        # types of the arguments (we don't want to be 'type unstable' for array inputs by
-        # returning a scalar from this branch but an array from the "reference_parameters"
-        # branch).
-        return @. 0.0 * n
-    else
-        error("Unrecognised option [krook_collisions] "
-              * "frequency_option=$(frequency_option)")
-    end
-end
 
 """
 Add collision operator
diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl
index 53ceda028..4db3fd743 100644
--- a/moment_kinetics/src/load_data.jl
+++ b/moment_kinetics/src/load_data.jl
@@ -26,6 +26,8 @@ using ..file_io: check_io_implementation, get_group, get_subgroup_keys, get_vari
 using ..input_structs
 using ..interpolation: interpolate_to_grid_1d!
 using ..krook_collisions
+using ..collision_frequencies: get_collision_frequency_ii, get_collision_frequency_ee,
+                                get_collision_frequency_ei
 using ..looping
 using ..moment_kinetics_input: mk_input
 using ..neutral_vz_advection: update_speed_neutral_vz!
@@ -775,7 +777,8 @@ function reload_evolving_fields!(pdf, moments, fields, boundary_distributions,
                 if length(moments.ion.external_source_controller_integral) == 1
                     moments.ion.external_source_controller_integral .=
                         load_slice(dynamic, "external_source_controller_integral", time_index)
-                else
+            elseif size(moments.ion.external_source_controller_integral)[1] > 1 ||
+                    size(moments.ion.external_source_controller_integral)[2] > 1 
                     moments.ion.external_source_controller_integral .=
                         reload_moment("external_source_controller_integral", dynamic,
                                       time_index, r, z, r_range, z_range, restart_r,
@@ -4211,6 +4214,88 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
     if variable_name == "temperature"
         vth = get_variable(run_info, "thermal_speed"; kwargs...)
         variable = vth.^2
+    elseif variable_name == "dT_dz"
+        T = get_variable(run_info, "temperature"; kwargs...)
+        variable = similar(T)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], T[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], T[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "dn_dz"
+        n = get_variable(run_info, "density"; kwargs...)
+        variable = similar(n)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], n[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], n[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "dupar_dz"
+        upar = get_variable(run_info, "parallel_flow"; kwargs...)
+        variable = similar(upar)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], upar[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], upar[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "mfp"
+        vth = get_variable(run_info, "thermal_speed"; kwargs...)
+        nu_ii = get_variable(run_info, "collision_frequency_ii"; kwargs...)
+        variable = vth ./ nu_ii
+    elseif variable_name == "L_T"
+        dT_dz = get_variable(run_info, "dT_dz"; kwargs...)
+        temp = get_variable(run_info, "temperature"; kwargs...)
+        # We define gradient lengthscale of T as LT^-1 = dln(T)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(temp .* dT_dz.^(-1))
+        # flat points in temperature have diverging LT, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
+    elseif variable_name == "L_n"
+        dn_dz = get_variable(run_info, "dn_dz"; kwargs...)
+        n = get_variable(run_info, "density"; kwargs...)
+        # We define gradient lengthscale of n as Ln^-1 = dln(n)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(n .* dn_dz.^(-1))
+        # flat points in temperature have diverging Ln, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
+    elseif variable_name == "L_upar"
+        dupar_dz = get_variable(run_info, "dupar_dz"; kwargs...)
+        upar = get_variable(run_info, "parallel_flow"; kwargs...)
+        # We define gradient lengthscale of upar as Lupar^-1 = dln(upar)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(upar .* dupar_dz.^(-1))
+        # flat points in temperature have diverging Lupar, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
+    elseif variable_name == "coll_krook_heat_flux"
+        n = get_variable(run_info, "density"; kwargs...)
+        vth = get_variable(run_info, "thermal_speed"; kwargs...)
+        dT_dz = get_variable(run_info, "dT_dz"; kwargs...)
+        nu_ii = get_variable(run_info, "collision_frequency_ii"; kwargs...)
+        variable = @. -(1/2) * 5/4 * n * vth^2 * dT_dz / nu_ii
     elseif variable_name == "collision_frequency_ii"
         n = get_variable(run_info, "density"; kwargs...)
         vth = get_variable(run_info, "thermal_speed"; kwargs...)
@@ -4322,7 +4407,7 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
             n = get_variable(run_info, "density"; kwargs...)
 
             # Note factor of 0.5 in front of qpar because the definition of qpar (see e.g.
-            # `update_qpar_species!()`) is unconventional (i.e. missing a factor of 0.5).
+            # `update_ion_qpar_species!()`) is unconventional (i.e. missing a factor of 0.5).
             # Factor of 3/2 in front of 1/2*n*vth^2*upar because this in 1V - would be 5/2
             # for 2V/3V cases.
             variable = @. 0.5*qpar + 0.75*n*vth^2*upar + 0.5*n*upar^3
@@ -4337,7 +4422,7 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
             n = get_variable(run_info, "density_neutral"; kwargs...)
 
             # Note factor of 0.5 in front of qpar because the definition of qpar (see e.g.
-            # `update_qpar_species!()`) is unconventional (i.e. missing a factor of 0.5).
+            # `update_ion_qpar_species!()`) is unconventional (i.e. missing a factor of 0.5).
             # Factor of 3/2 in front of 1/2*n*vth^2*upar because this in 1V - would be 5/2
             # for 2V/3V cases.
             variable = @. 0.5*qpar + 0.75*n*vth^2*upar + 0.5*n*upar^3
@@ -4409,9 +4494,6 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
         density = get_variable(run_info, "density")
         upar = get_variable(run_info, "parallel_flow")
         ppar = get_variable(run_info, "parallel_pressure")
-        density_neutral = get_variable(run_info, "density_neutral")
-        uz_neutral = get_variable(run_info, "uz_neutral")
-        pz_neutral = get_variable(run_info, "pz_neutral")
         vth = get_variable(run_info, "thermal_speed")
         dupar_dz = get_z_derivative(run_info, "parallel_flow")
         dppar_dz = get_z_derivative(run_info, "parallel_pressure")
@@ -4461,6 +4543,15 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
         setup_loop_ranges!(0, 1; s=nspecies, sn=run_info.n_neutral_species, r=nr, z=nz,
                            vperp=nvperp, vpa=nvpa, vzeta=run_info.vzeta.n,
                            vr=run_info.vr.n, vz=run_info.vz.n)
+        
+        # Use neutrals for fvec calculation in moment_kinetic version only when 
+        # n_neutrals != 0
+        if run_info.n_neutral_species != 0
+            density_neutral = get_variable(run_info, "density_neutral")
+            uz_neutral = get_variable(run_info, "uz_neutral")
+            pz_neutral = get_variable(run_info, "pz_neutral")
+        end
+
         for it ∈ 1:nt
             begin_serial_region()
             # Only need some struct with a 'speed' variable
@@ -4479,12 +4570,18 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
                              evolve_density=run_info.evolve_density,
                              evolve_upar=run_info.evolve_upar,
                              evolve_ppar=run_info.evolve_ppar)
-            @views fvec = (density=density[:,:,:,it],
-                           upar=upar[:,:,:,it],
-                           ppar=ppar[:,:,:,it],
-                           density_neutral=density_neutral[:,:,:,it],
-                           uz_neutral=uz_neutral[:,:,:,it],
-                           pz_neutral=pz_neutral[:,:,:,it])
+            if run_info.n_neutral_species != 0
+                @views fvec = (density=density[:,:,:,it],
+                            upar=upar[:,:,:,it],
+                            ppar=ppar[:,:,:,it],
+                            density_neutral=density_neutral[:,:,:,it],
+                            uz_neutral=uz_neutral[:,:,:,it],
+                            pz_neutral=pz_neutral[:,:,:,it])
+            else
+                @views fvec = (density=density[:,:,:,it],
+                            upar=upar[:,:,:,it],
+                            ppar=ppar[:,:,:,it])
+            end
             @views update_speed_vpa!(advect, fields, fvec, moments, run_info.vpa,
                                      run_info.vperp, run_info.z, run_info.r,
                                      run_info.composition, run_info.collisions,
diff --git a/moment_kinetics/src/moment_constraints.jl b/moment_kinetics/src/moment_constraints.jl
index ac91dd321..c881f4ca9 100644
--- a/moment_kinetics/src/moment_constraints.jl
+++ b/moment_kinetics/src/moment_constraints.jl
@@ -10,7 +10,7 @@ using ..communication: _block_synchronize
 using ..looping
 using ..timer_utils
 using ..type_definitions: mk_float
-using ..velocity_moments: integrate_over_vspace, update_qpar!
+using ..velocity_moments: integrate_over_vspace, update_ion_qpar!
 
 export hard_force_moment_constraints!, hard_force_moment_constraints_neutral!,
        electron_implicit_constraint_forcing!,
diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl
index 1f2c86593..089d8f130 100644
--- a/moment_kinetics/src/moment_kinetics.jl
+++ b/moment_kinetics/src/moment_kinetics.jl
@@ -38,9 +38,11 @@ include("nonlinear_solvers.jl")
 include("file_io.jl")
 include("geo.jl")
 include("gyroaverages.jl")
+include("collision_frequencies.jl")
 include("velocity_moments.jl")
 include("velocity_grid_transforms.jl")
 include("boundary_conditions.jl")
+include("krook_collisions.jl")
 include("electron_fluid_equations.jl")
 include("em_fields.jl")
 include("bgk.jl")
@@ -62,7 +64,6 @@ include("neutral_z_advection.jl")
 include("neutral_vz_advection.jl")
 include("charge_exchange.jl")
 include("ionization.jl")
-include("krook_collisions.jl")
 include("maxwell_diffusion.jl")
 include("continuity.jl")
 include("energy_equation.jl")
diff --git a/moment_kinetics/src/moment_kinetics_structs.jl b/moment_kinetics/src/moment_kinetics_structs.jl
index 18b5b0189..4df2043db 100644
--- a/moment_kinetics/src/moment_kinetics_structs.jl
+++ b/moment_kinetics/src/moment_kinetics_structs.jl
@@ -95,6 +95,8 @@ struct moments_ion_substruct
     qpar_updated::Vector{Bool}
     # this is the thermal speed based on the parallel temperature Tpar = ppar/dens: vth = sqrt(2*Tpar/m)
     vth::MPISharedArray{mk_float,3}
+    # this is the temperature, calculated from 2ppar/dens (the comment above may be out of date?)
+    temp::MPISharedArray{mk_float,3}
     # generalised Chodura integrals for the lower and upper plates
     chodura_integral_lower::MPISharedArray{mk_float,2}
     chodura_integral_upper::MPISharedArray{mk_float,2}
@@ -125,6 +127,8 @@ struct moments_ion_substruct
     dqpar_dz::Union{MPISharedArray{mk_float,3},Nothing}
     # this is the z-derivative of the thermal speed based on the parallel temperature Tpar = ppar/dens: vth = sqrt(2*Tpar/m)
     dvth_dz::Union{MPISharedArray{mk_float,3},Nothing}
+    # this is the z-derivative of the temperature
+    dT_dz::Union{MPISharedArray{mk_float,3},Nothing}
     # this is the entropy production dS/dt = - int (ln f sum_s' C_ss' [f_s,f_s']) d^3 v
     dSdt::MPISharedArray{mk_float,3}
     # Spatially varying amplitude of the external source term (third index is for different sources)
diff --git a/moment_kinetics/src/species_input.jl b/moment_kinetics/src/species_input.jl
index 59d7043fe..a4b1ef641 100644
--- a/moment_kinetics/src/species_input.jl
+++ b/moment_kinetics/src/species_input.jl
@@ -11,6 +11,7 @@ using ..input_structs: set_defaults_and_check_section!
 using ..input_structs: species_composition, ion_species_parameters, neutral_species_parameters
 using ..input_structs: spatial_initial_condition_input, velocity_initial_condition_input
 using ..input_structs: boltzmann_electron_response, boltzmann_electron_response_with_simple_sheath
+using ..input_structs: drift_kinetic_ions
 using ..reference_parameters: setup_reference_parameters
 
 function get_species_input(toml_input)
@@ -29,6 +30,15 @@ function get_species_input(toml_input)
         #   electron density is fixed to be N_e*(eϕ/T_e) and N_e is calculated w.r.t a
         #   reference value using J_||e + J_||i = 0 at z = 0
         electron_physics = boltzmann_electron_response,
+        # If ion_physics=drift_kinetic_ions, the ion distribution function is advanced in 
+        # time in the drift kinetic approximation like usual. 
+        # If ion_physics=gyrokinetic_ions, the ion distribution function is
+        # advanced in time using gyroaveraged fields at fixed guiding centre and moments of the
+        # pdf computed at fixed r
+        # If ion_physics=coll_krook_ions, there is no need for a shape function to evolve, and the code 
+        # only evolves ions in a fluid sense (i.e. all evolve_moments are set to true), with a 
+        # coll_krook closure for the ion heat flux.
+        ion_physics = drift_kinetic_ions,
         # initial Tₑ = 1
         T_e = 1.0,
         # wall temperature T_wall = Tw/Te
@@ -43,10 +53,7 @@ function get_species_input(toml_input)
         use_test_neutral_wall_pdf = false,
         # The ion flux reaching the wall that is recycled as neutrals is reduced by
         # `recycling_fraction` to account for ions absorbed by the wall.
-        recycling_fraction = 1.0,
-        # gyrokinetic_ions = True -> use gyroaveraged fields at fixed guiding centre and moments of the pdf computed at fixed r
-        # gyrokinetic_ions = False -> use drift kinetic approximation
-        gyrokinetic_ions = false)
+        recycling_fraction = 1.0)
 
     nspec_ion = composition_section["n_ion_species"]
     nspec_neutral = composition_section["n_neutral_species"]
diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl
index b9d85420d..5005ef240 100644
--- a/moment_kinetics/src/time_advance.jl
+++ b/moment_kinetics/src/time_advance.jl
@@ -20,7 +20,7 @@ using ..initial_conditions: initialize_electrons!
 using ..looping
 using ..moment_kinetics_structs: scratch_pdf, scratch_electron_pdf
 using ..velocity_moments: update_moments!, update_moments_neutral!, reset_moments_status!, update_derived_moments!, update_derived_moments_neutral!
-using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_qpar!, update_vth!
+using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_ion_qpar!, update_vth!
 using ..velocity_moments: update_neutral_density!, update_neutral_qz!
 using ..velocity_moments: update_neutral_uzeta!, update_neutral_uz!, update_neutral_ur!
 using ..velocity_moments: update_neutral_pzeta!, update_neutral_pz!, update_neutral_pr!
@@ -948,7 +948,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop
         # constraints to the pdf
         reset_moments_status!(moments)
         update_moments!(moments, pdf.ion.norm, gyroavs, vpa, vperp, z, r, composition,
-           r_spectral,geometry,scratch_dummy,z_advect)
+           r_spectral,geometry,scratch_dummy,z_advect, collisions)
         # enforce boundary conditions in r and z on the neutral particle distribution function
         if n_neutral_species > 0
             # Note, so far vr and vzeta do not need advect objects, so pass `nothing` for
@@ -2354,8 +2354,13 @@ moments and moment derivatives
     # Note these may be needed for the boundary condition on the neutrals, so must be
     # calculated before that is applied. Also may be needed to calculate advection speeds
     # for for CFL stability limit calculations in adaptive_timestep_update!().
-    update_derived_moments!(this_scratch, moments, vpa, vperp, z, r, composition,
-        r_spectral, geometry, gyroavs, scratch_dummy, z_advect, diagnostic_moments)
+    if composition.ion_physics ∈ (drift_kinetic_ions, gyrokinetic_ions)
+        update_derived_moments!(this_scratch, moments, vpa, vperp, z, r, composition,
+            r_spectral, geometry, gyroavs, scratch_dummy, z_advect, collisions, diagnostic_moments)
+    else
+        update_derived_moments!(this_scratch, moments, vpa, vperp, z, r, composition,
+            r_spectral, geometry, gyroavs, scratch_dummy, z_advect, collisions, false)
+    end
 
     calculate_ion_moment_derivatives!(moments, this_scratch, scratch_dummy, z, z_spectral,
                                       num_diss_params.ion.moment_dissipation_coefficient)
@@ -3254,167 +3259,169 @@ with fvec_in an input and fvec_out the output
                                             external_source_settings.neutral, r, z, dt)
     end
 
-    if advance.vpa_advection
-        vpa_advection!(fvec_out.pdf, fvec_in, fields, moments, vpa_advect, vpa, vperp, z, r, dt, t,
-            vpa_spectral, composition, collisions, external_source_settings.ion, geometry)
-    end
-
-    # z_advection! advances 1D advection equation in z
-    # apply z-advection operation to ion species
-
-    if advance.z_advection
-        z_advection!(fvec_out.pdf, fvec_in, moments, fields, z_advect, z, vpa, vperp, r,
-                     dt, t, z_spectral, composition, geometry, scratch_dummy)
-    end
-
-    # r advection relies on derivatives in z to get ExB
-    if advance.r_advection
-        r_advection!(fvec_out.pdf, fvec_in, moments, fields, r_advect, r, z, vperp, vpa,
-                     dt, r_spectral, composition, geometry, scratch_dummy)
-    end
-    # vperp_advection requires information about z and r advection
-    # so call vperp_advection! only after z and r advection routines
-    if advance.vperp_advection
-        vperp_advection!(fvec_out.pdf, fvec_in, vperp_advect, r, z, vperp, vpa,
-                      dt, vperp_spectral, composition, z_advect, r_advect, geometry,
-                      moments, fields, t)
-    end
+    if composition.ion_physics ∈ (drift_kinetic_ions, gyrokinetic_ions)
+        if advance.vpa_advection
+            vpa_advection!(fvec_out.pdf, fvec_in, fields, moments, vpa_advect, vpa, vperp, z, r, dt, t,
+                vpa_spectral, composition, collisions, external_source_settings.ion, geometry)
+        end
 
-    if advance.source_terms
-        source_terms!(fvec_out.pdf, fvec_in, moments, vpa, z, r, dt, z_spectral,
-                      composition, collisions, external_source_settings.ion)
-    end
+        # z_advection! advances 1D advection equation in z
+        # apply z-advection operation to ion species
 
-    if advance.neutral_z_advection
-        neutral_advection_z!(fvec_out.pdf_neutral, fvec_in, moments, neutral_z_advect,
-            r, z, vzeta, vr, vz, dt, t, z_spectral, composition, scratch_dummy)
-    end
+        if advance.z_advection
+            z_advection!(fvec_out.pdf, fvec_in, moments, fields, z_advect, z, vpa, vperp, r,
+                        dt, t, z_spectral, composition, geometry, scratch_dummy)
+        end
 
-    if advance.neutral_r_advection
-        neutral_advection_r!(fvec_out.pdf_neutral, fvec_in, neutral_r_advect,
-            r, z, vzeta, vr, vz, dt, r_spectral, composition, geometry, scratch_dummy)
-    end
+        # r advection relies on derivatives in z to get ExB
+        if advance.r_advection
+            r_advection!(fvec_out.pdf, fvec_in, moments, fields, r_advect, r, z, vperp, vpa,
+                        dt, r_spectral, composition, geometry, scratch_dummy)
+        end
+        # vperp_advection requires information about z and r advection
+        # so call vperp_advection! only after z and r advection routines
+        if advance.vperp_advection
+            vperp_advection!(fvec_out.pdf, fvec_in, vperp_advect, r, z, vperp, vpa,
+                        dt, vperp_spectral, composition, z_advect, r_advect, geometry,
+                        moments, fields, t)
+        end
 
-    if advance.neutral_vz_advection
-        neutral_advection_vz!(fvec_out.pdf_neutral, fvec_in, fields, moments,
-                              neutral_vz_advect, vz, vr, vzeta, z, r, dt, vz_spectral,
-                              composition, collisions, external_source_settings.neutral)
-    end
+        if advance.source_terms
+            source_terms!(fvec_out.pdf, fvec_in, moments, vpa, z, r, dt, z_spectral,
+                        composition, collisions, external_source_settings.ion)
+        end
 
-    if advance.neutral_source_terms
-        source_terms_neutral!(fvec_out.pdf_neutral, fvec_in, moments, vpa, z, r, dt, z_spectral,
-                      composition, collisions, external_source_settings.neutral)
-    end
+        if advance.neutral_z_advection
+            neutral_advection_z!(fvec_out.pdf_neutral, fvec_in, moments, neutral_z_advect,
+                r, z, vzeta, vr, vz, dt, t, z_spectral, composition, scratch_dummy)
+        end
 
-    if advance.manufactured_solns_test
-        source_terms_manufactured!(fvec_out.pdf, fvec_out.pdf_neutral, vz, vr, vzeta, vpa, vperp, z, r, t, dt, composition, manufactured_source_list)
-    end
+        if advance.neutral_r_advection
+            neutral_advection_r!(fvec_out.pdf_neutral, fvec_in, neutral_r_advect,
+                r, z, vzeta, vr, vz, dt, r_spectral, composition, geometry, scratch_dummy)
+        end
 
-    if advance.ion_cx_collisions || advance.ion_ionization_collisions
-        # gyroaverage neutral dfn and place it in the ion.buffer array for use in the collisions step
-        vzvrvzeta_to_vpavperp!(pdf.ion.buffer, fvec_in.pdf_neutral, vz, vr, vzeta, vpa, vperp, gyrophase, z, r, geometry, composition)
-    end
-    if advance.neutral_cx_collisions || advance.neutral_ionization_collisions
-        # interpolate ion particle dfn and place it in the neutral.buffer array for use in the collisions step
-        vpavperp_to_vzvrvzeta!(pdf.neutral.buffer, fvec_in.pdf, vz, vr, vzeta, vpa, vperp, z, r, geometry, composition)
-    end
+        if advance.neutral_vz_advection
+            neutral_advection_vz!(fvec_out.pdf_neutral, fvec_in, fields, moments,
+                                neutral_vz_advect, vz, vr, vzeta, z, r, dt, vz_spectral,
+                                composition, collisions, external_source_settings.neutral)
+        end
 
-    # account for charge exchange collisions between ions and neutrals
-    if advance.ion_cx_collisions_1V
-        ion_charge_exchange_collisions_1V!(fvec_out.pdf, fvec_in, moments, composition,
-                                           vpa, vz, collisions.reactions.charge_exchange_frequency,
-                                           vpa_spectral, vz_spectral, dt)
-    elseif advance.ion_cx_collisions
-        ion_charge_exchange_collisions_3V!(fvec_out.pdf, pdf.ion.buffer, fvec_in,
-                                           composition, vz, vr, vzeta, vpa, vperp, z, r,
-                                           collisions.reactions.charge_exchange_frequency, dt)
-    end
-    if advance.neutral_cx_collisions_1V
-        neutral_charge_exchange_collisions_1V!(fvec_out.pdf_neutral, fvec_in, moments,
-                                               composition, vpa, vz,
-                                               collisions.reactions.charge_exchange_frequency, vpa_spectral,
-                                               vz_spectral, dt)
-    elseif advance.neutral_cx_collisions
-        neutral_charge_exchange_collisions_3V!(fvec_out.pdf_neutral, pdf.neutral.buffer,
-                                               fvec_in, composition, vz, vr, vzeta, vpa,
-                                               vperp, z, r, collisions.reactions.charge_exchange_frequency,
-                                               dt)
-    end
-    # account for ionization collisions between ions and neutrals
-    if advance.ion_ionization_collisions_1V
-        ion_ionization_collisions_1V!(fvec_out.pdf, fvec_in, vz, vpa, vperp, z, r,
-                                      vz_spectral, moments, composition, collisions, dt)
-    elseif advance.ion_ionization_collisions
-        ion_ionization_collisions_3V!(fvec_out.pdf, pdf.ion.buffer, fvec_in, composition,
-                                      vz, vr, vzeta, vpa, vperp, z, r, collisions, dt)
-    end
-    if advance.neutral_ionization_collisions_1V
-        neutral_ionization_collisions_1V!(fvec_out.pdf_neutral, fvec_in, vz, vpa, vperp,
-                                          z, r, vz_spectral, moments, composition,
-                                          collisions, dt)
-    elseif advance.neutral_ionization_collisions
-        neutral_ionization_collisions_3V!(fvec_out.pdf_neutral, fvec_in, composition, vz,
-                                          vr, vzeta, vpa, vperp, z, r, collisions, dt)
-    end
+        if advance.neutral_source_terms
+            source_terms_neutral!(fvec_out.pdf_neutral, fvec_in, moments, vpa, z, r, dt, z_spectral,
+                        composition, collisions, external_source_settings.neutral)
+        end
 
-    # Add Krook collision operator for ions
-    if advance.krook_collisions_ii
-        krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions,
-                          vperp, vpa, dt)
-    end
-    # Add maxwellian diffusion collision operator for ions
-    if advance.mxwl_diff_collisions_ii
-        ion_vpa_maxwell_diffusion!(fvec_out.pdf, fvec_in, moments, vpa, vperp, vpa_spectral, 
-                                   dt, collisions.mxwl_diff.D_ii)
-    end
-    # Add maxwellian diffusion collision operator for neutrals
-    if advance.mxwl_diff_collisions_nn
-        neutral_vz_maxwell_diffusion!(fvec_out.pdf_neutral, fvec_in, moments, vzeta, vr, vz, vz_spectral, 
-                                   dt, collisions.mxwl_diff.D_nn)
-    end
+        if advance.manufactured_solns_test
+            source_terms_manufactured!(fvec_out.pdf, fvec_out.pdf_neutral, vz, vr, vzeta, vpa, vperp, z, r, t, dt, composition, manufactured_source_list)
+        end
 
-    if advance.external_source
-        total_external_ion_sources!(fvec_out.pdf, fvec_in, moments, external_source_settings.ion,
-                            vperp, vpa, dt, scratch_dummy)
-    end
-    if advance.neutral_external_source
-        total_external_neutral_sources!(fvec_out.pdf_neutral, fvec_in, moments,
-                                external_source_settings.neutral, vzeta, vr, vz, dt)
-    end
-
-    # add numerical dissipation
-    if advance.ion_numerical_dissipation
-        vpa_dissipation!(fvec_out.pdf, fvec_in.pdf, vpa, vpa_spectral, dt,
-                         num_diss_params.ion.vpa_dissipation_coefficient)
-        vperp_dissipation!(fvec_out.pdf, fvec_in.pdf, vperp, vperp_spectral, dt,
-                         num_diss_params.ion.vperp_dissipation_coefficient)
-        z_dissipation!(fvec_out.pdf, fvec_in.pdf, z, z_spectral, dt,
-                       num_diss_params.ion.z_dissipation_coefficient, scratch_dummy)
-        r_dissipation!(fvec_out.pdf, fvec_in.pdf, r, r_spectral, dt,
-                       num_diss_params.ion.r_dissipation_coefficient, scratch_dummy)
-    end
-    if advance.neutral_numerical_dissipation
-        vz_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, vz,
-                                vz_spectral, dt, num_diss_params.neutral.vz_dissipation_coefficient)
-        z_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, z, z_spectral,
-                               dt, num_diss_params.neutral.z_dissipation_coefficient, scratch_dummy)
-        r_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, r, r_spectral,
-                               dt, num_diss_params.neutral.r_dissipation_coefficient, scratch_dummy)
-    end
-    # advance with the Fokker-Planck self-collision operator
-    if advance.explicit_weakform_fp_collisions
-        update_entropy_diagnostic = (istage == 1)
-        if collisions.fkpl.self_collisions
-            # self collisions for each species
-            explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,composition,
-                                 collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy,
-                                                     diagnose_entropy_production = update_entropy_diagnostic)
+        if advance.ion_cx_collisions || advance.ion_ionization_collisions
+            # gyroaverage neutral dfn and place it in the ion.buffer array for use in the collisions step
+            vzvrvzeta_to_vpavperp!(pdf.ion.buffer, fvec_in.pdf_neutral, vz, vr, vzeta, vpa, vperp, gyrophase, z, r, geometry, composition)
         end
-        if collisions.fkpl.slowing_down_test
-        # include cross-collsions with fixed Maxwellian backgrounds
-            explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,
-                             composition,collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,
-                                             diagnose_entropy_production = update_entropy_diagnostic)
+        if advance.neutral_cx_collisions || advance.neutral_ionization_collisions
+            # interpolate ion particle dfn and place it in the neutral.buffer array for use in the collisions step
+            vpavperp_to_vzvrvzeta!(pdf.neutral.buffer, fvec_in.pdf, vz, vr, vzeta, vpa, vperp, z, r, geometry, composition)
+        end
+
+        # account for charge exchange collisions between ions and neutrals
+        if advance.ion_cx_collisions_1V
+            ion_charge_exchange_collisions_1V!(fvec_out.pdf, fvec_in, moments, composition,
+                                            vpa, vz, collisions.reactions.charge_exchange_frequency,
+                                            vpa_spectral, vz_spectral, dt)
+        elseif advance.ion_cx_collisions
+            ion_charge_exchange_collisions_3V!(fvec_out.pdf, pdf.ion.buffer, fvec_in,
+                                            composition, vz, vr, vzeta, vpa, vperp, z, r,
+                                            collisions.reactions.charge_exchange_frequency, dt)
+        end
+        if advance.neutral_cx_collisions_1V
+            neutral_charge_exchange_collisions_1V!(fvec_out.pdf_neutral, fvec_in, moments,
+                                                composition, vpa, vz,
+                                                collisions.reactions.charge_exchange_frequency, vpa_spectral,
+                                                vz_spectral, dt)
+        elseif advance.neutral_cx_collisions
+            neutral_charge_exchange_collisions_3V!(fvec_out.pdf_neutral, pdf.neutral.buffer,
+                                                fvec_in, composition, vz, vr, vzeta, vpa,
+                                                vperp, z, r, collisions.reactions.charge_exchange_frequency,
+                                                dt)
+        end
+        # account for ionization collisions between ions and neutrals
+        if advance.ion_ionization_collisions_1V
+            ion_ionization_collisions_1V!(fvec_out.pdf, fvec_in, vz, vpa, vperp, z, r,
+                                        vz_spectral, moments, composition, collisions, dt)
+        elseif advance.ion_ionization_collisions
+            ion_ionization_collisions_3V!(fvec_out.pdf, pdf.ion.buffer, fvec_in, composition,
+                                        vz, vr, vzeta, vpa, vperp, z, r, collisions, dt)
+        end
+        if advance.neutral_ionization_collisions_1V
+            neutral_ionization_collisions_1V!(fvec_out.pdf_neutral, fvec_in, vz, vpa, vperp,
+                                            z, r, vz_spectral, moments, composition,
+                                            collisions, dt)
+        elseif advance.neutral_ionization_collisions
+            neutral_ionization_collisions_3V!(fvec_out.pdf_neutral, fvec_in, composition, vz,
+                                            vr, vzeta, vpa, vperp, z, r, collisions, dt)
+        end
+
+        # Add Krook collision operator for ions
+        if advance.krook_collisions_ii
+            krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions,
+                            vperp, vpa, dt)
+        end
+        # Add maxwellian diffusion collision operator for ions
+        if advance.mxwl_diff_collisions_ii
+            ion_vpa_maxwell_diffusion!(fvec_out.pdf, fvec_in, moments, vpa, vperp, vpa_spectral, 
+                                    dt, collisions.mxwl_diff.D_ii)
+        end
+        # Add maxwellian diffusion collision operator for neutrals
+        if advance.mxwl_diff_collisions_nn
+            neutral_vz_maxwell_diffusion!(fvec_out.pdf_neutral, fvec_in, moments, vzeta, vr, vz, vz_spectral, 
+                                    dt, collisions.mxwl_diff.D_nn)
+        end
+
+        if advance.external_source
+            total_external_ion_sources!(fvec_out.pdf, fvec_in, moments, external_source_settings.ion,
+                                vperp, vpa, dt, scratch_dummy)
+        end
+        if advance.neutral_external_source
+            total_external_neutral_sources!(fvec_out.pdf_neutral, fvec_in, moments,
+                                    external_source_settings.neutral, vzeta, vr, vz, dt)
+        end
+
+        # add numerical dissipation
+        if advance.ion_numerical_dissipation
+            vpa_dissipation!(fvec_out.pdf, fvec_in.pdf, vpa, vpa_spectral, dt,
+                            num_diss_params.ion.vpa_dissipation_coefficient)
+            vperp_dissipation!(fvec_out.pdf, fvec_in.pdf, vperp, vperp_spectral, dt,
+                            num_diss_params.ion.vperp_dissipation_coefficient)
+            z_dissipation!(fvec_out.pdf, fvec_in.pdf, z, z_spectral, dt,
+                        num_diss_params.ion.z_dissipation_coefficient, scratch_dummy)
+            r_dissipation!(fvec_out.pdf, fvec_in.pdf, r, r_spectral, dt,
+                        num_diss_params.ion.r_dissipation_coefficient, scratch_dummy)
+        end
+        if advance.neutral_numerical_dissipation
+            vz_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, vz,
+                                    vz_spectral, dt, num_diss_params.neutral.vz_dissipation_coefficient)
+            z_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, z, z_spectral,
+                                dt, num_diss_params.neutral.z_dissipation_coefficient, scratch_dummy)
+            r_dissipation_neutral!(fvec_out.pdf_neutral, fvec_in.pdf_neutral, r, r_spectral,
+                                dt, num_diss_params.neutral.r_dissipation_coefficient, scratch_dummy)
+        end
+        # advance with the Fokker-Planck self-collision operator
+        if advance.explicit_weakform_fp_collisions
+            update_entropy_diagnostic = (istage == 1)
+            if collisions.fkpl.self_collisions
+                # self collisions for each species
+                explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,composition,
+                                    collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy,
+                                                        diagnose_entropy_production = update_entropy_diagnostic)
+            end
+            if collisions.fkpl.slowing_down_test
+            # include cross-collsions with fixed Maxwellian backgrounds
+                explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,
+                                composition,collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,
+                                                diagnose_entropy_production = update_entropy_diagnostic)
+            end
         end
     end
     
@@ -3760,7 +3767,7 @@ Do a backward-Euler timestep for all terms in the ion kinetic equation.
         # Ensure moments are consistent with f_new
         update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition,
                                 r_spectral, geometry, gyroavs, scratch_dummy, z_advect,
-                                false)
+                                collisions, false)
         calculate_ion_moment_derivatives!(moments, new_scratch, scratch_dummy, z,
                                           z_spectral,
                                           num_diss_params.ion.moment_dissipation_coefficient)
diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl
index ebb1aeaa2..c58be5cbf 100644
--- a/moment_kinetics/src/velocity_moments.jl
+++ b/moment_kinetics/src/velocity_moments.jl
@@ -11,7 +11,7 @@ export update_density!
 export update_upar!
 export update_ppar!
 export update_pperp!
-export update_qpar!
+export update_ion_qpar!
 export update_vth!
 export reset_moments_status!
 export update_neutral_density!
@@ -41,10 +41,12 @@ using ..derivatives: derivative_z!, second_derivative_z!
 using ..derivatives: derivative_r!, second_derivative_r!
 using ..looping
 using ..gyroaverages: gyro_operators, gyroaverage_pdf!
+using ..collision_frequencies: get_collision_frequency_ii
 using ..input_structs
 using ..moment_kinetics_structs: moments_ion_substruct, moments_electron_substruct,
                                  moments_neutral_substruct
 
+
 #global tmpsum1 = 0.0
 #global tmpsum2 = 0.0
 #global dens_hist = zeros(17,1)
@@ -71,6 +73,8 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar,
     parallel_pressure_updated .= false
     # allocate array used for the perpendicular pressure
     perpendicular_pressure = allocate_shared_float(nz, nr, n_species)
+    # allocate array used for the temperature
+    temperature = allocate_shared_float(nz, nr, n_species)
     # allocate array used for the parallel flow
     parallel_heat_flux = allocate_shared_float(nz, nr, n_species)
     # allocate array of Bools that indicate if the parallel flow is updated for each species
@@ -129,11 +133,13 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar,
         d2ppar_dz2 = allocate_shared_float(nz, nr, n_species)
         dqpar_dz = allocate_shared_float(nz, nr, n_species)
         dvth_dz = allocate_shared_float(nz, nr, n_species)
+        dT_dz = allocate_shared_float(nz, nr, n_species)
     else
         dppar_dz_upwind = nothing
         d2ppar_dz2 = nothing
         dqpar_dz = nothing
         dvth_dz = nothing
+        dT_dz = nothing
     end
 
     entropy_production = allocate_shared_float(nz, nr, n_species)
@@ -187,10 +193,10 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar,
     # return struct containing arrays needed to update moments
     return moments_ion_substruct(density, density_updated, parallel_flow,
         parallel_flow_updated, parallel_pressure, parallel_pressure_updated,perpendicular_pressure,
-        parallel_heat_flux, parallel_heat_flux_updated, thermal_speed, 
+        parallel_heat_flux, parallel_heat_flux_updated, thermal_speed, temperature, 
         chodura_integral_lower, chodura_integral_upper, v_norm_fac,
         ddens_dz, ddens_dz_upwind, d2dens_dz2, dupar_dz, dupar_dz_upwind, d2upar_dz2,
-        dppar_dz, dppar_dz_upwind, d2ppar_dz2, dqpar_dz, dvth_dz, entropy_production,
+        dppar_dz, dppar_dz_upwind, d2ppar_dz2, dqpar_dz, dvth_dz, dT_dz, entropy_production,
         external_source_amplitude, external_source_density_amplitude,
         external_source_momentum_amplitude, external_source_pressure_amplitude,
         external_source_controller_integral, constraints_A_coefficient,
@@ -427,8 +433,8 @@ this function is only used once after initialisation
 the function used to update moments at run time is update_derived_moments! in time_advance.jl
 """
 function update_moments!(moments, ff_in, gyroavs::gyro_operators, vpa, vperp, z, r, composition,
-        r_spectral, geometry, scratch_dummy, z_advect)
-    if composition.gyrokinetic_ions
+        r_spectral, geometry, scratch_dummy, z_advect, collisions)
+    if composition.ion_physics == gyrokinetic_ions
         ff = scratch_dummy.buffer_vpavperpzrs_1 # the buffer array for the ion pdf -> make sure not to reuse this array below
         # fill buffer with ring-averaged F (gyroaverage at fixed position)
         gyroaverage_pdf!(ff,ff_in,gyroavs,vpa,vperp,z,r,composition)
@@ -466,12 +472,12 @@ function update_moments!(moments, ff_in, gyroavs::gyro_operators, vpa, vperp, z,
         end
         @views update_pperp_species!(moments.ion.pperp[:,:,is], ff[:,:,:,:,is], vpa, vperp, z, r)
         if moments.ion.qpar_updated[is] == false
-            @views update_qpar_species!(moments.ion.qpar[:,:,is],
+            @views update_ion_qpar_species!(moments.ion.qpar[:,:,is],
                                         moments.ion.dens[:,:,is],
                                         moments.ion.upar[:,:,is],
-                                        moments.ion.vth[:,:,is], ff[:,:,:,:,is], vpa,
+                                        moments.ion.vth[:,:,is],  moments.ion.dT_dz, ff[:,:,:,:,is], vpa,
                                         vperp, z, r, moments.evolve_density,
-                                        moments.evolve_upar, moments.evolve_ppar)
+                                        moments.evolve_upar, moments.evolve_ppar, composition.ion_physics, collisions)
             moments.ion.qpar_updated[is] = true
         end
     end
@@ -736,17 +742,18 @@ end
 """
 NB: the incoming pdf is the normalized pdf
 """
-function update_qpar!(qpar, qpar_updated, density, upar, vth, pdf, vpa, vperp, z, r,
-                      composition, evolve_density, evolve_upar, evolve_ppar)
+function update_ion_qpar!(qpar, qpar_updated, density, upar, vth, dT_dz, pdf, vpa, vperp, z, r,
+                          composition, ion_physics, collisions, evolve_density, evolve_upar, evolve_ppar)
     @boundscheck composition.n_ion_species == size(qpar,3) || throw(BoundsError(qpar))
 
     begin_s_r_z_region()
 
     @loop_s is begin
         if qpar_updated[is] == false
-            @views update_qpar_species!(qpar[:,:,is], density[:,:,is], upar[:,:,is],
-                                        vth[:,:,is], pdf[:,:,:,:,is], vpa, vperp, z, r,
-                                        evolve_density, evolve_upar, evolve_ppar)
+            @views update_ion_qpar_species!(qpar[:,:,is], density[:,:,is], upar[:,:,is],
+                                        vth[:,:,is], dT_dz, pdf[:,:,:,:,is], vpa, vperp, z, r,
+                                        evolve_density, evolve_upar, evolve_ppar, 
+                                        ion_physics, collisions)
             qpar_updated[is] = true
         end
     end
@@ -755,8 +762,25 @@ end
 """
 calculate the updated parallel heat flux (qpar) for a given species
 """
-function update_qpar_species!(qpar, density, upar, vth, ff, vpa, vperp, z, r, evolve_density,
-                              evolve_upar, evolve_ppar)
+function update_ion_qpar_species!(qpar, density, upar, vth, dT_dz, ff, vpa, vperp, z, r, evolve_density,
+                                  evolve_upar, evolve_ppar, ion_physics, collisions)
+    if ion_physics ∈ (drift_kinetic_ions, gyrokinetic_ions)
+        calculate_ion_qpar_from_pdf!(qpar, density, upar, vth, ff, vpa, vperp, z, r, evolve_density,
+                                     evolve_upar, evolve_ppar)
+    elseif ion_physics == coll_krook_ions
+        calculate_ion_qpar_from_coll_krook!(qpar, density, upar, vth, dT_dz, z, r, vperp, collisions, evolve_density, 
+                                            evolve_upar, evolve_ppar)
+    else
+        throw(ArgumentError("ion model $ion_physics not implemented for qpar calculation"))
+    end
+    return nothing
+end
+
+"""
+calculate parallel heat flux if ion composition flag is kinetic ions
+"""
+function calculate_ion_qpar_from_pdf!(qpar, density, upar, vth, ff, vpa, vperp, z, r, evolve_density, 
+                                      evolve_upar, evolve_ppar)
     @boundscheck r.n == size(ff, 4) || throw(BoundsError(ff))
     @boundscheck z.n == size(ff, 3) || throw(BoundsError(ff))
     @boundscheck vperp.n == size(ff, 2) || throw(BoundsError(ff))
@@ -793,7 +817,82 @@ function update_qpar_species!(qpar, density, upar, vth, ff, vpa, vperp, z, r, ev
     end
     return nothing
 end
+"""
+calculate parallel heat flux if ion composition flag is coll_krook fluid ions
+"""
+function calculate_ion_qpar_from_coll_krook!(qpar, density, upar, vth, dT_dz, z, r, vperp, collisions, evolve_density, evolve_upar, evolve_ppar)
+    # Note that this is a braginskii heat flux for ions using the krook operator. The full Fokker-Planck operator
+    # Braginskii heat flux is different! This also assumes one ion species, and so no friction between ions.
+    @boundscheck r.n == size(qpar, 2) || throw(BoundsError(qpar))
+    @boundscheck z.n == size(qpar, 1) || throw(BoundsError(qpar))
+
+    # calculate coll_krook heat flux. Currently only works for one ion species! (hence the 1 in dT_dz[iz,ir,1])
+    if evolve_density && evolve_upar && evolve_ppar
+        begin_r_z_region()
+        @loop_r_z ir iz begin
+            nu_ii = get_collision_frequency_ii(collisions, density[iz,ir], vth[iz,ir])
+            qpar[iz,ir] = -(1/2) * 5/4 * density[iz,ir] * vth[iz,ir]^2 /nu_ii * dT_dz[iz,ir,1]
+        end
+    else
+        throw(ArgumentError("coll_krook heat flux simulation requires evolve_density, 
+              evolve_upar and evolve_ppar to be true, since it is a purely fluid simulation"))
+    end
 
+    # add boundary condition to the heat flux, since now there is no distribution function 
+    # (in this case shape function) whose cutoff boundary condition can hold the parallel heat
+    # flux in check. See Stangeby textbook, equations (2.92) and (2.93), and the paragraph between.
+
+    if z.bc == "periodic"
+        # There's no wall boundary condition here, do nothing (qpar can be what it wants)
+        return nothing
+    end
+
+    begin_r_region()
+
+    if z.irank == 0 && (z.irank == z.nrank - 1)
+        z_indices = (1, z.n)
+    elseif z.irank == 0
+        z_indices = (1,)
+    elseif z.irank == z.nrank - 1
+        z_indices = (z.n,)
+    else
+        return nothing
+    end
+    # Stangeby (25.2) suggests that, when including kinetic effects, a value
+    # for gamma_i of around 2.5 is sensible.
+    # However, maybe for the purposes of this coll_krook scan, at very high
+    # collisionality we expect the distribution function of the ions at the
+    # sheath entrance to be close to a drifting maxwellian, in which case
+    # the original Stangeby (2.92) would be more appropriate. However, this
+    # also depends on whether we're 1V or 2V - as in 1V gamma_i = 2.5, 
+    # in 2V gamma_i = 3.5.
+
+    if vperp.n == 1
+        gamma_i = 2.5
+    else
+        gamma_i = 3.5
+    end
+    @loop_r ir begin
+        for iz ∈ z_indices
+            this_ppar = vth[iz,ir]^2 * density[iz,ir]/2.0
+            this_upar = upar[iz,ir]
+            this_dens = density[iz,ir]
+            particle_flux = this_dens * this_upar
+            T_i = vth[iz,ir]^2
+
+            # Stangeby (2.92)
+            total_heat_flux = gamma_i * T_i * particle_flux
+
+            # E.g. Helander&Sigmar (2.14), but in 1D we have no viscosity and only 3/2
+            # rather than 5/2.
+            conductive_heat_flux = total_heat_flux - 1.5 * this_ppar * this_upar - 
+                                                     0.5 * this_dens * this_upar^3
+
+            qpar[iz,ir] = conductive_heat_flux
+        end
+    end
+    return nothing
+end
 """
 runtime diagnostic routine for computing the Chodura ratio
 in a single species plasma with Z = 1
@@ -959,6 +1058,14 @@ function calculate_ion_moment_derivatives!(moments, scratch, scratch_dummy, z, z
                              buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z)
         @views derivative_z!(moments.ion.dvth_dz, vth, buffer_r_1,
                              buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z)
+
+        # calculate the z derivative of the ion temperature
+        @loop_s_r_z is ir iz begin
+            # store the temperature in dummy_zrs
+            dummy_zrs[iz,ir,is] = 2*ppar[iz,ir,is]/density[iz,ir,is]
+        end
+        @views derivative_z!(moments.ion.dT_dz, dummy_zrs, buffer_r_1,
+                            buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z)
     end
 end
 
@@ -1627,9 +1734,9 @@ end
 update velocity moments that are calculable from the evolved ion pdf
 """
 function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition,
-    r_spectral, geometry, gyroavs, scratch_dummy, z_advect, diagnostic_moments)
+    r_spectral, geometry, gyroavs, scratch_dummy, z_advect, collisions, diagnostic_moments)
 
-    if composition.gyrokinetic_ions
+    if composition.ion_physics == gyrokinetic_ions
         ff = scratch_dummy.buffer_vpavperpzrs_1
         # fill buffer with ring-averaged F (gyroaverage at fixed position)
         gyroaverage_pdf!(ff,new_scratch.pdf,gyroavs,vpa,vperp,z,r,composition)
@@ -1675,9 +1782,9 @@ function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composi
         rethrow(e)
     end
     # update the parallel heat flux
-    update_qpar!(moments.ion.qpar, moments.ion.qpar_updated, new_scratch.density,
-                 new_scratch.upar, moments.ion.vth, ff, vpa, vperp, z, r,
-                 composition, moments.evolve_density, moments.evolve_upar,
+    update_ion_qpar!(moments.ion.qpar, moments.ion.qpar_updated, new_scratch.density,
+                 new_scratch.upar, moments.ion.vth, moments.ion.dT_dz, ff, vpa, vperp, z, r,
+                 composition, composition.ion_physics, collisions, moments.evolve_density, moments.evolve_upar,
                  moments.evolve_ppar)
     # add further moments to be computed here
 
diff --git a/moment_kinetics/src/vpa_advection.jl b/moment_kinetics/src/vpa_advection.jl
index 8a04e4936..fb5ba59f7 100644
--- a/moment_kinetics/src/vpa_advection.jl
+++ b/moment_kinetics/src/vpa_advection.jl
@@ -67,7 +67,7 @@ end
                               fvec_in.density_neutral, fvec_in.uz_neutral,
                               fvec_in.pz_neutral)
     update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition,
-                            r_spectral, geometry, gyroavs, scratch_dummy, z_advect, false)
+                            r_spectral, geometry, gyroavs, scratch_dummy, z_advect, collisions, false)
     calculate_ion_moment_derivatives!(moments, new_scratch, scratch_dummy, z,
                                       z_spectral,
                                       num_diss_params.ion.moment_dissipation_coefficient)
diff --git a/moment_kinetics/test/gyroaverage_tests.jl b/moment_kinetics/test/gyroaverage_tests.jl
index b42564709..4f89f46d8 100644
--- a/moment_kinetics/test/gyroaverage_tests.jl
+++ b/moment_kinetics/test/gyroaverage_tests.jl
@@ -252,7 +252,7 @@ function gyroaverage_test(absolute_error; rhostar=0.1, pitch=0.5, ngrid=5, kr=2,
 end
 
 function create_test_composition()
-    input_dict = OptionsDict("composition" => OptionsDict("n_ion_species" => 1, "n_neutral_species" => 0, "gyrokinetic_ions" => true ) )
+    input_dict = OptionsDict("composition" => OptionsDict("n_ion_species" => 1, "n_neutral_species" => 0, "ion_physics" => "gyrokinetic_ions") )
     #println(input_dict)
     return get_species_input(input_dict)
 end
diff --git a/test_scripts/gyroaverage_test.jl b/test_scripts/gyroaverage_test.jl
index 5550e33b4..c1dcbb489 100644
--- a/test_scripts/gyroaverage_test.jl
+++ b/test_scripts/gyroaverage_test.jl
@@ -1,302 +1,302 @@
-export gyroaverage_test
-
-using Printf
-using Plots
-using LaTeXStrings
-using MPI
-using Measures
-using SpecialFunctions: besselj0
-
-import moment_kinetics
-using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.input_structs
-using moment_kinetics.geo: init_magnetic_geometry
-using moment_kinetics.communication
-using moment_kinetics.looping
-using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
-using moment_kinetics.gyroaverages: gyroaverage_pdf!
-using moment_kinetics.gyroaverages: gyroaverage_field!, init_gyro_operators
-using moment_kinetics.species_input: get_species_input
-using moment_kinetics.type_definitions: mk_float, mk_int
-
-function print_matrix(matrix,name::String,n::mk_int,m::mk_int)
-        println("\n ",name," \n")
-        for i in 1:n
-            for j in 1:m
-                @printf("%.2f ", matrix[i,j])
-            end
-            println("")
-        end
-        println("\n")
-    end
-    
-    function print_vector(vector,name::String,m::mk_int)
-        println("\n ",name," \n")
-        for j in 1:m
-            @printf("%.3f ", vector[j])
-        end
-        println("")
-        println("\n")
-    end
-
-
-function gyroaverage_test(;rhostar=0.1, pitch=0.5, ngrid=5, kr=2, kz=2, phaser=0.0, phasez=0.0, nelement=4, ngrid_vperp=3, nelement_vperp=1, Lvperp=3.0, ngrid_gyrophase=100, discretization="chebyshev_pseudospectral", r_bc="periodic", z_bc = "wall", plot_test_results=false)
-
-        #ngrid = 17
-        #nelement = 4
-        r_ngrid = ngrid #number of points per element 
-        r_nelement_local = nelement # number of elements per rank
-        r_nelement_global = r_nelement_local # total number of elements 
-        r_L = 1.0
-
-        z_ngrid = ngrid #number of points per element 
-        z_nelement_local = nelement # number of elements per rank
-        z_nelement_global = z_nelement_local # total number of elements 
-        z_L = 1.0
-
-        vperp_ngrid = ngrid_vperp #number of points per element 
-        vperp_nelement_local = nelement_vperp # number of elements per rank
-        vperp_nelement_global = vperp_nelement_local # total number of elements 
-        vperp_L = Lvperp
-        vperp_bc = "zero"
-        
-        vpa_ngrid = 1 #number of points per element 
-        vpa_nelement_local = 1 # number of elements per rank
-        vpa_nelement_global = vpa_nelement_local # total number of elements 
-        vpa_L = 1.0
-        vpa_bc = "" # should not be used
-        
-        gyrophase_ngrid = ngrid_gyrophase #number of points per element 
-        gyrophase_nelement_local = 1 # number of elements per rank
-        gyrophase_nelement_global = gyrophase_nelement_local # total number of elements 
-        gyrophase_discretization = "finite_difference"
-        gyrophase_L = 2.0*pi
-        
-        fd_option = "fourth_order_centered"
-        cheb_option = "matrix"
-        element_spacing_option = "uniform"
-        # create the 'input' struct containing input info needed to create a
-        # coordinate
-        
-        coords_input = OptionsDict(
-            "r"=>OptionsDict("ngrid"=>r_ngrid, "nelement"=>r_nelement_global,
-                             "nelement_local"=>r_nelement_local, "L"=>r_L,
-                             "discretization"=>discretization,
-                             "finite_difference_option"=>fd_option,
-                             "cheb_option"=>cheb_option, "bc"=>r_bc,
-                             "element_spacing_option"=>element_spacing_option),
-            "z"=>OptionsDict("ngrid"=>z_ngrid, "nelement"=>z_nelement_global,
-                             "nelement_local"=>z_nelement_local, "L"=>z_L,
-                             "discretization"=>discretization,
-                             "finite_difference_option"=>fd_option,
-                             "cheb_option"=>cheb_option, "bc"=>z_bc,
-                             "element_spacing_option"=>element_spacing_option),
-            "vperp"=>OptionsDict("ngrid"=>vperp_ngrid, "nelement"=>vperp_nelement_global,
-                                 "nelement_local"=>vperp_nelement_local, "L"=>vperp_L,
-                                 "discretization"=>discretization,
-                                 "finite_difference_option"=>fd_option,
-                                 "cheb_option"=>cheb_option, "bc"=>vperp_bc,
-                                 "element_spacing_option"=>element_spacing_option),
-            "vpa"=>OptionsDict("ngrid"=>vpa_ngrid, "nelement"=>vpa_nelement_global,
-                               "nelement_local"=>vpa_nelement_local, "L"=>vpa_L,
-                               "discretization"=>discretization,
-                               "finite_difference_option"=>fd_option,
-                               "cheb_option"=>cheb_option, "bc"=>vpa_bc,
-                               "element_spacing_option"=>element_spacing_option),
-            "gyrophase"=>OptionsDict("ngrid"=>gyrophase_ngrid,
-                                     "nelement"=>gyrophase_nelement_global,
-                                     "nelement_local"=>gyrophase_nelement_local,
-                                     "L"=>gyrophase_L, "discretization"=>discretization,
-                                     "finite_difference_option"=>fd_option,
-                                     "cheb_option"=>cheb_option, "bc"=>"periodic",
-                                     "element_spacing_option"=>element_spacing_option),
-        )
-        
-        # create the coordinate structs
-        r, r_spectral = define_coordinate(coords_input, "r"; collision_operator_dim=false)
-        z, z_spectral = define_coordinate(coords_input, "z"; collision_operator_dim=false)
-        vperp, vperp_spectral = define_coordinate(coords_input, "vperp";
-                                                  collision_operator_dim=false)
-        vpa, vpa_spectral = define_coordinate(coords_input, "vpa";
-                                              collision_operator_dim=false)
-        gyrophase, gyrophase_spectral = define_coordinate(coords_input, "gyrophase";
-                                                          collision_operator_dim=false)
-        
-        # create test geometry
-        #rhostar = 0.1 #rhostar of ions for ExB drift
-        option = "constant-helical"
-        #pitch = 1.0
-        DeltaB = 1.0
-        geometry_in = geometry_input(rhostar,option,pitch,DeltaB,0.0,0.0,0.0,0.0)
-        geometry = init_magnetic_geometry(geometry_in,z,r)
-        
-        # create test composition
-        composition = create_test_composition()
-            
-        # Set up MPI
-        initialize_comms!()
-        setup_distributed_memory_MPI(1,1,1,1)
-        looping.setup_loop_ranges!(block_rank[], block_size[];
-                                       s=composition.n_ion_species, sn=1,
-                                       r=r.n, z=z.n, vperp=vperp.n, vpa=vpa.n,
-                                       vzeta=1, vr=1, vz=1)
-                                       
-        # initialise the matrix for the gyroaverages
-        gyro = init_gyro_operators(vperp,z,r,gyrophase,geometry,composition)
-        # initialise a test field
-        phi = allocate_shared_float(z.n,r.n)
-        gphi = allocate_shared_float(vperp.n,z.n,r.n)
-        gphi_exact = allocate_float(vperp.n,z.n,r.n)
-        gphi_err = allocate_float(vperp.n,z.n,r.n)
-        begin_serial_region()
-        @serial_region begin 
-            fill_test_arrays!(phi,gphi_exact,vperp,z,r,geometry,kz,kr,phasez,phaser)
-        end
-        
-        # gyroaverage phi
-        gyroaverage_field!(gphi,phi,gyro,vperp,z,r,composition)
-        
-        # compute errors
-        begin_serial_region()
-        @serial_region begin
-            @. gphi_err = abs(gphi - gphi_exact)
-            println("Test gyroaverage_field!()")
-            for ivperp in 1:vperp.n
-                println("ivperp: ",ivperp," max(abs(gphi_err)): ",maximum(gphi_err[ivperp,:,:])," max(abs(gphi)): ",maximum(gphi[ivperp,:,:]))
-            end
-            println("")
-            if plot_test_results
-                @views heatmap(r.grid, z.grid, phi[:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
-                    windowsize = (360,240), margin = 15pt)
-                outfile = "phi_vs_r_z.pdf"
-                savefig(outfile)
-                println("Saved outfile: "*outfile)
-                for ivperp in 1:vperp.n
-                    @views heatmap(r.grid, z.grid, gphi[ivperp,:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
-                        windowsize = (360,240), margin = 15pt)
-                    outfile = "gphi_ivperp_"*string(ivperp)*"_vs_r_z.pdf"
-                    savefig(outfile)
-                    println("Saved outfile: "*outfile)
-                end
-            end
-        end
-        
-        # repeat the test for a pdf
-        # initialise a test field
-        nvpa = 1
-        n_ion_species = composition.n_ion_species
-        pdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
-        gpdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
-        gpdf_exact = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
-        gpdf_err = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
-        begin_serial_region()
-        @serial_region begin
-            fill_pdf_test_arrays!(pdf,gpdf_exact,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
-        end
-        
-        gyroaverage_pdf!(gpdf,pdf,gyro,vpa,vperp,z,r,composition)
-        # compute errors
-        begin_serial_region()
-        @serial_region begin
-            @. gpdf_err = abs(gpdf - gpdf_exact)
-            println("Test gyroaverage_pdf!()")
-            for is in 1:n_ion_species
-                for ivperp in 1:vperp.n
-                    for ivpa in 1:nvpa
-                        println("ivpa: ",ivpa," ivperp: ",ivperp," is: ",is," max(abs(gphi_err)): ",maximum(gpdf_err[ivpa,ivperp,:,:,is]),
-                         " max(abs(gpdf)): ",maximum(gpdf[ivpa,ivperp,:,:,is]))
-                    end
-                end
-            end
-            println("")
-        end
-        
-        finalize_comms!()
-end
-
-function create_test_composition()
-    electron_physics = boltzmann_electron_response
-    n_ion_species = 1
-    n_neutral_species = 0
-    n_species = n_ion_species + n_neutral_species
-    use_test_neutral_wall_pdf = false
-    # electron temperature over reference temperature
-    T_e = 1.0
-    # temperature at the entrance to the wall in terms of the electron temperature
-    T_wall = 1.0
-    # wall potential at z = 0
-    phi_wall = 0.0
-    # constant to test nonzero Er
-    Er_constant = 0.0
-    # ratio of the neutral particle mass to the ion particle mass
-    mn_over_mi = 1.0
-    # ratio of the electron particle mass to the ion particle mass
-    me_over_mi = 1.0/1836.0
-    # The ion flux reaching the wall that is recycled as neutrals is reduced by
-    # `recycling_fraction` to account for ions absorbed by the wall.
-    recycling_fraction = 1.0
-    gyrokinetic_ions = true
-    species_opts = OptionsDict("n_ion_species" => n_ion_species,
-                               "n_neutral_species" => n_neutral_species, "T_e" => T_e,
-                               "T_wall" => T_wall, "phi_wall" => phi_wall,
-                               "mn_over_mi" => mn_over_mi, "me_over_mi" => me_over_mi,
-                               "recycling_fraction" => recycling_fraction,
-                               "gyrokinetic_ions" => gyrokinetic_ions)
-    return get_species_input(OptionsDict("composition" => species_opts))
-end
-
-function fill_test_arrays!(phi,gphi,vperp,z,r,geometry,kz,kr,phasez,phaser)
-   for ir in 1:r.n
-      for iz in 1:z.n
-         Bmag = geometry.Bmag[iz,ir] 
-         bzeta = geometry.bzeta[iz,ir] 
-         rhostar = geometry.rhostar 
-         # convert integer "wavenumbers" to actual wavenumbers 
-         kkr = 2.0*pi*kr/r.L
-         kkz = 2.0*pi*kz/z.L
-         kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
-         
-         phi[iz,ir] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
-         for ivperp in 1:vperp.n
-            krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
-            # use that phi is a sum of Kronecker deltas in k space to write gphi
-            gphi[ivperp,iz,ir] = besselj0(krho)*phi[iz,ir]
-         end
-      end
-   end  
-   return nothing
-end
-
-function fill_pdf_test_arrays!(pdf,gpdf,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
-   for is in 1:composition.n_ion_species
-       for ir in 1:r.n
-          for iz in 1:z.n
-             Bmag = geometry.Bmag[iz,ir] 
-             bzeta = geometry.bzeta[iz,ir] 
-             rhostar = geometry.rhostar 
-             # convert integer "wavenumbers" to actual wavenumbers 
-             kkr = 2.0*pi*kr/r.L
-             kkz = 2.0*pi*kz/z.L
-             kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
-             
-             for ivperp in 1:vperp.n
-                for ivpa in 1:vpa.n
-                    pdf[ivpa,ivperp,iz,ir,is] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
-                    krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
-                    # use that pdf is a sum of Kronecker deltas in k space to write gpdf
-                    gpdf[ivpa,ivperp,iz,ir,is] = besselj0(krho)*pdf[ivpa,ivperp,iz,ir,is]
-                end
-             end
-          end
-       end
-   end  
-   return nothing
-end
-
-if abspath(PROGRAM_FILE) == @__FILE__
-    using Pkg
-    Pkg.activate(".")
-    # example function call with arguments for a successful test
-    #gyroaverage_test()
-    gyroaverage_test(rhostar=0.01,pitch=0.5,kr=2,kz=3,phaser=0.25*pi,phasez=0.5*pi,ngrid=9,nelement=8,ngrid_vperp=5,nelement_vperp=3,Lvperp=18.0,ngrid_gyrophase=100,r_bc="periodic",z_bc="periodic")
-end
+export gyroaverage_test
+
+using Printf
+using Plots
+using LaTeXStrings
+using MPI
+using Measures
+using SpecialFunctions: besselj0
+
+import moment_kinetics
+using moment_kinetics.coordinates: define_coordinate
+using moment_kinetics.input_structs
+using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.communication
+using moment_kinetics.looping
+using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
+using moment_kinetics.gyroaverages: gyroaverage_pdf!
+using moment_kinetics.gyroaverages: gyroaverage_field!, init_gyro_operators
+using moment_kinetics.species_input: get_species_input
+using moment_kinetics.type_definitions: mk_float, mk_int
+
+function print_matrix(matrix,name::String,n::mk_int,m::mk_int)
+        println("\n ",name," \n")
+        for i in 1:n
+            for j in 1:m
+                @printf("%.2f ", matrix[i,j])
+            end
+            println("")
+        end
+        println("\n")
+    end
+    
+    function print_vector(vector,name::String,m::mk_int)
+        println("\n ",name," \n")
+        for j in 1:m
+            @printf("%.3f ", vector[j])
+        end
+        println("")
+        println("\n")
+    end
+
+
+function gyroaverage_test(;rhostar=0.1, pitch=0.5, ngrid=5, kr=2, kz=2, phaser=0.0, phasez=0.0, nelement=4, ngrid_vperp=3, nelement_vperp=1, Lvperp=3.0, ngrid_gyrophase=100, discretization="chebyshev_pseudospectral", r_bc="periodic", z_bc = "wall", plot_test_results=false)
+
+        #ngrid = 17
+        #nelement = 4
+        r_ngrid = ngrid #number of points per element 
+        r_nelement_local = nelement # number of elements per rank
+        r_nelement_global = r_nelement_local # total number of elements 
+        r_L = 1.0
+
+        z_ngrid = ngrid #number of points per element 
+        z_nelement_local = nelement # number of elements per rank
+        z_nelement_global = z_nelement_local # total number of elements 
+        z_L = 1.0
+
+        vperp_ngrid = ngrid_vperp #number of points per element 
+        vperp_nelement_local = nelement_vperp # number of elements per rank
+        vperp_nelement_global = vperp_nelement_local # total number of elements 
+        vperp_L = Lvperp
+        vperp_bc = "zero"
+        
+        vpa_ngrid = 1 #number of points per element 
+        vpa_nelement_local = 1 # number of elements per rank
+        vpa_nelement_global = vpa_nelement_local # total number of elements 
+        vpa_L = 1.0
+        vpa_bc = "" # should not be used
+        
+        gyrophase_ngrid = ngrid_gyrophase #number of points per element 
+        gyrophase_nelement_local = 1 # number of elements per rank
+        gyrophase_nelement_global = gyrophase_nelement_local # total number of elements 
+        gyrophase_discretization = "finite_difference"
+        gyrophase_L = 2.0*pi
+        
+        fd_option = "fourth_order_centered"
+        cheb_option = "matrix"
+        element_spacing_option = "uniform"
+        # create the 'input' struct containing input info needed to create a
+        # coordinate
+        
+        coords_input = OptionsDict(
+            "r"=>OptionsDict("ngrid"=>r_ngrid, "nelement"=>r_nelement_global,
+                             "nelement_local"=>r_nelement_local, "L"=>r_L,
+                             "discretization"=>discretization,
+                             "finite_difference_option"=>fd_option,
+                             "cheb_option"=>cheb_option, "bc"=>r_bc,
+                             "element_spacing_option"=>element_spacing_option),
+            "z"=>OptionsDict("ngrid"=>z_ngrid, "nelement"=>z_nelement_global,
+                             "nelement_local"=>z_nelement_local, "L"=>z_L,
+                             "discretization"=>discretization,
+                             "finite_difference_option"=>fd_option,
+                             "cheb_option"=>cheb_option, "bc"=>z_bc,
+                             "element_spacing_option"=>element_spacing_option),
+            "vperp"=>OptionsDict("ngrid"=>vperp_ngrid, "nelement"=>vperp_nelement_global,
+                                 "nelement_local"=>vperp_nelement_local, "L"=>vperp_L,
+                                 "discretization"=>discretization,
+                                 "finite_difference_option"=>fd_option,
+                                 "cheb_option"=>cheb_option, "bc"=>vperp_bc,
+                                 "element_spacing_option"=>element_spacing_option),
+            "vpa"=>OptionsDict("ngrid"=>vpa_ngrid, "nelement"=>vpa_nelement_global,
+                               "nelement_local"=>vpa_nelement_local, "L"=>vpa_L,
+                               "discretization"=>discretization,
+                               "finite_difference_option"=>fd_option,
+                               "cheb_option"=>cheb_option, "bc"=>vpa_bc,
+                               "element_spacing_option"=>element_spacing_option),
+            "gyrophase"=>OptionsDict("ngrid"=>gyrophase_ngrid,
+                                     "nelement"=>gyrophase_nelement_global,
+                                     "nelement_local"=>gyrophase_nelement_local,
+                                     "L"=>gyrophase_L, "discretization"=>discretization,
+                                     "finite_difference_option"=>fd_option,
+                                     "cheb_option"=>cheb_option, "bc"=>"periodic",
+                                     "element_spacing_option"=>element_spacing_option),
+        )
+        
+        # create the coordinate structs
+        r, r_spectral = define_coordinate(coords_input, "r"; collision_operator_dim=false)
+        z, z_spectral = define_coordinate(coords_input, "z"; collision_operator_dim=false)
+        vperp, vperp_spectral = define_coordinate(coords_input, "vperp";
+                                                  collision_operator_dim=false)
+        vpa, vpa_spectral = define_coordinate(coords_input, "vpa";
+                                              collision_operator_dim=false)
+        gyrophase, gyrophase_spectral = define_coordinate(coords_input, "gyrophase";
+                                                          collision_operator_dim=false)
+        
+        # create test geometry
+        #rhostar = 0.1 #rhostar of ions for ExB drift
+        option = "constant-helical"
+        #pitch = 1.0
+        DeltaB = 1.0
+        geometry_in = geometry_input(rhostar,option,pitch,DeltaB,0.0,0.0,0.0,0.0)
+        geometry = init_magnetic_geometry(geometry_in,z,r)
+        
+        # create test composition
+        composition = create_test_composition()
+            
+        # Set up MPI
+        initialize_comms!()
+        setup_distributed_memory_MPI(1,1,1,1)
+        looping.setup_loop_ranges!(block_rank[], block_size[];
+                                       s=composition.n_ion_species, sn=1,
+                                       r=r.n, z=z.n, vperp=vperp.n, vpa=vpa.n,
+                                       vzeta=1, vr=1, vz=1)
+                                       
+        # initialise the matrix for the gyroaverages
+        gyro = init_gyro_operators(vperp,z,r,gyrophase,geometry,composition)
+        # initialise a test field
+        phi = allocate_shared_float(z.n,r.n)
+        gphi = allocate_shared_float(vperp.n,z.n,r.n)
+        gphi_exact = allocate_float(vperp.n,z.n,r.n)
+        gphi_err = allocate_float(vperp.n,z.n,r.n)
+        begin_serial_region()
+        @serial_region begin 
+            fill_test_arrays!(phi,gphi_exact,vperp,z,r,geometry,kz,kr,phasez,phaser)
+        end
+        
+        # gyroaverage phi
+        gyroaverage_field!(gphi,phi,gyro,vperp,z,r,composition)
+        
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gphi_err = abs(gphi - gphi_exact)
+            println("Test gyroaverage_field!()")
+            for ivperp in 1:vperp.n
+                println("ivperp: ",ivperp," max(abs(gphi_err)): ",maximum(gphi_err[ivperp,:,:])," max(abs(gphi)): ",maximum(gphi[ivperp,:,:]))
+            end
+            println("")
+            if plot_test_results
+                @views heatmap(r.grid, z.grid, phi[:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
+                    windowsize = (360,240), margin = 15pt)
+                outfile = "phi_vs_r_z.pdf"
+                savefig(outfile)
+                println("Saved outfile: "*outfile)
+                for ivperp in 1:vperp.n
+                    @views heatmap(r.grid, z.grid, gphi[ivperp,:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
+                        windowsize = (360,240), margin = 15pt)
+                    outfile = "gphi_ivperp_"*string(ivperp)*"_vs_r_z.pdf"
+                    savefig(outfile)
+                    println("Saved outfile: "*outfile)
+                end
+            end
+        end
+        
+        # repeat the test for a pdf
+        # initialise a test field
+        nvpa = 1
+        n_ion_species = composition.n_ion_species
+        pdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_exact = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_err = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        begin_serial_region()
+        @serial_region begin
+            fill_pdf_test_arrays!(pdf,gpdf_exact,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+        end
+        
+        gyroaverage_pdf!(gpdf,pdf,gyro,vpa,vperp,z,r,composition)
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gpdf_err = abs(gpdf - gpdf_exact)
+            println("Test gyroaverage_pdf!()")
+            for is in 1:n_ion_species
+                for ivperp in 1:vperp.n
+                    for ivpa in 1:nvpa
+                        println("ivpa: ",ivpa," ivperp: ",ivperp," is: ",is," max(abs(gphi_err)): ",maximum(gpdf_err[ivpa,ivperp,:,:,is]),
+                         " max(abs(gpdf)): ",maximum(gpdf[ivpa,ivperp,:,:,is]))
+                    end
+                end
+            end
+            println("")
+        end
+        
+        finalize_comms!()
+end
+
+function create_test_composition()
+    electron_physics = boltzmann_electron_response
+    n_ion_species = 1
+    n_neutral_species = 0
+    n_species = n_ion_species + n_neutral_species
+    use_test_neutral_wall_pdf = false
+    # electron temperature over reference temperature
+    T_e = 1.0
+    # temperature at the entrance to the wall in terms of the electron temperature
+    T_wall = 1.0
+    # wall potential at z = 0
+    phi_wall = 0.0
+    # constant to test nonzero Er
+    Er_constant = 0.0
+    # ratio of the neutral particle mass to the ion particle mass
+    mn_over_mi = 1.0
+    # ratio of the electron particle mass to the ion particle mass
+    me_over_mi = 1.0/1836.0
+    # The ion flux reaching the wall that is recycled as neutrals is reduced by
+    # `recycling_fraction` to account for ions absorbed by the wall.
+    recycling_fraction = 1.0
+    ion_physics = gyrokinetic_ions
+    species_opts = OptionsDict("n_ion_species" => n_ion_species,
+                               "n_neutral_species" => n_neutral_species, "T_e" => T_e,
+                               "T_wall" => T_wall, "phi_wall" => phi_wall,
+                               "mn_over_mi" => mn_over_mi, "me_over_mi" => me_over_mi,
+                               "recycling_fraction" => recycling_fraction,
+                               "ion_physics" => ion_physics)
+    return get_species_input(OptionsDict("composition" => species_opts))
+end
+
+function fill_test_arrays!(phi,gphi,vperp,z,r,geometry,kz,kr,phasez,phaser)
+   for ir in 1:r.n
+      for iz in 1:z.n
+         Bmag = geometry.Bmag[iz,ir] 
+         bzeta = geometry.bzeta[iz,ir] 
+         rhostar = geometry.rhostar 
+         # convert integer "wavenumbers" to actual wavenumbers 
+         kkr = 2.0*pi*kr/r.L
+         kkz = 2.0*pi*kz/z.L
+         kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+         
+         phi[iz,ir] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+         for ivperp in 1:vperp.n
+            krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+            # use that phi is a sum of Kronecker deltas in k space to write gphi
+            gphi[ivperp,iz,ir] = besselj0(krho)*phi[iz,ir]
+         end
+      end
+   end  
+   return nothing
+end
+
+function fill_pdf_test_arrays!(pdf,gpdf,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+   for is in 1:composition.n_ion_species
+       for ir in 1:r.n
+          for iz in 1:z.n
+             Bmag = geometry.Bmag[iz,ir] 
+             bzeta = geometry.bzeta[iz,ir] 
+             rhostar = geometry.rhostar 
+             # convert integer "wavenumbers" to actual wavenumbers 
+             kkr = 2.0*pi*kr/r.L
+             kkz = 2.0*pi*kz/z.L
+             kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+             
+             for ivperp in 1:vperp.n
+                for ivpa in 1:vpa.n
+                    pdf[ivpa,ivperp,iz,ir,is] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+                    krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+                    # use that pdf is a sum of Kronecker deltas in k space to write gpdf
+                    gpdf[ivpa,ivperp,iz,ir,is] = besselj0(krho)*pdf[ivpa,ivperp,iz,ir,is]
+                end
+             end
+          end
+       end
+   end  
+   return nothing
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    using Pkg
+    Pkg.activate(".")
+    # example function call with arguments for a successful test
+    #gyroaverage_test()
+    gyroaverage_test(rhostar=0.01,pitch=0.5,kr=2,kz=3,phaser=0.25*pi,phasez=0.5*pi,ngrid=9,nelement=8,ngrid_vperp=5,nelement_vperp=3,Lvperp=18.0,ngrid_gyrophase=100,r_bc="periodic",z_bc="periodic")
+end