  Multistabilities were found in the ocean-atmosphere flow, in a reduced order
ocean-atmosphere coupled model, when the non-linear temperature equations were
solved numerically. In this paper we explain how the full non-linear
Stefan-Bolzmann law was numerically implemented, and the resulting change to
the system dynamics compared to the original model where these terms were
linearised. Multiple stable solutions were found that display distinct
ocean-atmosphere flows, as well as different Lyapunov stability properties. In
addition, distinct Low Frequency Variability (LFV) behaviour was observed in
stable attractors. We investigated the impact on these solutions of changing
the magnitude of the ocean-atmospheric coupling, as well as the atmospheric
emissivity to simulate an increasing green-house effect. Where multistabilities
exist for fixed parameters, the possibility for tipping between solutions was
investigated, but tipping did not occur in this version of the model where
there is a constant solar forcing. This study was undertaken using a
reduced-order quasi-geostrophic ocean-atmosphere model, consisting of two
atmosphere layers, and one ocean layer, implemented in the Python programming
language.
