  We present a detailed study of the charmonium spectrum using anisotropic
lattice QCD. We first derive a tree-level improved clover quark action on the
anisotropic lattice for arbitrary quark mass. The heavy quark mass dependences
of the improvement coefficients, i.e. the ratio of the hopping parameters
$\zeta=K_t/K_s$ and the clover coefficients $c_{s,t}$, are examined at the tree
level. We then compute the charmonium spectrum in the quenched approximation
employing $\xi = a_s/a_t = 3$ anisotropic lattices. Simulations are made with
the standard anisotropic gauge action and the anisotropic clover quark action
at four lattice spacings in the range $a_s$=0.07-0.2 fm. The clover
coefficients $c_{s,t}$ are estimated from tree-level tadpole improvement. On
the other hand, for the ratio of the hopping parameters $\zeta$, we adopt both
the tree-level tadpole-improved value and a non-perturbative one. We calculate
the spectrum of S- and P-states and their excitations. The results largely
depend on the scale input even in the continuum limit, showing a quenching
effect. When the lattice spacing is determined from the $1P-1S$ splitting, the
deviation from the experimental value is estimated to be $\sim$30% for the
S-state hyperfine splitting and $\sim$20% for the P-state fine structure. Our
results are consistent with previous results at $\xi = 2$ obtained by Chen when
the lattice spacing is determined from the Sommer scale $r_0$. We also address
the problem with the hyperfine splitting that different choices of the clover
coefficients lead to disagreeing results in the continuum limit.
