An efficient method for calculating impedance matrix elements is proposed for analysis of microstrip structures with an arbitrary substrate thickness. Closed-form Green's functions are derived by applying the GPOF method to the remaining function after the extraction of the contributions of the surface wave pole, source dipole itself, and quasi-static (i.e.real images) from a spectral domain Green's function. When closed-form Green's functions are used in conjunction with rooftop-pulse subsectional basis functions and the razor testing function in an MoM with an MPIE formulation, the integrals appearing in the calculation procedure of the diagonal matrix elements are of two types. The first is x0n [e^(-jk0(x02 + y02 +a2)1/2)/(x02 + y02 +a2)1/2)]dx0dy0 (where n=0, 1) for the contribution of both the source dipole itself or real images where a=0 and complex images where a=complex constant, while the other is x0n H0(2)(kρp (x02 + y02)1/2)dx0dy0 for the contribution of the surface wave pole where kρp is a real pole due to the surface wave. Adopting a polar coordinate for the integral for both cases of n=0 and n=1 and performing analytical integrations for n=1 with respect to the variable x0 for both types not only removes the singularities but also drastically reduces the evaluation time for the numerical integration. In addition, the above numerical efficiency is also retained for the off-diagonal elements. To validate the proposed method, several numerical examples are presented.