Interaction of shock waves with turbulent boundary layers can enhance the surface heat flux dramatically. Reynolds-averaged Navier-Stokes simulations based on constant turbulent Prandtl number often give grossly erroneous heat transfer predictions in SBLI flows. This is due to the fact that the underlying Morkovin's hypothesis breaks down in the presence of shock waves; thus, the turbulent Prandtl number cannot be assumed to be a constant. In this paper, we develop a new variable turbulent Prandtl number model based on linearized Rankine-Hugoniot conditions applied to shock-turbulence interaction. The turbulent Prandtl number is a function of the shock strength and we propose a shock function to identify the location and strength of shock waves. The shock function also simulates the post-shock relaxation of the turbulent heat flux, akin to that observed in canonical shock-turbulence interaction. The model is combined with the wellvalidated shock-unsteadiness k-ω model and is applied to the complex shock topology observed in oblique shock-turbulent boundary layer interactions. Comparison with experimental data shows significant improvement in the surface heat transfer rate in the interaction region, both for attached and separated SBLI cases. The shock function is also used to propose a robust form of the existing shock-unsteadiness k-ω model that simplifies the numerical implementation enormously.