Freshwater flux used as a natural boundary condition for the salinity balance is applied to a primitive equation model of the oceanic general circulation. Instead of the relaxation condition or the virtual salt flux boundary conditions used in many existing models, the real freshwater flux across the upper surface is specified as the vertical velocity boundary condition for the continuity equation, and the salinity flux is set to identically zero at the sea surface. Numerical experiments show that a model with the natural boundary conditions runs smoothly. Much important physics involving the freshwater flux emerge from the new model. The barotropic Goldsbrough–Stommel gyres driven by the precipitation and evaporation, which were excluded in the previous numerical models, are reproduced. In addition, the model's results reveal extremely complex structure of the three-dimensional circulation driven by the freshwater flux. In fact, a relatively small amount of freshwater flux drives very strong meridional and zonal cells and baroclinic gyres, which are 100 times stronger than the driving freshwater flux. Most importantly, the model provides an accurate description of the meridional salt fluxes and their roles in setting up the thermohaline circulation. It is suggested that, with or without the rigid-lid approximation, the real freshwater flux can be used as the upper boundary condition in oceanic general circulation models, including the mixed-layer models, the ice–ocean coupling models, and atmosphere–ocean coupling models.