  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.57\pm0.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.
