The balance between the need of improving horizontal resolution in simulating local small-scale ocean processes and computational costs makes it desirable to refine model mesh locally. A three-dimensional, two-way nesting unstructured quadrilateral grid, primitive equations, finite-differencing, estuarine and coastal ocean model is developed for multi-scale modeling. Because the model grid is capable of multi-area nesting and multi-level refinement at each subdomain, the model is highly compatible with simulations involved in complex topography and studies of local small-scale ocean processes. The two-way information exchange is achieved by a virtual grid method, and its basic idea is to implement numerical integrations of variables at nesting interfaces with the support of virtual grid variables, which are interpolated or updated from actual grid variables. The model is novel for two interpolation schemes: the high-order spatial interpolation at the middle temporal level (HSIMT) parabolic interpolation scheme and HSIMT advection-equivalent interpolation scheme, and they have high-order accuracy and good consistency with the advection scheme applied to solving the tracer equations. The conservation of both volume and tracer contents is ensured via a flux correction algorithm. The two original interpolation schemes are examined in an ideal salinity advection experiment in the peak preservation skill, stability, and conservation properties. A realistic application to the Deep Waterway Project area in the Changjiang Estuary showed that the nested grid model can reproduce the hydrodynamic processes at the observed sites successfully while it failed to maintain the performance with the structured grid model in simulating the variance of salinity, for which the enforced conservation had primary responsibility. The HSIMT parabolic interpolation scheme was distinguished from other schemes for its outstanding performances in simulating salinity.