The Modular Arbitrary-Order Ocean-Atmosphere Model: MAOOAM v1.0

This paper describes a reduced-order quasi-geostrophic coupled ocean-atmosphere model that allows for an arbitrary number of atmospheric and oceanic modes to be retained in the spectral decomposition. The modularity of this new model allows one to easily modify the model physics. Using this new model, coined the"Modular Arbitrary-Order Ocean-Atmosphere Model"(MAOOAM), we analyse the dependence of the model dynamics on the truncation level of the spectral expansion, and unveil spurious behaviour that may exist at low resolution by a comparison with the higher-resolution configurations. In particular, we assess the robustness of the coupled low-frequency variability when the number of modes is increased. An"optimal"configuration is proposed for which the ocean resolution is sufficiently high, while the total number of modes is small enough to allow for a tractable and extensive analysis of the dynamics.

for the atmospheric streamfunction fields ψ 1 a at 250 hPa and ψ 3 a at 750 hPa, and the vertical velocity ω = dp/dt, read The Coriolis parameter f is linearized around a value f 0 estimated at latitude φ 0 = 45 • N, f = f 0 + βy, with β = df /dy.
The parameters k d and k d quantify the friction between the two atmospheric layers and between the ocean and the atmosphere, respectively, and ∆p = 500 hPa is the pressure difference between the atmospheric layers.
The equation of motion for the streamfunction ψ o of the ocean layer reads (cf. Pierini, 2011) L R is the reduced Rossby deformation radius, ρ the density, h the depth, and r the friction at the bottom of the active ocean layer. The rightmost term represents the impact of the wind stress, and is modulated by the drag coefficient of the mechanical ocean-atmosphere coupling, d = C/(ρh).
The time evolution of the atmosphere and ocean temperatures T a and T o obeys the following equations: Here, γ a and γ o are the heat capacities of the atmosphere and the active ocean layer. ψ a = (ψ 1 a + ψ 3 a )/2 is the atmospheric barotropic streamfunction. λ is the heat transfer coefficient at the ocean-atmosphere interface, and σ is the static stability of the atmosphere, taken to be constant. The quartic terms represent the long-wave radiation fluxes between the ocean, the atmosphere, and outer space, with a the emissivity of the grey-body atmosphere and σ B the Stefan-Boltzmann constant. By decomposing the temperatures as T a = T 0 a + δT a and T o = T 0 o + δT o , the quartic terms are linearized around spatially uniform temperatures T 0 a and T 0 o , as detailed in Appendix B of Vannitsem et al. (2015). R a and R o are the short-wave radiation fluxes entering the atmosphere and the ocean that are also decomposed as R a = R 0 a + δR a and R o = R 0 o + δR o . The hydrostatic relation in pressure coordinates (∂Φ/∂p) = −1/ρ a where the geopotential height Φ i = f 0 ψ i a and the ideal gas relation p = ρ a RT a allow one to write the spatially dependent atmospheric temperature anomaly δT a = 2f 0 θ a /R, with θ a ≡ (ψ 1 a − ψ 3 a )/2 often referred to as the baroclinic streamfunction. R is the ideal gas constant. This can be used to eliminate the vertical velocity ω from Eqs. (1)-(2) and (4). This reduces the independent dynamical fields to the streamfunction fields ψ a and ψ o , and the spatially dependent temperatures δT a and δT o .
The prognostic equations for these four fields are then non-dimensionalized by dividing time by f −1 0 , distance by a characteristic length scale L, pressure by the difference ∆p, temperature by f 2 0 L 2 /R, and streamfunction by L 2 f 0 . A more detailed discussion of the model equations and their non-dimensionalization can be found in Vannitsem and De Cruz (2014) and Vannitsem et al. (2015).
All the parameters of the model equations used in the present work are listed in Table 1.

Expansion of the dynamical fields
In non-dimensionalized coordinates x = x/L and y = y/L, the domain is defined by (0 ≤ x ≤ 2π n , 0 ≤ y ≤ π), with n = 2L y /L x the aspect ratio between its meridional and zonal extents (see Table 1 for the value used here). The atmospheric flow is defined in a zonally periodic channel with no-flux boundary conditions in the meridional direction (∂ · a /∂x ≡ 0 at y = 0, π), whereas the oceanic flow is confined within an ocean basin by imposing no-flux boundaries in both the meridional (∂ · o /∂x ≡ 0 at y = 0, π) and zonal (∂ · o /∂y ≡ 0 at x = 0, 2π/n) directions. These boundary conditions limit the functions used in the Fourier expansion of the dynamical fields. With the proper normalization, the basis functions for the atmosphere must be of the following form, following the nomenclature of Cehelsky and Tung (1987): F K M,P (x , y ) = 2 cos(M nx ) sin(P y ) F L H,P (x , y ) = 2 sin(Hnx ) sin(P y ).
Analogously, the oceanic basis functions must be of the form with integer values of M , H, P , H o , and P o .
For example, the spectral truncation used by Charney and Straus (1980) can be specified as Eqs. (6) Ordering the basis functions as in Eqs. (6)-(8), along increasing values of M = H (o) and then P (o) , allows one to write the set The dynamical fields can then be written as the following truncated series expansions: Furthermore, the short-wave radiation or insolation is determined by δR a = C a F 1 ; δR o = C o F 1 . In Eq. (13), a term φ j is added to the oceanic basis function φ j (x , y ) in order to give it a vanishing spatial average. This is required to guarantee mass conservation in the ocean (Cessi and Primeau, 2001;McWilliams, 1977), but otherwise does not affect the dynamics. Indeed, it can be added a posteriori when plotting the field ψ o (x , y , t). This term is non-zero for odd P o and H o , The mass conservation is automatically satisfied for ψ a (x , y , t), as the spatial averages of the atmospheric basis functions Substituting the fields in Eqs.
(1)-(5) and projecting on the different basis functions yield 2(n a + n o ) ordinary differential equations (ODEs) for as many variables. Due to the linearization of the quartic temperature fields in Eqs. (4) and (5), these equations are at most bilinear (due to the advection term) in the variables ψ a,i , θ a,i , ψ o,j , and δT o,j , which will henceforth jointly be referred to as η i , the components of the state vector η.
To construct the dynamical equations of these variables, one has to compute the various projections or inner products with the basis functions, for which the following shorthand notation will be used: As described by Cehelsky and Tung (1987), the inner products for the atmosphere can be computed as purely algebraic formulae of the wave numbers P , M , and H. We reiterate these algebraic formulae in Sect. A1 of Appendix A and extend them with the formulae for both the ocean-atmosphere coupling terms and the ocean inner products in Sect. A2. The inner products can be represented as either two-dimensional or three-dimensional tensors, which are sparse but generally not diagonal.

Technical implementation
Substituting the fields by Eqs. (11)-(14) and calculating the coefficients using the expressions for the inner products as in Appendix A yields a set of N ≡ 2(n a + n o ) prognostic ordinary differential equations. These equations are at most bilinear in the variables η i (1 ≤ i ≤ N ) due to the linearization of the radiative terms around a reference temperature present in Eqs. (4)-(5). This system of ODEs can therefore be most generically expressed as the sum of a constant, a matrix multiplication, and a tensor contraction: This expression can be further simplified by adding a dummy variable that is identically equal to one: η 0 ≡ 1. This extra variable allows one to merge c i , m i,j , and t i,j,k into the tensor T i,j,k , in which the linear terms are represented by T i,j,0 and the constant term by T i,0,0 : The elements of the tensor T i,j,k are specified in Appendix B. Recasting the system of ordinary differential equations for η i in the form of a tensor contraction has certain advantages, as we will clarify below. The symmetry of Eq. (18) allows for a unique representation of T i,j,k , if it is taken to be upper triangular in the last two indices (T i,j,k ≡ 0 if j > k). Since T i,j,k is known to be sparse, it is stored using the coordinate list representation, i.e. a list of tuples (i, j, k, T i,j,k ). This representation renders the computation of the tendencies dη i /dt computationally very efficient as well as conveniently parallelizable.
Two implementations of MAOOAM are provided as a Supplement: one in Lua and one in Fortran. The Lua code is optimized for LuaJIT, a just-in-time compiler for Lua (Pall, 2015), and runs about 20 % slower than the Fortran version. By default, the model equations are numerically integrated using the Heun method. We have tested higher-accuracy methods, but these did not significantly change the results. The integration method can easily be changed; as an example, a fourth-order Runge-Kutta integrator is also included in the Lua implementation.

Derivation of Jacobian, tangent linear, and adjoint models
The form of Eq. (18) allows one to easily compute the Jacobian matrix of this system of ODEs. Indeed, denoting the right-hand side of Eq. (18) as dη i /dt = f i , the expression reduces to The differential form of the tangent linear (TL) model for a small perturbation δη TL of a trajectory η * is then simply (Kalnay, 2003) dδη TL To obtain the differential form of the adjoint model along the trajectory η * , the Jacobian is transposed to yield the following equations for the adjoint variable δη AD : This section details some key results obtained with the model for various levels of spectral truncation, with the set of parameter values given in Table 1. The parameter values for L, L R , λ, r, d, C o , C a , k d , and k d were selected as detailed in Vannitsem et al. (2015). The same value was chosen for k d and k d , as was done in Charney and Straus (1980); see also Vannitsem and De Cruz (2014). Unless otherwise stated, all the following results are obtained after first integrating the model for a transient period of 30 726.5 years. The model is subsequently integrated for another 92 179.6 years to obtain a sufficiently long trajectory from which good statistics can be extracted.
For the atmospheric part of the model, a previous study (Cehelsky and Tung, 1987), referred to as CT in the following, has shown that spurious chaos and a too large variability in the modes near the spectral cut-off could take place if the resolution is not high enough. These manifestations of spurious behaviour can lead to solutions that differ significantly from the solutions of the full partial differential equations (PDEs, here Eqs. 1-5). These findings lead us to the important question of convergence: to what degree has the solution of the truncated equations converged towards the solution of the PDEs? Although we do not have access to the latter, one can infer how the solutions are altered when the resolution is increased. Therefore, it cannot be asserted that convergence has been reached, and this point was also clearly stated in CT. However, we can reasonably suppose that when the solutions stabilize, they give an insight into the full dynamics.
This question is now addressed for the MAOOAM coupled atmosphere-ocean model.  Table 2. To alleviate the notation in the following, a model configuration denoted simply by H max x-P max y indicates that the resolution is the same in both components: The first panel of Fig. 1, with the atm. 2x-2y oc. 2x-4y resolution, shows the typical attractor geometry found in Vannitsem et al. (2015) and Vannitsem (2015) with a noisy, seemingly periodic orbit associated with the development of a large low-frequency signal. However, as the resolution is increased in both the ocean and atmosphere components, this structure destabilizes and we obtain more compact, noisy attractors in Figs. 1 and 2. The cause of this structural change is an interesting question in itself, which is worth exploring further in the future, as it is associated with the problem of structural stability of models, but is beyond the scope of the present work.
Regarding the question of convergence, the variability of the atmospheric variables becomes quite stable as the resolution increases beyond 6x-6y. Indeed, the bounds of the attractors on the vertical axis (ψ a,1 ) stabilize at this resolution. This result is in agreement with the findings of CT. On the other hand, the convergence is not yet reached for the oceanic variables whose variability is strongly affected by adding further modes as in the 7x-7y and 8x-8y resolutions.
The impact of the resolution on the solutions can also be examined by computing the variance of each variable of the barotropic and baroclinic streamfunctions, since these are associated with the kinetic and potential energy of the system (Yao, 1980). The presence of spurious behaviour can then be detected through substantial changes in this variability. The distributions of the total variance of the variables ψ a,i and ψ o,i are depicted in Figs. 3-6. The results show that the variance distribution does not change much beyond the 4x-4y resolution for the atmospheric component. However, for the oceanic component, the variance distribution is strongly modified when the resolution increases, and therefore one cannot conclude from Fig. 6 that some sort of convergence is reached at the 8x-8y resolution. To interpret this specific property, one must recall an important feature of two-dimensional quasi-geostrophic turbulence, namely the presence of a specific space scale, the Rhines scale, which delimits the two regimes associated with a wave-dominated dynamics and a turbulent dynamics. This space scale is given by where U represents the root-mean-square velocity of the energy-containing scales (Rhines, 1975;Vallis and Maltrud, 1993;Vallis, 2006) and β = df /dy is the meridional derivative of the Coriolis parameter f . If one takes the typical velocity of the order of a few metres per second and a few centimetres per second within the atmosphere and the ocean at large scales, the typical length scales will be of the order of 1000 and 100 km, respectively. Therefore the highest wave numbers necessary to resolve the wave-dominated part within the atmosphere and the ocean differ by a factor of 10. Coming back to our analysis, if this limit is reached for the atmosphere in our model at H, P = 4-5, we should suspect that a value of should be used for the ocean. This of course imposes strong constraints on our reduced-order model and would considerably limit its flexibility.
Let us now focus on the development of the LFV in these different model configurations, and let us define the geopotential height difference δz between the locations (π/n, π/4) and (π/n, 3π/4) of the model's non-dimensional domain: δz(t) = z(π/n, π/4, t) − z(π/n, 3π/4, t), where z is the geopotential height field, as in Vannitsem et al. (2015). The results shown in Figs. 7 and 8 indicate that the LFV, present for atm. 2x-2y oc. 2x-4y as in Vannitsem (2015), is a very weak signal at intermediate resolutions, but develops again when the number of modes is increased, as shown by the 1-year and 5-year running means. It suggests that the LFV previously found in low-resolution versions (see Fig. 7, panel atm. 2x-2y oc. 2x-4y) is a robust feature of the model. Moreover, at high resolutions this LFV is weaker than for the VDDG model version, but it seems closer to the actual dynamics found for the North Atlantic Oscillation (NAO) as discussed in Li and Wang (2003) and Stephenson et al. (2000).
The climatologies of the atmospheric barotropic streamfunction expressed in geopotential height further highlight the changes in the statistical properties of the model as a function of resolution. As shown in Figs. 9 and 10, the convergence is pretty fast toward an averaged zonal atmospheric circulation as the model resolution is increased. By contrast, the convergence for the oceanic streamfunction ψ o is less clear (Figs. 11 and 12), although a recurrent "global" double gyre is present for each resolution. As for the LFV, the topology of the gyres at high resolutions and their small-scale structures also seem to depend on whether H max , M max , H max o and P max o , P max are even or odd numbers.
The previous results point toward the important question of the optimal resolution of the oceanic component needed to get a sufficiently low-resolution model while keeping a dynamics with strong similarities to a very high-resolution model. To answer this question, we have performed some higher-resolution integrations, but on shorter time spans. The time span for each integration is given in Table 2.
The variance distributions of the oceanic streamfunction variables (see Fig. 13) have decreased at the spectral cut-off's edges compared to the distributions of the lower-resolution model configurations shown in Fig. 5. However, this decrease is not sufficient, and apparently spurious effects are still present. For instance, the decay is not identical in both directions, with a slower decay rate as the zonal wave number H o increases. We can even notice a peak in the distribution around H 0 = H max 0 and P 0 = 2 for all these higher model resolutions. This indicates that in fact we are still far from a quantitatively representative solution in the ocean. It confirms that, as stated previously, a resolution of the order of the Rhines scale is needed to achieve a good convergence. For the ocean, it corresponds to a 100 km resolution which would then require roughly 2000 modes. Such a model will of course be very computationally expensive and cannot be considered a "reduced"-order model anymore.
However, the comparison between the atm. 5x-5y oc. 12x-12y model configuration and the 10x-10y or 12x-12y model configurations shows that the former displays a large-scale behaviour close to the latter two, but with a reduced complexity and computational cost. This similarity can be assessed by considering the climatologies of these higher-resolution runs displayed in Fig. 14 and by watching the corresponding videos (see below). We therefore believe that the atm. 5x-5y oc. 12x-12y model configuration is a good candidate when investigating more realistic dynamics than the one presented in VDDG. It must however be stressed that the VDDG model is still an important tool in this hierarchy of models since it already contains the basic mechanisms leading to low-frequency variability. In addition, the climatologies shown in Fig. 14  In these videos, a striking feature is the presence of a westward wave propagation within the ocean while the LFV is developing in the coupled system. This feature has been associated with the propagation of Rossby-like waves (Vannitsem, 2015).

Conclusions
A new reduced-order coupled ocean-atmosphere model is presented, extending the low-resolution versions previously published (Vannitsem and De Cruz, 2014;Vannitsem et al., 2015). It is referred to as MAOOAM, for Modular Arbitrary-Order Ocean-Atmosphere Model. This new model retains the main features of the previous versions but allows for the selection of an arbitrary resolution within the ocean and the atmosphere. Besides the potential utility of this new functionality for evaluating the impact of the number of modes on the dynamics (as has been done here), it opens the possibility of addressing several new questions in a very flexible way, such as the development of a consistent stochastic parameterization scheme through scales, the understanding of the predictability problem at multiple scales and the role of model error, or the implementation of a data-assimilation scheme for the coupled ocean-atmosphere system.
In the present work, we have studied the impact of the resolution on the model solution's dynamics, by investigating the properties of the attractors and the variance distributions in both the oceanic and atmospheric components. The conclusion that can be drawn is that the convergence of the atmospheric component of the system is quite fast (as noted in Cehelsky and Tung, 1987), with variance distributions decreasing rapidly as a function of scale. However, the convergence of the oceanic component is much slower. Consequently, none of the solutions presented so far have satisfactorily converged toward a dynamics that correctly reflects the wave-dominated regime of the coupled ocean-atmosphere system. This regime corresponds to a resolution associated with the Rhines scale (which for the ocean is equal to 100 km or, equivalently, to wave numbers of the order . This stresses the need for high-resolution oceanic models to correctly represent the full coupled dynamics. One coupled model configuration which could, however, be recommended so far is the atm. 5x-5y oc. 12x-12y configuration, which seems to display some robustness in the ocean climatology as compared to the full 10x-10y and 12x-12y configurations. This conclusion requires further investigation with even higher resolutions, together with the use of more advanced tools of analysis like the computation of the Lyapunov exponents as in Vannitsem and Lucarini (2016). These can be computed using the tangent linear model version for which an implementation is also provided. This will be the subject of a future investigation.
The robustness of the LFV pattern, one of the most interesting features of the model, has also been explored. As it turns out, a LFV is still present in a large portion of the model configurations explored (not in 2x-4y, 3x-3y, and 4x-4y), but a weaker LFV signal is found when high-resolution configurations are used. A dominant signal is found with a wide variety of periods ranging from 1 to 100 years, depending on the model configuration. A more detailed analysis of the underlying structure of the system's attractor is needed to clarify the origin of this diversity, for instance through a bifurcation analysis as in Vannitsem et al. (2015). Note that the VDDG model is still an important tool in this hierarchy of models, since it already contains the basic mechanisms leading to the LFV.
Another interesting finding is the change of structure of the climatologies of the ocean gyres when choosing even or odd wave Appendix A: Formulae to compute the inner products In the formulae of the inner products of the atmospheric modes, Cehelsky and Tung (1987) use the following helper functions: λ(r) = 0 (r even) or 1 (r odd), (A3) The same notation will be used in this appendix. In what follows, δ ij is the Kronecker delta, so that δ ij = 1 if i = j, and 0 otherwise. Likewise, the function δ(x) used in this appendix is defined as Using these functions, the various coefficients of the model are calculated, starting with the internal atmosphere coefficients.

A1 Atmospheric coefficients
In the following, we consider the ordering of the basis function used in Eqs. (11)-(14). For the sake of clarity, we add an extra informative upper index specifying the type of the atmospheric function in the definitions below. However, the inner products are completely defined by the lower indices alone. The atmospheric functions are thus noted: and the oceanic functions

A1.1 The a i,j coefficients
These coefficients correspond to the eigenvalues of the Laplacian operator acting on the spectral expansion basis functions: which are given for each case by These coefficients are needed to evaluate the contribution of the β-terms, and only involve the Kand L-type base functions.
We have that These coefficients are given by and the non-zero ones are given by where we have used the functions defined at the beginning of this appendix. All the other permutations can be obtained thanks to g α,β,γ i,j,k = −g β,α,γ j,i,k = g γ,α,β k,i,j = g β,γ,α j,k,i . (A22)

A1.4 The b i,j,k coefficients
These coefficients are given by Therefore we obtain

A1.5 The s i,j coefficients
These coefficients encode the inner products between the atmospheric and oceanic basis functions: which gives These coefficients are related to the forcing of the ocean on the atmosphere. They are given by the formula where the M j,j are given by the eigenvalues of the Laplacian operator acting on the oceanic basis functions (see next section).

A2.1 The K i,j coefficients
These coefficients are related to the forcing of the atmosphere on the ocean. They are given by

A2.2 The M i,j coefficients
These coefficients identify with the eigenvalues of the Laplacian acting on the oceanic basis functions: These coefficients are needed to evaluate the contribution of the β-terms and are given by .

A2.4 The O i,j,k coefficients
These coefficients are given by A2.5 The C i,j,k coefficients These coefficients are given by These coefficients are related to the short-wave radiative forcing of the ocean and are given by Appendix B: Definition of the tensor T i,j,k The system of non-dimensionalized ODEs for the model variables is encoded in the model tensor T i,j,k , of which the complete definition is given in this appendix. Tensor elements that are not listed below are equal to zero. To alleviate the notations, we use a shorthand notation for the indices of the different variables, Furthermore, we suppress the upper indices which indicate the atmospheric function types but are otherwise not needed to unambiguously specify the inner products.

B1 Atmosphere equations
The components of the tensor for the atmosphere streamfunction are given by The atmospheric temperature equations are determined by the tensor elements where we used the non-dimenionalized quantities