A large‐scale Monte Carlo program has been written to calculate the complete four‐component Stokes vector at any region in a fully inhomogeneous atmosphere‐ocean system with inclusion of a stochastic interface. The program uses as input the Mueller matrices for the aerosols in the atmosphere as well as the hydrosols in the ocean. The Mueller matrix for the dielectric interface is also accurately accounted for. Rigorous testing was performed by comparing the Monte Carlo technique with a completely independent invariant imbedding approach; the agreement was excellent. To study the effect of the dielectric interface between atmosphere and ocean we set up a model system consisting of only Rayleigh scattering with a calm ocean, although the program can handle a stochastic boundary as well. Results are presented for both the radiance and degree of polarization at various locations in the atmosphere‐ocean system as the refractive index of the dielectric interface varies from 1.338, typical of ocean water, to a value of 1.0 which makes the interface disappear.