  A description is given for preserving ${\bmsy\nabla}\cdot{\vec B}=0$ in a
magnetohydrodynamic (MHD) code that employs the upwind, Total Variation
Diminishing (TVD) scheme and the Strang-type operator splitting for
multi-dimensionality. The method is based on the staggered mesh technique to
constrain the transport of magnetic field: the magnetic field components are
defined at grid interfaces with their advective fluxes on grid edges, while
other quantities are defined at grid centers. The magnetic field at grid
centers for the upwind step is calculated by interpolating the values from grid
interfaces. The advective fluxes on grid edges for the magnetic field evolution
are calculated from the upwind fluxes at grid interfaces. Then, the magnetic
field can be maintained with ${\bmsy\nabla}\cdot{\vec B}=0$ exactly, if this is
so initially, while the upwind scheme is used for the update of fluid
quantities. The correctness of the code is demonstrated through tests comparing
numerical solutions either with analytic solutions or with numerical solutions
from the code using an explicit divergence-cleaning method. Also the robustness
is shown through tests involving realistic astrophysical problems.
