Starting from the fully-compressible Euler equations, a two-way-coupled system governing the long-wave behaviour of thin layers (with respect to the radius of Earth) representing the ocean and atmosphere, under an isentropic constraint, was derived. This approach incorporates bathymetry and topographic features as well as three-dimensional atmospheric non-uniformities through their depth-average over a spherical shell. Linear analysis of the obtained system yields two pairs of gravito-acoustic waves which are found to be representative of the fast-travelling atmospheric wave (with a propagation speed mainly governed by the atmospheric-layer-averaged speed of sound) and the slower-travelling gravity waves in the ocean (with a propagation speed mainly governed by local water depth). Remarkably, the 'Proudman resonance', observed in the forced shallow-water equation framework and invoked to justify, in part, observed large wave-heights, vanishes in favour of a continuous transition past the critical water depth, occurring when the two wave propagation speeds are closest. Two-dimensional non-linear global simulations were performed, using atmospheric conditions on the day, showcasing the predictive ability of the model. Local maxima of water-height disturbance in the farfield from the volcano, linked to the atmospheric wave deformation over time, are observed, emphasising the importance of the atmospheric-layer modelling and two-way coupling for any daylong predictions. An efficient implementation of the modelling strategy was carried out in the open source computational framework dNami to demonstrate the ability to perform faster-than-real-time simulations despite the additional equations in the governing system. Future work would see the strategy extended to incorporate additional layers and physics e.g. ocean and atmosphere stratification, interaction with the upper atmosphere.