This thesis contains the formulation, development and initial tests of a computer model for the prediction of fully three dimensional turbulent free surface flows typically found at localised areas of river systems. It is the intention that the model will be used to predict flow situations which are fully three dimensional. The model is, therefore, tested against a fully three dimensional test case of flow in a two-stage meandering channel. However, the model is not intended simply to be for computing flows in meandering river channels. Rather the model is intended to be used in a variety of problems which are outlined in the thesis. The Reynolds Averaged Navier-Stokes equations form the basis of the physical system. The Reynolds stresses are represented by two different stress-strain relationships: (1) a linear relationship and (2) a non-linear relationship. These relationships rely on an eddy viscosity and a turbulence time-scale which are calculated from two characterising turbulence quantities, a velocity squared scale, k, and an inverse length scale, . These quantities are computed from differential transport equations. Non-linear stress-strain relationships are relatively new and, it has been argued by their originators, require application to several different problems to fully ascertain their potential for future use. The author addresses this demand by applying them to two new problems. These are flow in a plenum chamber and open channel flow over a backward facing step. The equations are solved by an operator splitting method which, it is argued, allows for an accurate and realistic treatment of the troublesome advection terms at low spatial resolutions. This is thought to be essential since for three dimensional problems owing to computer time limitations achieving grid independent solutions with low order schemes is at present very difficult. The advantage of the present approach is demonstrated with reference to a simple one dimensional analogue.