A three dimensional numerical model in the σ coordinate system is developed to study the problem of waves. Turbulence effects are modeled by a subgrid scale (SGS) model with the concept of large eddy simulation (LES). The σ coordinate transformation is introduced to map the irregular physical domain of the wavy free surface and uneven bottom onto the regular computational domain of the shape of rectangular prism. The operator splitting method, which splits the solution procedure into the advection, diffusion, and propagation steps, is used to solve the modified Navier Stokes Equation. The model is used to simulate the propagation of solitary wave and wave passing over a submerged breakwater. Numerical results are compared with available analytical solutions and experimental data in terms of velocity profiles, free surface displacement, and energy conservation. Good agreement is obtained. The method is proved to be of high accuracy and efficiency in simulating surface wave propagation and wave structure interaction. It is suitable for the large and irregular physical domain, and requiring the non uniform grid system. The present work provides a foundation for further studies of random waves, wave structure interaction, wave discharge interaction, etc. A three dimensional numerical model in the σ coordinate system is developed to study the problem of waves. Turbulence effects are modeled by a subgrid scale (SGS) model with the concept of large eddy simulation (LES). The σ coordinate transformation is introduced to map the irregular physical domain of the wavy free surface and uneven bottom onto the regular computational domain of the shape of rectangular prism. the operator splitting method, which splits the solution procedure into the advection, diffusion, and propagation steps, is used to solve the modified Navier Stokes Equation. The model is used to simulate the propagation of solitary wave and wave passing over a submerged breakwater. Numerical results are compared with available analytical solutions and experimental data in terms of velocity profiles, free surface displacement, and energy conservation. Good agreement is obtained. The method is proved to be of high accuracy and efficiency in simulating surface wav It is suitable for the large and irregular physical domain, and requires the non uniform grid system. The present work provides a foundation for further studies of random waves, wave structure interaction, wave discharge interaction, etc.