  In a hydrodynamic model, we have studied $J/\psi$ production in Au+Au/Cu+Cu
collisions at RHIC energy $\sqrt{s}$=200 GeV. At the initial time, $J/\psi$'s
are randomly distributed in the fluid. As the fluid evolve in time, the free
streaming $J/\psi$'s are dissolved if the local fluid temperature exceeds a
threshold temperature $T_{J/\psi}$. Sequential melting of charmonium states
($\chi_c$, $\psi\prime$ and $J/\psi$), with melting temperatures
$T_{\chi_c}=T_{\psi\prime} \approx 1.2T_c$, $T_{J/\psi} \approx2T_c$ and
feed-down fraction $F\approx 0.3$, explains the PHENIX data on the centrality
dependence of $J/\psi$ suppression in Au+Au collisions. $J/\psi$ $p_T$ spectra
and the nuclear modification factor in Au+Au collisions are also well explained
in the model. The model however over predict centrality dependence of $J/\psi$
suppression in Cu+Cu collisions by 20-30%. The $J/\psi$ $p_T$ spectra are under
predicted by 20-30%. The model predict that in central Pb+Pb collisions at LHC
energy, $\sqrt{s}$=5500 GeV, $J/\psi$'s are suppressed by a factor of $\sim$
10. The model predicted $J/\psi$ $p_T$ distribution in Pb+Pb collisions at LHC
is similar to that in Au+Au collisions at RHIC.
