An improved three-field gyrofluid model is proposed to numerically simulate ion-scale turbulence in tokamak plasmas, which includes the nonlinear evolution of perturbed electrostatic potential, parallel ion velocity and ion pressure with adiabatic electron response. It is benchmarked through advancing a gyrofluid toroidal global (GFT_G) code as well as the local version (GFT_L), with the emphasis of the collisionless damping of zonal flows. The nonlinear equations are solved by using Fourier decomposition in poloidal and toroidal directions and semi-implicit finite difference method along radial direction. The numerical implementation is briefly explained, especially on the periodic boundary condition in GFT_L version. As a numerical test and also practical application, the nonlinear excitation of geodesic acoustic mode (GAM), as well as its radial structure, is investigated in tokamak plasma turbulence.