We present a methodology based on proper orthogonal decomposition (POD). We have implemented the POD-based solver in the large eddy simulation code Nek5000 and used it to solve the advection-diffusion equation for temperature in cases where buoyancy is not present. POD allows for the identification of the most energetic modes of turbulence when applied to a sufficient set of snapshots generated through Nek5000. The Navier-Stokes equations are then reduced to a set of ordinary differential equations by Galerkin projection. The flow field is reconstructed and used to advect the temperature on longer time scales and potentially coarser grids. The methodology is validated and tested on two problems: two-dimensional flow past a cylinder and three-dimensional flow in T-junctions. For the latter case, the benchmark chosen corresponds to the experiments of Hirota et al., who performed particle image velocimetry on the flow in a counterflow T-junction. In both test problems the dynamics of the reduced-order model reproduce well the history of the projected modes when a sufficient number of equations are considered. The dynamics of flow evolution and the interaction of different modes are also studied in detail for the T-junction.