We develop a hydrodynamic theory for flowing inhomogeneous polymer-nanorod composites (PNCs) coupling the Smoluchowski transport equation for the distribution function of the nanorod dispersed in a polymer matrix and the transport equation for the distribution of the polymer in the host matrix. The polymer molecule phase is modeled by bead-spring Rouse chains while the nanorod phase is modeled as semiflexible rods. The polymer-nanorod surface contact interaction and the conformational entropy of semiflexible nanorods are incorporated, resulting in a coupled system of nonlinear, nonlocal Smoluchowski equations for the polymer and nanorod. We then implement the theory to infer rheological properties and predict mesoscale morphologies in fully coupled plane shear flows. Our numerical study focuses on the mesoscale morphology development with respect to the surface contact interaction due to the pretreated surface properties of the nanorods, extending our studies on monodomain polymer-nanorod composites [16]. We find that surface contact interaction dominates the mesoscopic morphology and thereby corresponding rheological properties. When the nanorod favors parallel alignment with the polymer in the host matrix, the only globally stable state is the flow-aligning steady state. When the nanorod prefers to align orthogonally to the polymer in the matrix, however, spatially inhomogeneous structures, time-dependent homogeneous structures, and various spatial-temporal structures emerge in different regimes of the model parameter space and versus strength of the bulk imposed shear. Effective rheological features of the inhomogeneous morphologies are also predicted by the theory.