Numerical approximations to the dynamic equations are given which allow basin-size ocean circulation models formulated in isopycnic coordinates to accommodate variable bottom topography and irregular coastlines. Emphasis is placed on computational economy through the use of a split-explicit time integration scheme, on the proper formulation of the advection and Coriolis terms in the momentum equation in case of strongly varying layer thickness, and on the correct estimation of the horizontal pressure gradient force in grid boxes truncated by steep bottom slopes. The algorithms are tested in a series of two- and three-layer double-gyre experiments. In cases of steady forcing leading to a steady circulation, we are able to reproduce the expected motionless final state in coordinate layers that are below the direct influence of the wind forcing. This includes layers intersecting topographic obstacles. A long-term (25--30 year) vacillation tentatively associated with the outcropping of isopycnals along the edge of the cyclonic gyre in a steadily forced model is documented. The paper forms the basis for a subsequent study in which circulation features obtained in a realistic North Atlantic setting will be discussed. ¿ American Geophysical Union 1990