A direct numerical simulation is used to analyze the dynamics of a turbulent jet in a pycnocline. A horizontal cylindrical jet with a Gaussian mean velocity cross profile is considered. The velocity maximum lies within the pycnocline, and the initial jet diameter is set equal to the pycnocline width. In the numerical simulation, the Reynolds and Froude numbers in the jet are set close to their typical values measured in far wakes in laboratory experiments. The results suggest that two stages can be distinguished in flow development. At the initial stage (at times Nt < 10, where N is the characteristic buoyancy frequency), the vertical mean-velocity profile contracts, and vertical turbulent velocity fluctuations are collapsed, which is accompanied by the generation of internal waves. The mean velocity considerably exceeds the velocity in an unstratified jet at the same times. At the next stage, occurring at 10 < Nt < 100, the jet flow becomes self-similar. Large-scale vortices with vertical vorticity of alternating sign are formed at this stage. They are arranged in a staggered pattern near the streamwise (x) centerline of the jet. When Nt is sufficiently large, the streamwise vorticity field consists of oppositely signed transverse layers. The time dependences of the maximum mean velocity and the spanwise and vertical jet widths are described by asymptotics that agree well with experimental data. The asymptotics of the jet velocity and width are estimated theoretically as based on the property that the flow is self-similar and quasi-two-dimensional at late times. The estimates are in good agreement with the numerical results and experimental data.