The brain is a high-dimensional directional network system as it consists of many regions as network nodes that exert influences on each other. The directional influence exerted by one region to another is referred to as directional connectivity. It is challenging to produce accurate and scientifically meaningful estimates of many subjects' whole-brain directional networks based on their functional magnetic resonance imaging (fMRI) data for several reasons. First, the numbers of possible patterns of the whole-brain directional brain network are enormous due to the many network nodes (brain regions). Second, the whole-brain directional network can vary across subjects. Third, analyzing many subjects' fMRI data can be computationally intensive. To address these challenges, we develop a new model for whole-brain directional networks using the principles of the brain's functional organization. Specifically, we assume modular brain networks which reflect functional specialization and integration of the brain. To accommodate the variation of brain networks across subjects, we allow each region to fall into different modules in different subjects' brain networks by using a mixed membership stochastic blockmodel as a prior. We develop a variational Bayes approach to estimate the model and simultaneously identify modules and directional connections in hundreds of subjects' whole-brain directional networks. We show that our method is computationally efficient and can successfully identify functionally specialized brain areas, such as the visual cortex and motor cortex. Our method also reveals complex functions of other regions. Overall, our new method brings new insights into the functional organization of the human brain.