A rigorous and modular Monte Carlo radiative transfer model MCRT has been developed to compute radiance for realistic cloudy atmosphere over sea surfaces in the visible and near-IR spectra. The simulation methods of the sample's trajectories are consistently described in the whole simulation spatial domain, which consists of a number of voxel cells for atmosphere and square cells for the sea surface. An inverse transformation method is applied in an atmospheric scattering event simulation and a rigorous accept/reject method is formulated for simulating a reflection event at Cox-Munk (CM) anisotropically and Nakajima-Tanaka (NT) isotropically roughened sea surface. Both methods have been implemented within the model and were tested as two individual modules achieving high accuracy. The model as a whole system was also validated by other codes. The mean difference is about 0.084% and 0.528% in comparison to the Spherical Harmonics code (SHARM) and Monte Carlo Atmospheric Radiative Transfer Simulator (MCARaTS) respectively for 1D atmosphere and the NT model. The comparisons of the Spherical Harmonic Discrete Ordinate Method (SHDOM) show that agreements are obtained in the sun glint regions for 1D atmosphere and the CM model and the mean differences are below 0.478% for a 3D cloud field with lambertian. In general, the major advantage of MCRT is that it could simulate more accurate reflection direction at sea surfaces in a shorter time and hence more accurate radiance in complex cloudy atmosphere over sea surfaces with high numerical efficiency especially on a small PC.