A numerical method in the frequency domain is developed for analyzing three-dimensional gratings using the concept of a double-periodic magnetodielectric layer. The method is based on the three-dimensional volume integral equations for the equivalent electric and magnetic polarization currents of the assumed periodic medium. The integral equations are solved by using the integral functionals related to the polarization current distributions and the technique of double Floquet-Fourier series expansion. Once the integral functionals are determined, the scattered fields outside the layer are calculated accordingly. The unit cell of the layer comprises several parallelepiped segments of materials characterized by the complex-valued relative permittivity and permeability of step function profiles. The arbitrary profiles of three-dimensional dielectric or metallic gratings can be flexibly modeled by adjusting the material parameters and sizes or locations of the parallelepiped segments in the unit cell. Numerical examples for various grating geometries and their comparisons with those presented in the literature demonstrate the accuracy and usefulness of the proposed method.