In order to clarify the role of neutral dynamics in the Jovian magnetosphere-ionosphere-thermosphere coupling system, we have developed a new numerical model that includes the effect of neutral dynamics on the coupling current. The model calculates axisymmetric thermospheric dynamics and ion composition by considering fundamental physical and chemical processes. The ionospheric Pedersen current is obtained from the thermospheric and ionospheric parameters. The model simultaneously solves the torque equations of the magnetospheric plasma due to radial currents flowing at the magnetospheric equator, which enables us to update the electric field projected onto the ionosphere and the field-aligned currents (FACs) depending upon the thermospheric dynamics. The self-consistently calculated temperature and ion velocity are consistent with observations. The estimated neutral wind field captures the zonally averaged characteristics in previous three-dimensional models. The energy extracted from the planetary rotation is mainly used for magnetospheric plasma acceleration below 73.5° latitude while consumed in the upper atmosphere, mainly by Joule heating at above 73.5° latitude. The neutral wind dynamics contributes to a reduction in the electric field of 22% compared with the case of neutral rigid corotation. About 90% of this reduction is attributable to neutral winds below the 550-km altitude in the auroral region. The calculated radial current in the equatorial magnetosphere is smaller than observations. This indicates that the enhancement of the background conductance and/or the additional radial current at the outer boundary would be expected to reproduce the observed current.