Classical approaches to model the turbulent heat flux rely on the gradient diffusion hypothesis which assumes an alignment between the temperature gradient and turbulent heat flux. This assumption is generally valid but breaks down in the vicinity of shocks. Particularly in shock-turbulence boundary layer interaction (STBLI), the anisotropy of the post-shock turbulence imparts an anisotropic heat flux to the wall that is not aligned with the mean thermal gradient. We derive a turbulent heat flux model that can be implemented into a Reynolds Stress Model (RSM). It models the turbulent heat flux with the information provided from the Reynolds stress components of RSM. The heat flux model is first validated against DNS of a planar shock turbulence interaction problem. Afterwards, the model is implemented into a full Navier-Stokes solver, CFL3D, in order to simulate a shock-boundary layer interaction problem.