A quasi three-dimensional (3D) subgrid technique for the geodesic finite-difference time-domain (FDTD) algorithm is proposed to increase local resolution level for solving extremely low frequency (ELF) electromagnetic wave propagation problems in the Earth-ionosphere system. Basic idea of this technique is to first split the spherical triangular grid into smaller grids and then split the corresponding spherical hexagonal grid with smaller grids, after that the finite difference technique is applied to these grids to complete the subgrid process. In order to keep long term numerical stable, virtual field quantities are introduced and interpolation is used to calculate them. We have finally given an simulation example to prove correctness of the technique.