Accurately simulating high-speed gas flows requires meshes of high quality. Post shock conditions depend strongly on mesh resolution and must be adequately resolved for accurate force and energy predictions at domain boundaries. Furthermore, the inclusion of multiple chemical constituents and thermochemical nonequilibrium increases the size of the linear system and introduces stiff source terms, placing a premium on efficient solution strategies. This paper proposes such a strategy via adjoint-based goal-oriented mesh adaptation. The adaptation criterion is tuned to a particular functional of interest, enabling optimal grid refinement for continuum, multi-species flow in thermochemical nonequilibrium. The formulation of the direct and continuous adjoint problems is presented, including derivation of the adjoint boundary conditions for pressure-based functionals in multi-species plasmas. Details of the numerical implementation in a general, unstructured CFD solver are also included. Results are presented for an entry velocity, blunt-body geometry in inviscid flow with a two species Nitrogen chemical model. Three adaptation schemes are compared: full grid, gradient-based, and adjoint-based. The goal-oriented adjoint approach exhibits the fastest grid convergence of all methods and higher solution quality when compared to gradient-based adaptation.