Bug in oderhs.m gas thermal conduction calculation for running with diffusive neutrals + impurities

[mwigram]

Describe the bug

In oderhs.m, lines 3102-3147 code that calculate the gas thermal conduction coefficients (populating hcxg/hcyg, appearing to be an approximation for molecules) uses terms hard-coded to fluid species ifld=2 (e.g. ni(ix,iy,2) for variables "naavex/y", intended to be for neutral atoms (e.g. lines 3118 and 3123).

This is correct for when inertial neutrals model is active, but is incorrect for diffusive neutrals (when atoms do not populate the ni variable). This has been accounted for in the code by the condition "if (nisp >= 2) ", line 3102.

However, this fix only works when there are no impurity species (or additional plasma species in general). For example, running deuterium plasma with diffusive neutrals + Ne impurity, ni(ix,iy,2) does not correspond to neutral atoms but instead to the first neon charge state. But the condition nisp >= 2 is satisfied, so the gas thermal conduction coefficients code runs and populates hcxg/hcyg, but using the neon density. The resulting error is mostly small (not approaching MW scale), but hard to say if it will always be.

To Reproduce
Steps to reproduce the behavior:

1. Run a case with multiple ion fluid species and diffusive neutrals
2. Check hcxg/hcyg to show they are none-zero

Suggested fix
Change line 3102 from "if (nisp >= 2) " to "if (isupgon(1) .eq. 1)" - i.e. only runs when inertial neutral model is active, rather than when there are more than 2 ion species present.

[tdrognlien]

I think this could be fixed by replacing ni(,,2) by ng(,,1) for naavex and replacing niy0,1(..,2) by ngy0,1(,,1) for naavey and then removing the if(nisp = 2) test. There is a duality in the variables ng(..1) = ni(,,2) if isupgon=1 even though then ng(,,1) is not a primary evolved variable. Note sure why this wasn't done before (i.e., the suggested fix), so should think a little more before making the change.