The nonlinear propagation of picosecond or femtosecond optical pulses in multimode fiber amplifiers underlies a variety of intriguing physical phenomena as well as the potential for scaling sources of ultrashort pulses to higher powers. However, existing theoretical models of ultrashort-pulse amplification do not include some critical processes, and, as a result, they fail to capture basic features of experiments. We introduce a numerical model that combines steady-state rate equations with the unidirectional pulse propagation equation, incorporating dispersion, Kerr and Raman nonlinearities, and gain/loss-spectral effects in a mode-resolved treatment that is computationally efficient. This model allows investigation of spatiotemporal processes that are strongly affected by gain dynamics. Its capabilities are illustrated through examinations of amplification in few-mode gain fiber, multimode nonlinear amplification, and beam cleaning in a multimode fiber amplifier.