Study of oceanic circulation and climate requires models which can simulate tracer eddy diffusion and ad-vection accurately. It is shown that the traditional Eulerian coordinates can introduce large artificial hori-zontal diffusivity/viscosity due to the incorrect alignment of the axis. Therefore, such models can smear sharp fronts and introduce other numerical artifacts. For simulation with relatively low resolution, large lateral diffusion was explicitly used in models;therefore, such numerical diffusion may not be a problem. However, with the increase of horizontal resolution, the artificial diffusivity/viscosity associated with hori-zontal advection in the commonly used Eulerian coordinates may become one of the most challenging ob-stacles for modeling the ocean circulation accurately. Isopycnal eddy diffusion (mixing) has been widely used in numerical models. The common wisdom is that mixing along isopycnal is energy free. However, a careful examination reveals that this is not the case. In fact, eddy diffusion can be conceptually separated into two steps:stirring and subscale diffusion. Due to the thermobaric effect, stirring, or exchanging water masses, along isopycnal surface is associated with the change of GPE in the mean state. This is a new type of instability, called the thermobaric instability. In addition, due to cabbeling subscale diffusion of water parcels always leads to the release of GPE. The release of GPE due to isopycnal stirring and subscale diffusion may lead to the thermobaric instability.