Subgrid modeling for magnetohydrodynamic (MHD) turbulent shear flows is studied with the aid of a statistical result from a multiscale direct‐interaction theory. The ensemble‐mean MHD correlation function is combined with a filtering procedure to give a subgrid‐scale MHD stress. This result leads to a subgrid model of Smagorinsky type, which is applicable to large eddy simulations of earth’s magnetic dynamo and reversed‐field pinches of plasma.