diff --git a/src/time_advance.jl b/src/time_advance.jl
index 55499df5af..35f453a13a 100644
--- a/src/time_advance.jl
+++ b/src/time_advance.jl
@@ -911,29 +911,39 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
                     layout = fig[irow,1] = GridLayout()
                     ax = Axis(layout[1,1], xlabel="z", ylabel="density", width=ax_width, height=ax_height)
                     plot_1d(z.grid, moments.charged.dens[:,1,1]; ax=ax, label="ion")
-                    plot_1d(z.grid, moments.neutral.dens[:,1,1]; ax=ax, label="neutral")
+                    if composition.n_neutral_species > 0
+                        plot_1d(z.grid, moments.neutral.dens[:,1,1]; ax=ax, label="neutral")
+                    end
                     #axislegend(ax)
                     ax = Axis(layout[1,2], xlabel="z", ylabel="upar", width=ax_width, height=ax_height)
                     plot_1d(z.grid, moments.charged.upar[:,1,1]; ax=ax, label="ion")
-                    plot_1d(z.grid, moments.neutral.uz[:,1,1]; ax=ax, label="neutral")
+                    if composition.n_neutral_species > 0
+                        plot_1d(z.grid, moments.neutral.uz[:,1,1]; ax=ax, label="neutral")
+                    end
                     #axislegend(ax)
 
                     irow += 1
                     layout = fig[irow,1] = GridLayout()
                     ax = Axis(layout[1,1], xlabel="z", ylabel="ppar", width=ax_width, height=ax_height)
                     plot_1d(z.grid, moments.charged.ppar[:,1,1]; ax=ax, label="ion")
-                    plot_1d(z.grid, moments.neutral.pz[:,1,1]; ax=ax, label="neutral")
+                    if composition.n_neutral_species > 0
+                        plot_1d(z.grid, moments.neutral.pz[:,1,1]; ax=ax, label="neutral")
+                    end
                     #axislegend(ax)
                     ax = Axis(layout[1,2], xlabel="z", ylabel="vth", width=ax_width, height=ax_height)
                     plot_1d(z.grid, moments.charged.vth[:,1,1]; ax=ax, label="ion")
-                    plot_1d(z.grid, moments.neutral.vth[:,1,1]; ax=ax, label="neutral")
+                    if composition.n_neutral_species > 0
+                        plot_1d(z.grid, moments.neutral.vth[:,1,1]; ax=ax, label="neutral")
+                    end
                     #axislegend(ax)
 
                     irow += 1
                     layout = fig[irow,1] = GridLayout()
                     ax = Axis(layout[1,1], xlabel="z", ylabel="qpar", width=ax_width, height=ax_height)
                     plot_1d(z.grid, moments.charged.qpar[:,1,1]; ax=ax, label="ion")
-                    plot_1d(z.grid, moments.neutral.qz[:,1,1]; ax=ax, label="neutral")
+                    if composition.n_neutral_species > 0
+                        plot_1d(z.grid, moments.neutral.qz[:,1,1]; ax=ax, label="neutral")
+                    end
                     #axislegend(ax)
 
                     resize_to_layout!(fig)