This paper presents a new method for solving the hydrodynamic (HD) model in submicron semiconductor device simulation. The main feature of this method is that the Poisson, current-continuity and energy-balance equations in the HD model are all expressed in self-adjoint forms through a set of new Slotboom-like variables. As a consequence, the discretization results in a system of finite-difference equations with a diagonally dominant coefficient matrix for each HD equation. The simultaneous HD equations are decoupled by using the Gummel block iteration method. To solve each equation, a fixed-point iteration technique is employed which explicitly updates the state variables at each spatial mesh-point. In addition to avoiding direct solution of large matrix equations, the diagonal dominance guarantees that each HD equation will converge for any initial value. We demonstrate the method by simulating a 2-D submicron MOSFET, and by comparison with Monte Carlo calculations. Excellent numerical convergence, stability, and efficiency are observed.