  Using 1000 ray-tracing simulations for a {\Lambda}-dominated cold dark model
in Sato et al. (2009), we study the covariance matrix of cosmic shear
correlation functions, which is the standard statistics used in the previous
measurements. The shear correlation function of a particular separation angle
is affected by Fourier modes over a wide range of multipoles, even beyond a
survey area, which complicates the analysis of the covariance matrix. To
overcome such obstacles we first construct Gaussian shear simulations from the
1000 realizations, and then use the Gaussian simulations to disentangle the
Gaussian covariance contribution to the covariance matrix we measured from the
original simulations. We found that an analytical formula of Gaussian
covariance overestimates the covariance amplitudes due to an effect of finite
survey area. Furthermore, the clean separation of the Gaussian covariance
allows to examine the non-Gaussian covariance contributions as a function of
separation angles and source redshifts. For upcoming surveys with typical
source redshifts of z_s=0.6 and 1.0, the non-Gaussian contribution to the
diagonal covariance components at 1 arcminute scales is greater than the
Gaussian contribution by a factor of 20 and 10, respectively. Predictions based
on the halo model qualitatively well reproduce the simulation results, however
show a sizable disagreement in the covariance amplitudes. By combining these
simulation results we develop a fitting formula to the covariance matrix for a
survey with arbitrary area coverage, taking into account effects of the
finiteness of survey area on the Gaussian covariance.
