Many mesoscale models (such as MM5) do not represent well the turbulence associated with tropopause jet streams. Shear instabilities in the free atmosphere generate energetic eddies on vertical scales O(10 m)-O(100 m), which lie in the subgrid scales of most mesoscale model applications. Hence, the results of any such modeling study, which crucially depends on physical and/or chemical processes in the upper troposphere-lower stratosphere region, will remain inconclusive until an appropriate parameterization scheme for the tropopausal clear air turbulence is incorporated. First, we will present results from 3D direct numerical simulations (DNS) of a nonuniformly stratified (with a doubling of buoyancy frequency (N) at the tropopause) model tropopause jet with 512 or 1024 vertical levels. The governing equations, used in the DNS, are the fully nonlinear, 3D incompressible, stratified Navier-Stokes equations, under the Boussinesq approximation, supplemented with horizontally homogeneous forcing terms. The vertical variability of important turbulence statistics and turbulence outer length scales are characterized. Detailed analyses of budgets of variances and vertical fluxes, quantifying all relative contributions from production, dissipation, pressure-strain and transport terms, for a quasi-equilibrium (i.e. nearly balanced budgets with negligible residual) DNS dataset will also be presented. Splitting of the transport (diffusion) terms revealed the existence of a two-way exchange of kinetic energy between the turbulence production zones, on either sides of the jet, and the innermost core of the jet, in addition to the one-way transport of kinetic energy from the production zones to the edges of the quasi-neutral mixing layer (in the mean potential temperature). The return transport from the centre of the jet to the production zones is accomplished through the turbulent (triple-moment) transport term in the equation for the vertical velocity variance. It is also found that scalings of several turbulent quantities (such as the flux Richardson number, Ri_f) with gradient Richardson number, Ri_g, exhibit multiple branches, as a consequence of the strongly nonhomogeneous stratification (due to the doubling of N at the model tropopause). This result implies that any unique Ri_f-Ri_g relationship, commonly used in many first-order turbulence closure models, becomes invalid under strong inhomogeneity in the stratification, such as that occuring near the tropopause. Based on the above promising quasi-balanced DNS dataset, we then develop a simple (homogeneous in horizontal and nonhomogeneous in vertical) second-order closure model which simulates all variances and vertical fluxes. The variances, fluxes, and budget terms simulated by the closure model match favorably with the vertical variability of most of the corresponding quantities in the original DNS dataset. It is also verified that the present closure model is capable of reproducing, qualitatively, the presence of multiple branches in the scaling of Ri_f with Ri_g.