The way we output components of contravariant/covariant vectors is confusing, and potentially incorrect in terms of units. Referring to Jbgradpsi in gpec_profile_output.nc, the contravariant representation of a Jacobian-weighted radial field is:
$$
\displaystyle J\vec{b} = (Jb)^i\vec{e}_{i} = (J\cdot \vec{b})\cdot \nabla \psi \frac{ \partial \vec{R} }{ \partial \psi }
$$
The output lists Jbgradpsi in terms of Tesla. We know that the vector $\vec{b}$ should have units of Tesla. To keep it simple, let's consider the default (hamada) coordinates. The hamada jacobian's units are
$$
J = \left( \vec{B}\cdot\nabla\theta \right)^{-1} \sim \left( [T][m^{-1}][] \right)^{-1} \sim [T^{-1}][m]
$$
so the quantity $J \vec{b}$ should have length units. The quantity $\nabla \psi$ may have units of $[T][m]$ or $[m^{-1}]$ depending on whether here $\psi$ refers to poloidal or normalized poloidal flux. Then comes the question of whether Jbgradpsi includes the units of its vector, $\frac{ \partial \vec{R} }{ \partial \psi }$, which also depends on which $\psi$ we refer to. Someone should go through and carefully determine what units these quantities have, and add comments to the code where they are initialized. This work needs to be done anyway for the Julia conversion. This would help us be confident in those quantitative predictions made by the code that depend on units. It is also relevant to determining the units of $\delta$ (see #189).
