A new formulation of the wellbore boundary condition in multiphase reservoir simulation is suggested, allowing automatic determination of the direction of the flow at any time instant. It is particularly useful in coupled transient modeling of the liquid loading phenomenon in gas wells. The traditional Inflow Performance Relationship (IPR) uses a common well-block productivity index, but determines the inflow/outflow of the individual phases through a wellbore connection in an otherwise detached manner. The new approach takes into account that all phases should flow in the same direction and should have a common zero flow condition. The restriction leads to the definition of a new state variable: the multi-phase zero-flow pressure (MPZFP, P). According to the new formulation, the flow direction at the connection is determined by the sign of the difference between this grid-block variable and the wellbore flowing pressure. Some advantages of the new formulation are: the flow is always co-current through a single connection (but multiple connections can lead to cross-flow); the meaning of well-block productivity index is preserved; the phase composition of the connection stream is determined by the upstream condition; when wellbore pressure is far from the MPZFP, the traditional (detached) IPRs are recovered; and - last but not least - previous numerical convergence problems are completely eliminated during coupled dynamic simulation of the wellbore/reservoir system in the presence of significant capillary pressure. Examples illustrate the application of the new wellbore boundary condition in coupled modeling of the complex phenomena under liquid loading conditions.