A coastal ocean hydrodynamics modelling system has been created, with potential for future improvement. The system consists of a wave generation-absorbing system which initiates the exact desired hydrodynamic condition for computation; Boussinesq-Geen-Naghdi wave-current models at different approximation levels; a surf zone wave breaking and energy dissipating model; and a moving multi-shoreline algorithm. Different boundary conditions may be applied such as no normal flow and free tangential slip; periodic boundary; moving shoreline boundary; changing water surface atmosphere pressure; air-water stress; bottom frictional stress. Fundamentally rotational nonlinear phase resolving wave-current models are developed at different approximation levels for coastal ocean dynamics such as highly dispersive nonlinear waves, near-shore currents, and depth-varying velocity field, from deep water through the surf zone to the coast. Without imposing constraints on rotationality, we take the benefit of both Boussinesq scaling for water wave transformation and the Green-Naghdi weighted residual. The system uses polynomial basis functions to reconstruct velocity vertical profiles which are carried through the basic equations of motion keeping terms up to the desired Boussinesq scaling order, and solved in a weighted residual form. The models show rapid convergence to exact solutions for the linear dispersion relation, shoaling gradient, and orbital velocities; however, properties may be substantially improved for a given order of approximation using asymptotic rearrangements as is commonly used in other Boussinesq systems. This improvement is achieved using large numbers of degrees of freedom inherent in the expression of the polynomial basis functions to either match up to certain terms in a Taylor series, or minimize errors over a depth-wave length ratio, kh, range. Explicit coefficients are given at O(mu^2) and O(mu^4), while more generalized basis functions are given at higher order. The computational cost, after using Gaussian Elimination to reduce the number of coupled equations, is similar to that of normal Boussinesq theories, even though there are more unknown variables to solve than in normal Boussinesq models. The number of coupled equations can be reduced by multiplying proper coefficients and subtracting from each other. Therefore, the size of the matrix to be solved is in similar to that of normal Boussinesq models. Wave breaking and moving multi-shoreline theories are developed to extend the present model to the surf zone dynamics. Because the model is fundamentally rotational, it uses fewer ad-hoc assumptions than are found in many Boussinesq breaking wave systems. Eddy viscosity is used to describe both breaking dissipation and bottom stress effects with extra specified bottom frictional dissipation. The one-equation k-l turbulent kinetic energy equation coupled with rotational wave models are adopted to describe the eddy viscosity. This method results in a greatly improved velocity profile under breaking and comparable computational efficiency when compared to most other Boussinesq-type models. Still, more complete and accurate eddy viscosity model may be desired in future. For the simulation, the accurate generation and absorption of water waves in the phase-resolving model are important issues in describing nearshore processes. In order to introduce the target wave signals into the computation, a source function method for combined wave generation and absorption within modified sponge layers is presented here. This approach could be easily adapted to wide variety of systems, and does not required the solution of complex Green's functions but rather the simpler knowledge of solutions for free waves. These solutions may be linear or nonlinear, regular or irregular. Further, generated waves can be made arbitrarily accurate through the selection of various sponge layer coefficients. One-dimensional and two-dimensional testing cases are demonstrated; these cases show either good agreement with experiments or satisfactory prediction of real-world simulation. The wave transformation over varying depth, wave breaking, wave runup, undertow, rotational velocity field, nearshore currents, inundation, and other desired features are validated. Because the depth-varying velocity field is reconstructed by the definition of velocity expansion, the 1D and 2D simulations are actually quasi-2D and quasi-3D modelling. The rotational modelling ability has been fully demonstrated through a variety of test cases. The moving shoreline scheme works well, while a more accurate wet-dry algorithm may be an alternate in further development. Several wave and surge simulations were also performed; these simulations may be seen at http://www.youtube.com/userzhangyaond/videos . Future work should be extended to applications of the model to geophysical and engineering problems.