The $gamma$-ray loud narrow line Seyfert 1 ($gamma$-NLSy1) galaxy PKS 1502+036 ($z=0.409$) exhibited its first $gamma$-ray outburst on 2015 December 20. In the energy range of 0.1-300 GeV, the highest flux measured by {it Fermi}-Large Area Telescope is (3.90 $pm$ 1.52) $times$ 10$^{-6}$ ph cm$^{-2}$ s$^{-1}$, which is the highest $gamma$-ray flux ever detected from this object. The associated spectral shape is soft ($Gamma_{0.1-300~{rm GeV}}=2.57pm0.17$) and this corresponds to an isotropic $gamma$-ray luminosity of (1.2 $pm$ 0.6) $times$ 10$^{48}$ erg s$^{-1}$. We generate the broadband spectral energy distribution (SED) during the GeV flare and reproduce it using a one zone leptonic emission model. The optical-UV spectrum can be explained by a combination of synchrotron and the accretion disk emission, whereas, the X-ray to $gamma$-ray SED can be satisfactorily reproduced by inverse-Compton scattering of thermal photons originated from the torus. The derived SED parameters hint for the increase in the bulk Lorentz factor as a major cause of the flare and the location of the emission region is estimated as outside the broad line region but still inside torus. A comparison of the GeV flaring SED of PKS 1502+036 with that of two other $gamma$-NLSy1 galaxies, namely, 1H 0323+342 ($z=0.061$) and PMN J0948+0022 ($z=0.585$), and also with FSRQ 3C 279 ($z=0.536$) has led to the conclusion that the GeV flaring SEDs of $gamma$-NLSy1 galaxies resemble with FSRQs and a major fraction of their bolometric luminosity is emitted at $gamma$-ray energies.