By a nonstationary conformal mapping, the principal equations are rewritten in a surfacefollowing coordinate system and reduced to two simple evolutionary equations for the elevation and the velocity potential of the surface; Fourier expansion is used to approximate these equations. High accuracy of the method was confirmed by validation of the nonstationary model against known solutions and by comparison between the results obtained with different resolution in the horizontal. The method developed is applied to simulation of evolution of wave fields with different initial conditions. Numerical experiments with initially monochromatic waves with different steepness show that the model is able to simulate breaking conditions when the surface becomes a multi-valued function of the horizontal coordinate; an estimate of the critical initial wave height that divides between non-breaking and eventually breaking waves is obtained. Simulations of nonlinear evolution of a wave field represented initially by two modes with close wave numbers (amplitude modulation) and a wave field with a phase modulation both result in appearance of large and very steep waves, which also break if the initial amplitudes are sufficiently large. Statistical properties of multimode linear and nonlinear wave fields have been compared.