  A promising route towards fault-tolerant quantum error correction is the
concatenation of a Gottesman-Kitaev-Preskill (GKP) code with a qubit code.
Development of such concatenated codes requires simulation tools which
realistically model noise, while being able to simulate the dynamics of many
modes. However, so far, large-scale simulation tools for concatenated GKP codes
have been limited to idealized noise models and GKP code implementations. Here,
we introduce the Bosonic Pauli+ model (BP+), which can be simulated efficiently
for a large number of modes, while capturing the rich dynamics in the bosonic
multi-mode Hilbert space. We demonstrate the method by simulating a hybrid
surface code, where the data qubits are finite-energy GKP qubits stabilized
using the small-Big-small (sBs) protocol, and the syndrome qubits are standard
two-level systems. Using BP+, we present logical error rates of such an
implementation. Confidence in the accuracy of the method is gained by comparing
its predictions with full time evolution simulations for several relevant
quantum circuits. While developed specifically for GKP qubits stabilized using
the sBs protocol, the mathematical structure of BP+ is generic and may be
applicable also to the simulation of concatenations using other bosonic codes.
