Monte Carlo (MC) is a significant technique for finding the radiative transfer equation (RTE) solution. Nowadays, the Henyey-Greenstein (HG) scattering phase function (spf) has been widely used in most studies during the core procedure of randomly choosing scattering angles in oceanographic lidar MC simulations. However, the HG phase function does not work well at small or large scattering angles. Other spfs work well, e.g., Fournier-Forand phase function (FF); however, solving the cumulative distribution function (cdf) of the scattering phase function (even if possible) would result in a complicated formula. To avoid the above-mentioned problems, we present a semi-analytic MC radiative transfer model in this paper, which uses the cdf equation to build up a lookup table (LUT) of ψ vs. P Ψ ( ψ ) to determine scattering angles for various spfs (e.g., FF, Petzold measured particle phase function, and so on). Moreover, a lidar geometric model for analytically estimating the probability of photon scatter back to a remote receiver was developed; in particular, inhomogeneous layers are divided into voxels with different optical properties; therefore, it is useful for inhomogeneous water. First, the simulations between the inverse function method for HG cdf and the LUT method for FF cdf were compared. Then, multiple scattering and wind-driven sea surface condition effects were studied. Finally, we compared our simulation results with measurements of airborne lidar. The mean relative errors between simulation and measurements in inhomogeneous water are within 14% for the LUT method and within 22% for the inverse cdf (ICDF) method. The results suggest feasibility and effectiveness of our simulation model.