In this work we present a new model of BGK type for polyatomic gases. The model incorporates the different relaxation rates of translational, rotational and/or vibrational modes characterizing polyatomic molecules using a BGK-type equation, and additional relaxation equations for the temperatures associated to each internal energy mode. We construct an efficient numerical scheme which is implicit in the relaxation terms, and test the model and the scheme on several problems, confirming the Asymptotic Preserving properties of the scheme, and comparing the results provided by the model with experimental and DSMC simulations, carried out on the full Boltzmann polyatomic equation.