diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl
index de20e7a44..400270234 100644
--- a/moment_kinetics/src/external_sources.jl
+++ b/moment_kinetics/src/external_sources.jl
@@ -1293,7 +1293,7 @@ source amplitude.
         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]) *
+                    (0.5 * ion_source_settings.source_T + upar[iz,ir]^2 - ppar[iz,ir]/density[iz,ir]) *
                     ion_source_settings.source_strength *
                     ion_source_settings.r_amplitude[ir] *
                     ion_source_settings.z_amplitude[iz]
@@ -1397,7 +1397,7 @@ source amplitude.
         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) *
+                    ((0.5 * ion_source_settings.source_T + 2 * upar[iz,ir]^2 - ppar[iz,ir]/density[iz,ir]) *
                     amplitude) * ion_source_settings.controller_source_profile[iz,ir]
             end
         end
diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl
index 4db3fd743..4330a682d 100644
--- a/moment_kinetics/src/load_data.jl
+++ b/moment_kinetics/src/load_data.jl
@@ -4295,7 +4295,7 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
         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
+        variable = @. -(1/2) * 3/2 * 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...)
diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl
index 5005ef240..413e33c01 100644
--- a/moment_kinetics/src/time_advance.jl
+++ b/moment_kinetics/src/time_advance.jl
@@ -1979,6 +1979,9 @@ function  time_advance!(pdf, scratch, scratch_implicit, scratch_electron, t_para
                                 "t = ", rpad(string(round(t_params.t[], sigdigits=6)), 7), "  ",
                                 "nstep = ", rpad(string(t_params.step_counter[]), 7), "  ",
                                 Dates.format(now(), dateformat"H:MM:SS"))
+                    println("ion temperature at middle: ", moments.ion.vth[65,1,1]^2)
+                    println("ion density at middle: ", moments.ion.dens[65,1,1])
+                    println("ion source amplitude: ", moments.ion.external_source_amplitude[65,1,:])
                         flush(stdout)
                     end
                 end
diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl
index c58be5cbf..770ecf6fa 100644
--- a/moment_kinetics/src/velocity_moments.jl
+++ b/moment_kinetics/src/velocity_moments.jl
@@ -831,7 +831,7 @@ function calculate_ion_qpar_from_coll_krook!(qpar, density, upar, vth, dT_dz, z,
         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]
+            qpar[iz,ir] = -(1/2) * 3/2 * 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,