We present the QCD predictions for the azimuthal $cos 2varphi$ asymmetry in charm leptoproduction for the kinematics of the COMPASS experiment at CERN. The asymmetry is predicted to be large, about 15%. The radiative corrections to the QCD predictions for the $cos 2varphi$ distribution are estimated to be small, less than 10%. Our calculations show that the azimuthal asymmetry in charm production is well defined in pQCD: it is stable both perturbatively and parametrically, and practically insensitive to theoretical uncertainties in the input parameters. We analyze the nonperturbative contributions to the $cos 2varphi$ distribution due to the gluon transverse motion in the target and the $c$-quark fragmentation. Because of the $c$-quark low mass, the nonperturbative contributions are expected to be sizable, about (30--40)%. We conclude that extraction of the azimuthal asymmetries from available COMPASS data will provide valuable information about the transverse momentum dependent distribution of the gluon in the proton and the $c$-quark hadronization mechanism. Finally, we discuss the $cos 2varphi$ asymmetry as a probe of the gluonic analogue of the Boer-Mulders function, $h_{1}^{perp g}$, describing the linear polarization of gluons inside unpolarized proton.