A three-dimensional primitive equation model Océan Parallélisé (OPA) was coupled to the biogeochemical model Pelagic Interaction Scheme for Carbon and Ecosystem Studies to simulate the ocean circulation and the marine biological productivity through the biogeochemical cycles of carbon and the main nutrients (P, N, Si, Fe). We focus on surface phytoplankton dynamics in the Indian Ocean extending from 30°S to 30°N and from 30°E to 120°E. The seasonal cycle of phytoplankton over the Indian Ocean is generally characterized by two blooms, one during the summer monsoon, the other one during the winter monsoon. Based on the method proposed by Lévy et al. (2007), different biogeochemical provinces can be defined during the summer and winter monsoons. The model performed relatively well by simulating the main features of the cumulated increase in chlorophyll, and the time of the bloom onsets is consistent with data. It also reproduced quite well the main biogeochemical provinces in good agreement with data in most of the Arabian Sea (except in the central part), the Bay of Bengal, and in the convergence zone south of the equator. The analysis of the modeled biogeochemical processes has shown that during the blooms onset periods, the most limiting nutrient was nitrogen except some areas around India and the eastern part of the Bay of Bengal where the ecosystem tends toward silicate limitation. These limitations change during the blooms development. The model also highlighted a variety of the critical physical processes (horizontal and vertical advection, turbulent diffusion, mixed layer depth) involved in each biogeochemical province bloom dynamics.