In this work, an improved method for the efficient automatic simulation of optical band shapes and resonance Raman (rR) intensities within the "independent mode displaced harmonic oscillator" is described. Despite the relative simplicity of this model, it is able to account for the intensity distribution in absorption (ABS), fluorescence, and rR spectra corresponding to strongly dipole allowed electronic transitions with high accuracy. In order to include temperature-induced effects, we propose a simple extension of the time dependent wavepacket formalism developed by Heller which enables one to derive analytical expressions for the intensities of hot bands in ABS and rR spectra from the dependence of the wavepacket evolution on its initial coordinate. We have also greatly optimized the computational procedures for numerical integration of complicated oscillating integrals. This is important for efficient simulations of higher-order rR spectra and excitation profiles, as well as for the fitting of experimental spectra of large molecules. In particular, the multimode damping mechanism is taken into account for efficient reduction of the upper time limit in the numerical integration. Excited state energy gradient as well as excited state geometry optimization calculations are employed in order to determine excited state dimensionless normal coordinate displacements. The gradient techniques are highly cost-effective provided that analytical excited state derivatives with respect to nuclear displacements are available. Through comparison with experimental spectra of some representative molecules, we illustrate that the gradient techniques can even outperform the geometry optimization method if the harmonic approximation becomes inadequate.