114 second order on the ? dimension. Temporal discretization is achieved through a second order implicit factored scheme for nonlinear terms and a Crank–Nicholson scheme for linear terms. Phase errors are minimized because all terms in the numerical scheme are centered at the same instant, t = (n + 1/2?t). Phase errors are prone to occur in numerical schemes in which all terms are not centered in the same instant. The 3-D + 2-DH computations are entirely coupled and thence performed at every time step. The 3-D model uses the free surface elevation computed in the 2-DH model, and has a matching condition: the vertically averaged 3-D velocity profiles need to match the 2-DH velocities. When coupled with the 3-D model, the bottom stress in the 2-DH model depends on the 3-D velocity profiles. In simulations in which only the 2-DH model is required, the bottom stress is computed by the usual quadratic law. The 2-DH solving engine is based on the generalized minimal residual method iterative solver. The 3-D solution is implicit along the vertical ? direction and explicit along ? planes, thence it becomes numerically a 1-D problem for every water column represented by a mesh node. The 1-D problems are solved in a finite di?erence scheme, using double sweep techniques. An optional numeric-analytic solution for the vertical profiles is also available. This numeric- analytic option renders very good results in regions where 3-D advection is of less importance. Due to the iterative method of the solver, in general when modelling bays or coastal waters the best gain, i.e. the ratio between real time and simulation time, is obtained with an average Courant number around 5. It is quite common though to have a max Courant > 20 in a given mesh. Depending on the size of the mesh and scales of the problem, it is customary in SisBAHIA to have gains over 120 for coastal area models running on a PC powered by processors above 3 GHz. That is, 2 months of real time 3-D circulation in a coastal area domain would be simulated in 12 hours. Along open boundaries, elevations can be prescribed in many di?erent ways, including synthetic tides generated by given harmonic constants, and data measured or provided at discrete times. A di?erent value, and/or phase shift, can be given for each node along any open boundary segment. Imposition of radiation condition, in a way similar to the one adopted in the Princeton Ocean Model (POM), can also be used. The angle of in?ow along open boundaries can be automatically computed, as a function of the geometry of the domain, or can be enforced. The angle of out?ow can be free or enforced. Moreover, prescription of 2-DH current velocities along open boundaries assimilated from ocean circulation models can also be applied. In the case of land boundaries, one can prescribe either normal or imposed directional ?uxes or velocities, both of which can be constant or variable in time (e.g. a river discharge curve). Leaky boundaries are allowed, as are slip and no-slip boundaries, and the equivalent roughness along each boundary node can be prescribed. For the 3-D model, zero velocity is the bottom boundary condition, and the wind stress is prescribed as the free surface condition. Dynamic drying and ?ooding areas can be modelled with proper discretization, or can be simulated in a virtual way through leaky lateral boundaries. In computing bottom stresses, the amplitude of the equivalent bottom roughness can be specified for each bottom node, re?ecting the type of material (i.e. rock, sand, mud, vegetation, etc.). The computed friction coe?cients of the bottom vary dynamically in time and space. When performing turbulence modelling, a multi-scale model is employed. Horizontal subgrid scale turbulent stresses are based on filtering techniques, also known as Large Eddy Simulation