The realism of large scale numerical ocean models has improved dra- matically in recent years, in part because modern computers permit a more faithful representation of the differential equations by their algebraic analogs. Equally significant, if not more so, has been the improved under- standing of physical processes on space and time scales smaller than those that can be represented in such models. Today, some of the most challeng- ing issues remaining in ocean modeling are associated with parameterizing the effects of these high-frequency, small-space scale processes. Accurate parameterizations are especially needed in long term integrations of coarse resolution ocean models that are designed to understand the ocean vari- ability within the climate system on seasonal to decadal time scales. Traditionally, parameterizations of subgrid-scale, high-frequency mo- tions in ocean modeling have been based on simple formulations, such as the Reynolds decomposition with constant diffusivity values. Until recently, modelers were concerned with first order issues such as a correct represen- tation of the basic features of the ocean circulation. As the numerical simu- lations become better and less dependent on the discretization choices, the focus is turning to the physics of the needed parameterizations and their numerical implementation. At the present time, the success of any large scale numerical simulation is directly dependent upon the choices that are made for the parameterization of various subgrid processes.