OpenArray v1.0: a simple operator library for the decoupling of ocean modeling and parallel computing

Rapidly evolving computational techniques are making a large gap between scientific aspiration and code implementation in climate modeling. In this work, we design a simple computing library to bridge the gap and decouple the work of ocean modeling from parallel computing. This library provides 12 basic operators that feature user-friendly interfaces, effective programming, and implicit parallelism. Several state-of-the-art computing techniques, including computing graph and just-in-time compiling, are employed to parallelize the seemingly serial code and speed up the ocean models. These operator interfaces are designed using native Fortran programming language to smooth the learning curve. We further implement a highly readable and efficient ocean model that contains only 1860 lines of code but achieves a 91 % parallel efficiency in strong scaling and 99 % parallel efficiency in weak scaling with 4096 Intel CPU cores. This ocean model also exhibits excellent scalability on the heterogeneous Sunway TaihuLight supercomputer. This work presents a promising alternative tool for the development of ocean models.

1 Introduction

Many earth system models have been developed in the past several decades to improve the predictive understanding of the earth system (Bonan and Doney, 2018; Collins et al., 2018; Taylor et al., 2012). These models are becoming increasingly complicated, and the amount of code has expanded from a few thousand lines to tens of thousands or even millions of lines. In terms of software engineering, an increase in code causes the models to be more difficult to develop and maintain.

The complexity of these models mainly originates from three aspects. First, more model components and physical processes have been embedded into the earth system models, leading to a 10-fold increase in the amount of code (e.g., Alexander and Easterbrook, 2015). Second, some heterogeneous and advanced computing platforms (e.g., Lawrence et al., 2018) have been widely used by the climate modeling community, resulting in a 5-fold increase in the amount of code (e.g., Xu et al., 2015). Last, most of the model programs need to be rewritten due to the continual development of novel numerical methods and meshes. The promotion of novel numerical methods and technologies produced in the fields of computational mathematics and computer science have been limited in climate science because of the extremely heavy burden caused by program rewriting and migration.

Over the next few decades, tremendous computing capacities will be accompanied by more heterogeneous architectures which are equipped with two or more kinds of cores or processing elements (Shan, 2006), thus making a much more sophisticated computing environment for climate modelers than ever before (National Research Council, 2012). Clearly, transiting the current earth system models to the next generation of computing environments will be highly challenging and disruptive. Overall, complex codes in earth system models combined with rapidly evolving computational techniques create a very large gap between scientific aspiration and code implementation in the climate modeling community.

To reduce the complexity of earth system models and bridge this gap, a universal and productive computing library is a promising solution. Through establishing an implicit parallel and platform-independent computing library, the complex models can be simplified and will no longer need explicit parallelism and transiting, thus effectively decoupling the development of ocean models from complicated parallel computing techniques and diverse heterogeneous computing platforms.

Many efforts have been made to address the complexity of parallel programming for numerical simulations, such as operator overloading, source-to-source translator, and domain-specific language (DSL). Operator overloading supports the customized data type and provides simple operators and function interfaces to implement the model algorithm. This technique is widely used because the implementation is straightforward and easy to understand (Corliss and Griewank, 1994; Walther et al., 2003). However, it is prone to work inefficiently because overloading execution induces numerous unnecessary intermediate variables, consuming valuable memory bandwidth resources. Using a source-to-source translator offers another solution. As indicated by the name, this method converts one language, which is usually strictly constrained by self-defined rules, to another (Bae et al., 2013; Lidman et al., 2012). It requires tremendous work to develop and maintain a robust source-to-source compiler. Furthermore, DSLs can provide high-level abstraction interfaces that use mathematical notations similar to those used by domain scientists so that they can write much more concise and more straightforward code. Some outstanding DSLs, such as ATMOL (van Engelen, 2001), ICON DSL (Torres et al., 2013), STELLA (Gysi et al., 2015), and ATLAS (Deconinck et al., 2017), are used by the numerical model community. Although they seem to be source-to-source techniques, DSLs are newly defined languages and produce executable programs instead of target languages. Therefore, the new syntax makes it difficult for the modelers to master the DSLs. In addition, most DSLs are not yet supported by robust compilers due to their relatively short history. Most of the source-to-source translators and DSLs still do not support the rapidly evolving heterogeneous computing platforms, such as the Chinese Sunway TaihuLight supercomputer which is based on the homegrown Sunway heterogeneous many-core processors and located at the National Supercomputing Center in Wuxi.

Other methods such as COARRAY Fortran and CPP templates provide alternative ways. Using COARRAY Fortran, a modeler has to control the reading and writing operation of each image (Mellor-Crummey et al., 2009). In a sense, one has to manipulate the images in parallel instead of writing serial code. In term of CPP templates, it is usually suitable for small amounts of code and difficult for debugging (Porkoláb et al., 2007).

Inspired by the philosophy of operator overloading, source-to-source translating and DSLs, we integrated the advantages of these three methods into a simple computing library which is called OpenArray. The main contributions of OpenArray are as follows:

• Easy to use. The modelers can write simple operator expressions in Fortran to solve partial differential equations (PDEs). The entire program appears to be serial and the modelers do not need to know any parallel computing techniques. We summarized 12 basic generalized operators to support whole calculations in a particular class of ocean models which use the finite difference method and staggered grid.

• High efficiency. We adopt some advanced techniques, including intermediate computation graphing, asynchronous communication, kernel fusion, loop optimization, and vectorization, to decrease the consumption of memory bandwidth and improve efficiency. Performance of the programs implemented by OpenArray is similar to the original but manually optimized parallel program.

• Portability. Currently OpenArray supports both CPU and Sunway platforms. More platforms including GPU will be supported in the future. The complexity of cross-platform migration is moved from the models to OpenArray. The applications based on OpenArray can then be migrated seamlessly to the supported platforms.

Furthermore, we developed a numerical ocean model based on the Princeton Ocean Model (POM; Blumberg and Mellor, 1987) to test the capability and efficiency of OpenArray. The new model is called the Generalized Operator Model of the Ocean (GOMO). Because the parallel computing details are completely hidden, GOMO consists of only 1860 lines of Fortran code and is more easily understood and maintained than the original POM. Moreover, GOMO exhibits excellent scalability and portability on both central processing unit (CPU) and Sunway platforms.

The remainder of this paper is organized as follows. Section 2 introduces some concepts and presents the detailed mathematical descriptions of formulating the PDEs into operator expressions. Section 3 describes the detailed design and optimization techniques of OpenArray. The implementation of GOMO is described in Sect. 4. Section 5 evaluates the performances of OpenArray and GOMO. Finally, discussion and conclusion are given in Sects. 6 and 7, respectively.

2 Concepts of the array, operator, and abstract staggered grid

In this section, we introduce three important concepts in OpenArray: array, operator, and abstract staggered grid to illustrate the design of OpenArray.

2.1 Array

To achieve this simplicity, we designed a derived data type, Array, which inspired our project name, OpenArray. The new Array data type comprises a series of information, including a 3-dimensional (3-D) array to store data, a pointer to the computational grid, a Message Passing Interface (MPI) communicator, the size of the halo region, and other information about the data distribution. All the information is used to manipulate the Array as an object to simplify the parallel computing. In traditional ocean models, calculations for each grid point and the i, j, and k loops in the horizontal and vertical directions are unavoidable. The advantage of taking the Array as an object is the significant reduction in the number of loop operations in the models, making the code more intuitive and readable. When using the OpenArray library in a program, one can use type(Array) to declare new variables.

2.2 Operator

To illustrate the concept of an operator, we first take a 2-dimensional (2-D) continuous equation solving sea surface elevation as an example:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial DU}{\partial x}}+{\displaystyle \frac{\partial DV}{\partial y}}=\mathrm{0},\end{array}$$

where η is the surface elevation, U and V are the zonal and meridional velocities, and D is the depth of the fluid column. We choose the finite difference method and staggered Arakawa C grid scheme, which are adopted by most regional ocean models. In Arakawa C grid, D is calculated at the centers, U component is calculated at the left and right side of the variable D, and V component is calculated at the lower and upper side of the variable D (Fig. 1). Variables ($D,U,V$) located at different positions own different sets of gird increments. Taking the term $\frac{\partial DU}{\partial x}$ as an example, we firstly apply linear interpolation to obtain the D^′s value at U point represented by tmpD. Through a backward difference to the product of tmpD and U, the discrete expression of $\frac{\partial DU}{\partial x}$ can be obtained.

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{tmpD}(i+\mathrm{1},j)=& \mathrm{0.5}\cdot \left(D\right(i+\mathrm{1},j)+D(i,j\left)\right)\\ & \cdot U(i+\mathrm{1},j),\end{array}\end{array}$$

and

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\displaystyle \frac{\partial DU}{\partial x}}={\displaystyle \frac{\mathrm{tmpD}\left(i+\mathrm{1},j\right)-\mathrm{tmpD}\left(i,j\right)}{\mathrm{d}x{\left(i,j\right)}^{\ast }}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot \mathrm{d}x(i,j{)}^{*}}}\cdot \left(\left(D(i+\mathrm{1},j)+D(i,j)\right)\cdot U(i+\mathrm{1},j)\right.\\ & \left.-\left(D(i,j)+D(i-\mathrm{1},j)\right)\cdot U(i,j)\right),\end{array}\end{array}$$

where $\mathrm{d}x{\left(i,j\right)}^{\ast }=\mathrm{0.5}\cdot \left(\mathrm{d}x\right(i,j)+\mathrm{d}x(i-\mathrm{1},j\left)\right)$.

Figure 1 Arrangement of variables in the staggered Arakawa C grid.

In this way, the above continuous equation can be discretized into the following form.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{{\mathit{\eta }}_{t+\mathrm{1}}\left(i,j\right)-{\mathit{\eta }}_{t-\mathrm{1}}\left(i,j\right)}{\mathrm{2}\cdot \mathrm{d}t}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot \mathrm{d}x{\left(i,j\right)}^{*}}}\cdot \left(\left(D\left(i+\mathrm{1},j\right)+D\left(i,j\right)\right)\cdot U\left(i+\mathrm{1},j\right)\right.\\ & \left.-\left(D\left(i,j\right)+D\left(i-\mathrm{1},j\right)\right)\cdot U\left(i,j\right)\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\cdot \mathrm{d}y{\left(i,j\right)}^{*}}}\left(\left(D\left(i,j+\mathrm{1}\right)+D\left(i,j\right)\right)\cdot V\left(i,j+\mathrm{1}\right)\right.\\ & \left.-\left(D\left(i,j\right)+D\left(i,j-\mathrm{1}\right)\right)\cdot V\left(i,j\right)\right)=\mathrm{0},\end{array}\end{array}$$

where $\mathrm{d}x{\left(i,j\right)}^{\ast }=\mathrm{0.5}\cdot \left(\mathrm{d}x\right(i,j)+\mathrm{d}x(i-\mathrm{1},j\left)\right),\mathrm{d}y{\left(i,j\right)}^{\ast }=\mathrm{0.5}\cdot \left(\mathrm{d}y\right(i,j)+\mathrm{d}y(i,j-\mathrm{1}\left)\right)$, and subscripts η_t+1 and η_t−1 denote the surface elevations at the (t+1) time step and (t−1) time step. To simplify the discrete form, we introduce some notation for the differentiation (${\mathit{\delta }}_{\mathrm{f}}^x$, ${\mathit{\delta }}_{\mathrm{b}}^y$) and interpolation (${\stackrel{\mathrm{‾}}{\left(\right)}}_{\mathrm{f}}^x$, ${\stackrel{\mathrm{‾}}{\left(\right)}}_{\mathrm{b}}^y)$. The δ and overbar symbols define the differential operator and average operator, respectively. The subscript x or y denotes that the operation acts in the x or y direction, and the superscript f or b denotes that the approximation operation is forward or backward.

Table 1 lists the detailed definitions of the 12 basic operators. The term var denotes a 3-D model variable. All 12 operators for the finite difference calculations are named using three letters in the form [A|D][$X\left|Y\right|Z$][F|B]. The first letter contains two options, A or D, indicating an average or a differential operator. The second letter contains three options, X, Y, or Z, representing the direction of the operation. The last letter contains two options, F or B, representing forward or backward operation. The dx, dy, and dz are the distances between two adjacent grid points along the x, y, and z directions.

Table 1 Definitions of the 12 basic operators.

Using the basic operators, Eq. (4) is expressed as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{{\mathit{\eta }}_{t+\mathrm{1}}-{\mathit{\eta }}_{t-\mathrm{1}}}{\mathrm{2}\cdot \mathrm{d}t}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot U\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot V\right)=\mathrm{0}.\end{array}$$

Thus,

$$\begin{array}{}\text{(6)}& {\mathit{\eta }}_{t+\mathrm{1}}={\mathit{\eta }}_{t-\mathrm{1}}-\mathrm{2}\cdot \mathrm{d}t\cdot \left({\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot U\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot V\right)\right).\end{array}$$

Then, Eq. (6) can be easily translated into a line of code using operators (the bottom left part in Fig. 2). Compared with the pseudo-codes (the right part), the corresponding implementation by operators is more straightforward and more consistent with the equations.

Figure 2 Implementation of Eq. (6) by basic operators. The elf and elb variables are the surface elevations at times (t+1) and (t−1), respectively.

Next, we will use the operators in shallow water equations, which are more complicated than those in the previous case. Assuming that the flow is in hydrostatic balance and that the density and viscosity coefficients are constant, and neglecting the molecular friction, the shallow water equations are

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial DU}{\partial x}}+{\displaystyle \frac{\partial DV}{\partial y}}=\mathrm{0},\text{(8)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial DU}{\partial t}}+{\displaystyle \frac{\partial DUU}{\partial x}}+{\displaystyle \frac{\partial DVU}{\partial y}}-fVD=-gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}\\ & +\mathit{\mu }D\left({\displaystyle \frac{{\partial }^{\mathrm{2}}U}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}U}{\partial y^{\mathrm{2}}}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial DV}{\partial t}}+{\displaystyle \frac{\partial DUV}{\partial x}}+{\displaystyle \frac{\partial DVV}{\partial y}}+fUD=-gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial y}}\\ & +\mathit{\mu }D\left({\displaystyle \frac{{\partial }^{\mathrm{2}}V}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}V}{\partial y^{\mathrm{2}}}}\right),\end{array}\end{array}$$

where f is the Coriolis parameter, g is the gravitational acceleration, and μ is the coefficient of kinematic viscosity. Using the Arakawa C grid and leapfrog time difference scheme, the discrete forms represented by operators are shown in Eqs. (10)–(12).

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle \frac{{\mathit{\eta }}_{t+\mathrm{1}}-{\mathit{\eta }}_{t-\mathrm{1}}}{\mathrm{2}\cdot \mathrm{d}t}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot U\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot V\right)=\mathrm{0},\text{(11)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{D_{t+\mathrm{1}}{U}_{t+\mathrm{1}}-D_{t-\mathrm{1}}{U}_{t-\mathrm{1}}}{\mathrm{2}\cdot \mathrm{d}t}}+{\mathit{\delta }}_{\mathrm{b}}^x\left({\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot U}}_{\mathrm{f}}^x\cdot {\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^x\right)\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot V}}_{\mathrm{b}}^x\cdot {\stackrel{\mathrm{‾}}{U}}_{\mathrm{b}}^y\right)-{\stackrel{\mathrm{‾}}{f{\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^y\cdot D}}_{\mathrm{b}}^x=-g\cdot {\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot {\mathit{\delta }}_{\mathrm{b}}^x\left(\mathit{\eta }\right)\\ & +\mathit{\mu }\cdot {\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot \left({\mathit{\delta }}_{\mathrm{b}}^x\left({\mathit{\delta }}_{\mathrm{f}}^x\left(U_{t-\mathrm{1}}\right)\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\mathit{\delta }}_{\mathrm{b}}^y\left(U_{t-\mathrm{1}}\right)\right)\right),\end{array}\text{(12)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{D_{t+\mathrm{1}}{V}_{t+\mathrm{1}}-D_{t-\mathrm{1}}{V}_{t-\mathrm{1}}}{\mathrm{2}\cdot \mathrm{d}t}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot U}}_{\mathrm{b}}^y\cdot {\stackrel{\mathrm{‾}}{V}}_{\mathrm{b}}^x\right)\\ & +{\mathit{\delta }}_{\mathrm{b}}^y\left(\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot V}{}_{\mathrm{f}}^y\cdot {\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^y\right)+{\stackrel{\mathrm{‾}}{f{\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^x\cdot D}}_{\mathrm{b}}^y=-g\cdot {\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot {\mathit{\delta }}_{\mathrm{b}}^y\left(\mathit{\eta }\right)\\ & +\mathit{\mu }\cdot {\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot \left({\mathit{\delta }}_{\mathrm{f}}^x\left({\mathit{\delta }}_{\mathrm{b}}^x\left(V_{t-\mathrm{1}}\right)\right)+{\mathit{\delta }}_{\mathrm{b}}^y\left({\mathit{\delta }}_{\mathrm{f}}^y\left(V_{t-\mathrm{1}}\right)\right)\right).\end{array}\end{array}$$

As the shallow water equations are solved, spatial average and differential operations are called repeatedly. Implementing these operations is troublesome and thus it is favorable to abstract these common operations from PDEs and encapsulate them into user-friendly, platform-independent, and implicit parallel operators. As shown in Fig. 3, we require only three lines of code to solve the shallow water equations. This more realistic case suggests that even more complex PDEs can be constructed and solved by following this elegant approach.

Figure 3 Implementation of the shallow water equations by basic operators. elf, el, and elb denote sea surface elevations at times (t+1), t, and (t−1), respectively. Uf, U and Ub denote the zonal velocity at times (t+1), t, and (t−1), respectively. Vf, V, and Vb denote the meridional velocity at times (t+1), t, and (t−1), respectively. Variable aam denotes the viscosity coefficient.

2.3 Abstract staggered grid

Most ocean models are implemented based on the staggered Arakawa grids (Arakawa and Lamb, 1981; Griffies et al., 2000). The variables in ocean models are allocated at different grid points. The calculations that use these variables are performed after several reasonable interpolations or differences. When we call the differential operations on a staggered grid, the difference value between adjacent points should be divided by the grid increment to obtain the final result. Setting the correct grid increment for modelers is troublesome work that is extremely prone to error, especially when the grid is nonuniform. Therefore, we propose an abstract staggered grid to support flexible switching of operator calculations among different staggered grids. When the grid information is provided at the initialization phase of OpenArray, a set of grid increments, including horizontal increments (dx(i,j), dy(i,j)), and vertical increment (dz(k)), will be combined with each corresponding physical variable through grid binding. Thus, the operators can implicitly set the correct grid increments for different Array variables, even if the grid is nonuniform.

As shown in Fig. 4, the cubes in the (a), (b), (c), and (d) panels are the minimum abstract grid accounting for 1∕8 of the volume of the cube in panel (e). The eight points of each cube are numbered sequentially from 0 to 7, and each point has a set of grid increments, i.e., dx, dy, and dz. For example, all the variables of an abstract Arakawa A grid are located at point 3. For the Arakawa B grid, the horizontal velocity Array (U, V) is located at point 0, the temperature (T), the salinity (S), and the depth (D) are located at point 3, and the vertical velocity Array (W) is located at point 7. For the Arakawa C grid, Array U is located at point 2 and Array V is located at point 1. In contrast, for the Arakawa D grid, Array U is located at point 1 and Array V is located at point 2.

Figure 4 The schematic diagram of the relative positions of the variables on the abstract staggered grid and the jumping procedures among the grid points.

When we call the average and differential operators mentioned in Table 1, e.g., on the abstract Arakawa C grid, the position of Array D is point 3, and the average AXB operator acting on Array D will change the position from point 3 to point 1. Since Array U is also allocated at point 1, the operation AXB(D)⋅U is allowed. In addition, the subsequent differential operator on Array AXB(D)⋅U will change the position of Array DXF(AXB(D)⋅U) from point 1 to point 3.

The jumping rules of different operators are given in Table 2. Due to the design of the abstract staggered grids, the jumping rules for the Arakawa A, B, C, and D grids are fixed. A change in the position of an array is determined only by the direction of a certain operator acting on that array.

Table 2 The jumping rules of an operator acting on an Array.

The grid type can be changed conveniently. For instance, if one would like to do so, only the following steps need to be taken. First the position information of each physical variable needs to be reset (shown in Fig. 4). Then the discrete form of each equation needs to be redesigned. We take the Eq. (1) switching from Arakawa C grid to Arakawa B grid as an example. The positions of the horizontal velocity Array U and Array V are changed to point 0; Array η and Array D stay the same. The discrete form is changed from Eq. (4) into Eq. (13), and the corresponding implementation by operators is changed from Eq. (6) into Eq. (14). By doing so, the transformation of grid types can be easily implemented.

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathit{\eta }}_{t+\mathrm{1}}\left(i,j\right)-{\mathit{\eta }}_{t-\mathrm{1}}\left(i,j\right)}{\mathrm{2}\cdot \mathrm{d}t}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\cdot \mathrm{d}x(i,j{)}^{*}}}\left(\right(D(i+\mathrm{1},j)+D(i,j))\cdot (U(i+\mathrm{1},j)\\ & +U(i+\mathrm{1},j+\mathrm{1}))-(D(i,j)+D(i-\mathrm{1},j))\\ & \cdot \left(U\right(i,j)+U(i,j+\mathrm{1}\left)\right))\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\cdot \mathrm{d}y(i,j{)}^{*}}}\left(\right(D(i,j+\mathrm{1})+D(i,j))\\ & \cdot \left(V\right(i,j+\mathrm{1})+V(i+\mathrm{1},j+\mathrm{1}\left)\right)\\ & -\left(D\right(i,j)+D(i,j-\mathrm{1}\left)\right)\\ & \cdot \left(V\right(i,j)+V(i+\mathrm{1},j\left)\right))=\mathrm{0},\end{array}\text{(14)}& {\displaystyle }{\mathit{\eta }}_{t+\mathrm{1}}={\mathit{\eta }}_{t-\mathrm{1}}-\mathrm{2}\cdot \mathrm{d}t\cdot \left({\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\cdot {\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^y\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\cdot {\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^x\right)\right).\end{array}$$

The position information and jumping rules are used to implicitly check whether the discrete form of an equation is correct. The grid increments are hidden by all the differential operators, thus it makes the code simple and clean. In addition, since the rules are suitable for multiple staggered Arakawa grids, the modelers can flexibly switch the ocean model between different Arakawa grids. Notably, the users of OpenArray should input the correct positions of each array in the initialization phase. The value of the position is an input parameter when declaring an Array. An error will be reported if an operation is performed between misplaced points.

Although most of the existing ocean models use finite difference or finite volume methods on structured or semi-structured meshes (e.g., Blumberg and Mellor, 1987; Shchepetkin and McWilliams, 2005), there are still some ocean models using unstructured meshes (e.g., Chen et al., 2003; Korn, 2017), and even the spectral element method (e.g., Levin et al., 2000). In our current work, we design the basic operators only for finite difference and finite volume methods with structured grids. More customized operators for the other numerical methods and meshes will be implemented in our future work.

3 Design of OpenArray

Through the above operator notations in Table 1, ocean modelers can quickly convert the discrete PDEs into the corresponding operator expression forms. The main purpose of OpenArray is to make complex parallel programming transparent to the modelers. As illustrated in Fig. 5, we use a computation graph as an intermediate representation, meaning that the operator expression forms written in Fortran will be translated into a computation graph with a particular data structure. In addition, OpenArray will use the intermediate computation graph to analyze the dependency of the distributed data and produce the underlying parallel code. Finally, we use stable and mature compilers, such as the GNU Compiler Collection (GCC), Intel compiler (ICC), and Sunway compiler (SWACC), to generate the executable programs according to different back-end platforms. These four steps and some related techniques are described in detail in this section.

Figure 5 The workflow of OpenArray.

3.1 Operator expression

Although the basic generalized operators listed in Table 1 are only suitable to execute 1st-order difference, other high-order difference or even more complicated operations can be combined by these basic operators. For example, a 2nd-order difference operation can be expressed as ${\mathit{\delta }}_{\mathrm{f}}^x\left({\mathit{\delta }}_{\mathrm{b}}^x\right(\mathrm{var}\left)\right)$. Supposing the grid distance is uniform, the corresponding discrete form is $\left[\mathrm{var}\right(i+\mathrm{1},j,k)+\mathrm{var}(i-\mathrm{1},j,k)-\mathrm{2}\cdot \mathrm{var}(i,j,k\left)\right]/\mathrm{d}{x}^{\mathrm{2}}$. In addition, the central difference operation can be expressed as $({\mathit{\delta }}_{\mathrm{f}}^x\left(\mathrm{var}\right)+{\mathit{\delta }}_{\mathrm{b}}^x(\mathrm{var}\left)\right)/\mathrm{2}$ since the corresponding discrete form is $\left[\mathrm{var}\right(i+\mathrm{1},j,k)-\mathrm{var}(i-\mathrm{1},j,k\left)\right]/\mathrm{2}\mathrm{d}x$.

Using these operators to express the discrete PDE, the code and formula are very similar. We call this effect “the self-documenting code is the formula”. Figure 6 shows the one-to-one correspondence of each item in the code and the items in the sea surface elevation equation. The code is very easy to program and understand. Clearly, the basic operators and the combined operators greatly simplify the development and maintenance of ocean models. The complicated parallel and optimization techniques are hidden behind these operators. Modelers no longer need to care about details and can escape from the “parallelism swamp”, and can therefore concentrate on the scientific issues.

Figure 6 The effect of “the self-documenting code is the formula” illustrated by the sea surface elevation equation.

3.2 Intermediate computation graph

Considering the example mentioned in Fig. 6, if one needs to compute the term DXF(AXB(D)⋅u) with the traditional operator overloading method, one first computes AXB(D) and stores the result into a temporary array (named tmp1), and then executes (tmp1⋅u) and stores the result into a new array, tmp2. The last step is to compute DXF(tmp2) and store the result in a new array, tmp3. Numerous temporary arrays consume a considerable amount of memory, making the efficiency of operator overloading poor.

To solve this problem, we convert an operator expression form into a directed and acyclic graph, which consists of basic data and function nodes, to implement a so-called lazy expression evaluation (Bloss et al., 1988; Reynolds, 1999). Unlike the traditional operator overloading method, we overload all arithmetic functions to generate an intermediate computation graph rather than to obtain the result of each function. This method is widely used in deep-learning frameworks, e.g., TensorFlow (Abadi et al., 2016) and Theano (Bastien et al., 2012), to improve computing efficiency. Figure 7 shows the procedure of parsing the operator expression form of the sea level elevation equation into a computation graph. The input variables in the square boxes include the sea surface elevation (elb), the zonal velocity (u), the meridional velocity (v), and the depth (D). Variable dt2 is a constant equal to 2⋅dt. The final output is the sea surface elevation at the next time step (elf). The operators in the circles have been overloaded in OpenArray. In summary, all the operators provided by OpenArray are functions for the Array calculation, in which the “=” notation is the assignment function, the “–” notation is the subtraction function, the “*” notation is the multiplication function, the “+” notation is the addition function, DXF and DYF are the differential functions, and AXF and AYF are the average functions.

Figure 7 Parsing the operator expression form into the computation graph.

3.3 Code generation

Given a computation graph, we design a lightweight engine to generate the corresponding source code (Fig. 8). Each operator node in the computation graph is called a kernel. The sequence of all kernels in a graph is usually fused into a large kernel function. Therefore, the underlying engine schedules and executes the fused kernel once and obtains the final result directly without any auxiliary or temporary variables. Simultaneously, the scheduling overhead of the computation graph and the startup overhead of the basic kernels can be reduced.

Figure 8 The schematic diagram of kernel fusion.

Most of the scientific computational applications are limited by the memory bandwidth and cannot fully exploit the computing power of a processor. Fortunately, kernel fusion is an effective optimization method to improve memory locality. When two kernels need to process some data, their fusion holds shared data in the memory. Prior to the kernel fusion, the computation graph is analyzed to find the operator nodes that can be fused, and the analysis results are stored in several subgraphs. Users can access any individual subgraph by assigning the subgraph to an intermediate variable for diagnostic purposes. After being given a series of subgraphs, the underlying engine dynamically generates the corresponding kernel function in C$++$ using just-in-time (JIT) compilation techniques (Suganuma and Yasue, 2005). The JIT compiler used in OpenArray can fuse numbers of operators into a large compiled kernel. The benefit of fusing operators is to alleviate memory bandwidth limitations and improve performance compared with executing operators one by one. In order to generate a kernel function based on a subgraph, we first add the function header and variable definitions according to the name and type in the Array structure. And then we add the loop head through the dimension information. Finally, we perform a depth-first walk on the expression tree to convert data, operators, and assignment nodes into a complete expression including load variables, arithmetic operation, and equal symbol with C$++$ language.

Notably, the time to compile a single kernel function is short but practical applications usually need to be run for thousands of time steps and the overhead of generating and compiling the kernel functions for the computation graph is extremely high. Therefore, we generate a fusion kernel function only once for each subgraph, and put it into a function pool. Later, when facing the same computation subgraph, we fetch the corresponding fusion kernel function directly from the pool.

Since the arrays in OpenArray are distributed among different processing units, and the operator needs to use the data in the neighboring points, to ensure the correctness it is necessary to check the data consistency before fusion. The use of different data-splitting methods for distributed arrays can greatly affect computing performance. The current data-splitting method in OpenArray is the widely used block-based strategy. Solving PDEs on structured grids often divides the simulated domain into blocks that are distributed to different processing units. However, the differential and average operators always require their neighboring points to perform array computations. Clearly, ocean modelers have to frequently call corresponding functions to carefully control the communication of the local boundary region.

Therefore, we implemented a general boundary management module to implicitly maintain and update the local boundary information so that the modelers no longer need to address the message communication. The boundary management module uses asynchronous communication to update and maintain the data of the boundary region, which is useful for simultaneous computing and communication. These procedures of asynchronous communication are implicitly invoked when calling the basic kernel or the fused kernel to ensure that the parallel details are completely transparent to the modelers. For the global boundary conditions of the limited physical domains, the values at the physical border are always set to zero within the operators and operator expressions. In realistic cases, the global boundary conditions are set by a series of functions (e.g., radiation, wall) provided by OpenArray.

3.4 Portable program for different back-end platforms

With the help of dynamic code generation and JIT compilation technology, OpenArray can be migrated to different back-end platforms. Several basic libraries, including Boost C$++$ libraries and Armadillo library, are required. The JIT compilation module is based on low-level virtual machine (LLVM), thus theoretically the module can only be ported to platforms supporting LLVM. If LLVM is not supported, as on the Sunway platform, one can generate the fusion kernels in advance by running the ocean model on an X86 platform. If the target platform uses CPUs with acceleration cards, such as GPU clusters, it is necessary to add control statements in the CPU code, including data transmission, calculation, synchronous, and asynchronous statements. In addition, the accelerating solution should involve the selection of the best parameters, e.g., “blockDim” and “gridDim” on GPU platforms. In short, the code generation module of OpenArray also needs to be refactored to be able to generate codes for different back-end platforms. The application based on OpenArray can then be migrated seamlessly to the target platform. Currently, we have designed the corresponding source code generation module for Intel CPU and Sunway processors in OpenArray.

According to the TOP500 list released in November 2018, the Sunway TaihuLight is ranked third in the world, with a LINPACK benchmark rating of 93 petaflops provided by Sunway many-core processors (or Sunway CPUs). As shown in Fig. 9, every Sunway CPU includes 260 processing elements (or cores) that are divided into four core groups. Each core group consists of 64 computing processing elements (CPEs) and a management processing element (MPE) (Qiao et al., 2017). CPEs handle large-scale computing tasks and the MPE is responsible for the task scheduling and communication. The relationship between MPE and CPE is like that between CPU and many-core accelerator, except that they are fused into a single Sunway processor sharing a unified memory space. To make the most of the computing resources of the Sunway TaihuLight, we generate kernel functions for the MPE, which is responsible for the thread control, and CPE, which performs the computations. The kernel functions are fully optimized with several code optimization techniques (Pugh, 1991) such as loop tiling, loop aligning, single-instruction multiple-date (SIMD) vectorization, and function inline. In addition, due to the high memory access latency of CPEs, we accelerate data access by providing instructions for direct memory access in the kernel to transfer data between the main memory and local memory (Fu et al., 2017).

Figure 9 The MPE–CPE hybrid architecture of the Sunway processor. Every Sunway processor includes 4 core groups (CGs) connected by the network-on-chip (NoC). Each CG consists of a management processing element (MPE), 64 computing processing elements (CPEs) and a memory controller (MC). The Sunway processor uses the system interface (SI) to connect with outside devices.

4 Implementation of GOMO

In this section, we introduce how to implement a numerical ocean model using OpenArray. The most important step is to derive the primitive discrete governing equations in operator expression form, and then the following work is completed by OpenArray.

The fundamental equations of GOMO are derived from POM. GOMO features a bottom-following free-surface staggered Arakawa C grid. To effectively evolve the rapid surface fluctuations, GOMO uses the mode-splitting algorithm inherited from POM to address the fast-propagating surface gravity waves and slow-propagating internal waves in barotropic (external) and baroclinic (internal) modes, respectively. The details of the continuous governing equations, the corresponding operator expression form, and the descriptions of all the variables used in GOMO are listed in the appendices A–C, respectively.

Figure 10 shows the basic flow diagram of GOMO. At the beginning, we initialize OpenArray to make all operators suitable for GOMO. After loading the initial values and the model parameters, the distance information is input into the differential operators through grid binding. In the external mode, the main consumption is computing the 2-D sea surface elevation η and column-averaged velocity (U_A, V_A). In the internal mode, 3-D array computations predominate in order to calculate baroclinic motions ($U,V,W$), tracers ($T,S,\mathit{\rho }$), and the turbulence closure scheme (q²,q²l) (Mellor and Yamada, 1982), where $U,V,$ and W are the velocity fields in the x, y, and σ directions and $T,S,$ and ρ are the potential temperature, the salinity, and the density, respectively. $q^{\mathrm{2}}/\mathrm{2},$ and q²l∕2 are the turbulence kinetic energy and production of turbulence kinetic energy with turbulence length scale, respectively.

Figure 10 Flow diagram of GOMO.

Table 3 Comparing GOMO with several variations in the POM.

When the user dives into the GOMO code, the main time stepping loop in GOMO appears to run on a single processor. However, as described above, implicit parallelism is the most prominent feature of the program using OpenArray. The operators in OpenArray, not only the difference and average operators, but also the “+”, “–”, “*”, “/” and “=” operators in the Fortran code, are all overloaded for the special data structure “Array”. The seemingly serial Fortran code is implicitly converted to parallel C$++$ code by OpenArray, and the parallelization is hidden from the modelers.

Because the complicated parallel optimization and tuning processes are decoupled from the ocean modeling, we completely implemented GOMO based on OpenArray in only 4 weeks, whereas implementation may take several months or even longer when using the MPI or CUDA library.

In comparison with the existing POM and its multiple variations, including (to name a few) Stony Brook Parallel Ocean Model (sbPOM), mpiPOM, and POMgpu, GOMO has less code but is more powerful in terms of compatibility. As shown in Table 3, the serial version of POM (POM2k) contains 3521 lines of code. Models sbPOM and mpiPOM are parallelized using MPI, while POMgpu is based on MPI and CUDA-C. The codes of sbPOM, mpiPOM, and POMgpu are extended to 4801, 9680, and 30 443 lines, respectively. In contrast, the code of GOMO is decreased to 1860 lines. Moreover, GOMO completes the same function as the other approaches while using the least amount of code (Table 4), since the complexity has been transferred to OpenArray, which includes about 11 800 lines of codes.

Table 4 Comparison of the amount of code for different functions.

In addition, poor portability considerably restricts the use of advanced hardware in oceanography. With the advantages of OpenArray, GOMO is adaptable to different hardware architectures, such as the Sunway processor. The modelers do not need to modify any code when changing platforms, eliminating the heavy burden of transmitting code. As computing platforms become increasingly diverse and complex, GOMO becomes more powerful and attractive than the machine-dependent models.

5 Results

In this section, we first evaluate the basic performance of OpenArray using benchmark tests on a single CPU platform. After checking the correctness of GOMO through an ideal seamount test case, we use GOMO to further test the scalability and efficiency of OpenArray.

5.1 Benchmark testing

We choose two typical PDEs and their implementations from Rodinia v3.1, which is a benchmark suite for heterogeneous computing (Che et al., 2009), as the original version. For comparison, we re-implement these two PDEs using OpenArray. In addition, we added two other test cases. As shown in Table 5, the 2-D continuity equation is used to solve sea surface height, and its continuous form is shown in Eq. (1). The 2-D heat diffusion equation is a parabolic PDE that describes the distribution of heat over time in a given region. Hotspot is a thermal simulation used for estimating processor temperature on structured grids (Che et al., 2009; Huang et al., 2006). We tested one 2-D case (Hotspot2D) and one 3-D case (Hotspot3D) of this program. The average runtime for 100 iterations is taken as the performance metric. All tests are executed on a single workstation with an Intel Xeon E5-2650 CPU. The experimental results show that the performance of OpenArray versions is comparable to the original versions.

Table 5 Four benchmark tests.

5.2 Validation tests of GOMO

The seamount problem proposed by Beckmann and Haidvogel is a widely used ideal test case for regional ocean models (Beckmann and Haidvogel, 1993). It is a stratified Taylor column problem which simulates the flow over an isolated seamount with a constant salinity and a reference vertical temperature stratification. An eastward horizontal current of 0.1 m s⁻¹ is added at model initialization. The southern and northern boundaries are closed. If the Rossby number is small, an obvious anticyclonic circulation is trapped by the mount in the deep water.

Using the seamount test case, we compare GOMO and sbPOM results. The configurations of both models are exactly the same. Figure 11 shows that GOMO and sbPOM both capture the anticyclonic circulation at 3500 m depth. The shaded plot shows the surface elevation, and the array plot shows the current at 3500 m. Panels 11a, b, and c show the results of GOMO, sbPOM, and the difference (GOMO-sbPOM), respectively. The differences in the surface elevation and deep currents between the two models are negligible (Fig. 11c).

Figure 11 Comparison of the surface elevation (shaded) and currents at 3500 m depth (vector) between GOMO and sbPOM on the 4th model day. (a) GOMO, (b) sbPOM, and (c) GOMO-sbPOM.

5.3 The weak and strong scalability of GOMO

The seamount test case is used to compare the performance of sbPOM and GOMO in a parallel environment. We use the X86 cluster at the National Supercomputing Center in Wuxi, China, which provides 5000 Intel Xeon E5-2650 v2 CPUs for our account at most. Figure 12a shows the result of a strong scaling evaluation, in which the model size is fixed at $\mathrm{2048}\times \mathrm{2048}\times \mathrm{50}$. The dashed line indicates the ideal speedup. For the largest parallelism with 4096 processes, GOMO and sbPOM achieve 91 % and 92 % parallel efficiency, respectively. Figure 12b shows the weak scalability of sbPOM and GOMO. In the weak scaling test, the model size for each process is fixed at $\mathrm{128}\times \mathrm{128}\times \mathrm{50}$, and the number of processes is gradually increased from 16 to 4096. Taking the performance of 16 processes as a baseline, we determine that the parallel efficiencies of GOMO and sbPOM using 4096 processes are 99.0 % and 99.2 %, respectively.

Figure 12 Performance comparison between sbPOM and GOMO on the X86 cluster. Panels (a) shows the strong scaling result; vertical axis denotes the speedup relative to 16 processes in a single node. Panel (b) shows the weak scaling result.

5.4 Testing on the Sunway platform

We also test the scalability of GOMO on the Sunway platform. Supposing that the baseline is the runtime of GOMO at 10 000 Sunway cores with a grid size of $\mathrm{4096}\times \mathrm{4096}\times \mathrm{50}$, the parallel efficiency of GOMO can still reach 85 % at 150 000 cores, as shown in Fig. 13. However, we notice that the scalability declines sharply when the number of cores exceeds 150 000. There are two reasons leading to this decline. First, the block size assigned to each core decreases as the number of cores increases, causing more communication during boundary region updating. Second, some processes cannot be accelerated even though more computing resources are available; for example, the time spent on creating the computation graph, generating the fusion kernels, and compiling the JIT cannot be reduced. Even though the fusion-kernel codes are generated and compiled only once at the beginning of a job, it consumes about 2 min. In a sense, OpenArray performs better when processing large-scale data, and GOMO is more suitable for high-resolution scenarios. In the future, we will further optimize the communication and graph-creating modules to improve the efficiency for large-scale cores.

Figure 13 Parallel efficiency of GOMO on the Sunway platform.

6 Discussion

As we mentioned in Sect. 1, the advantages of OpenArray are its ease of use, high efficiency, and portability. Using OpenArray, modelers without any parallel computing skill or experience can write simple operator expressions in Fortran to implement complex ocean models. The ocean models can be run on any CPU or Sunway platforms which have deployed the OpenArray library. We call this effect “write once, run everywhere”. Other similar libraries (e.g., ATMOL, ICON DSL, STELLA, and COARRAY) require users to manually control the boundary communication and task scheduling to some extent. In contrast, OpenArray implements completely implicit parallelism with user-friendly interfaces and programming languages.

However, there are still several problems to be solved in the development of OpenArray. The first issue is computational efficiency. Once a variable is in one of the processor registers or in the highest speed cache, it should be used as much as possible before being replaced. In fact, we should never move variables more than once each time step. The memory consumption brought by overloading techniques is usually high due to the unnecessary variable moving and unavoidable cache missing. The current efficiency and scalability of GOMO are close to sbPOM since we have adopted a series of optimization methods, such as memory pool, graph computing, JIT compilation, and vectorization, to alleviate the requirement of memory bandwidth. However, we have to admit that we cannot fully solve the memory bandwidth limit problem at present. We think that time skewing is a cache-oblivious algorithm for stencil computations (Frigo and Strumpen, 2005) since it can exploit temporal locality optimally throughout the entire memory hierarchy. In addition, the polyhedral model may be another potential approach, which uses an abstract mathematical representation based on integer polyhedral to analyze and optimize the memory access pattern of a program.

The second issue is that the current OpenArray version cannot support customized operators. When modelers try out another higher-order advection or any other numerical scheme, the 12 basic operators provided by OpenArray are not abundant. We considered using a template mechanism to support the customized operators. The rules of operations are defined in a template file, where the calculation form of each customized operator is described by a regular expression. If users want to add a customized operator, they only need to append a regular expression into the template file.

OpenArray and GOMO will continue to be developed, and the following three key improvements are planned for the following years.

First, we are developing the GPU version of OpenArray. During the development, the principle is to keep hot data in GPU memory or directly swapping between GPUs and avoid returning data to the main CPU memory. NVLink provides high bandwidth and outstanding scalability for GPU-to-CPU or GPU-to-GPU communication and addresses the interconnect issue for multi-GPU and multi-GPU/CPU systems.

Second, the data input/output is becoming a bottleneck of earth system models as the resolution increases rapidly. At present we encapsulate the PnetCDF library to provide simple I/O interfaces, such as load operation and store operation. A climate fast input/output (CFIO) library (Huang et al., 2014) will be implemented into OpenArray in the next few years. The performance of CFIO is approximately 220 % faster than PnetCDF because of the overlapping of I/O and computing. CFIO will be merged into the future version of OpenArray and the performance is expected to be further improved.

Finally, as for most of the ocean models, GOMO also faces the load imbalance issue. We are adding the more effective load balance schemes, including space-filling curve (Dennis, 2007) and curvilinear orthogonal grids, into OpenArray in order to reduce the computational cost on land points.

OpenArray is a product of the collaboration between oceanographers and computer scientists. It plays an important role in simplify the porting work on the Sunway TaihuLight supercomputer. We believe that OpenArray and GOMO will continue to be maintained and upgraded. We aim to promote it to the model community as a development tool for future numerical models.

7 Conclusions

In this paper, we design a simple computing library (OpenArray) to decouple ocean modeling and parallel computing. OpenArray provides 12 basic operators that are abstracted from PDEs and extended to ocean model governing equations. These operators feature user-friendly interfaces and an implicit parallelization ability. Furthermore, some state-of-the-art optimization mechanisms, including computation graphing, kernel fusion, dynamic source code generation, and JIT compiling, are applied to boost the performance. The experimental results prove that the performance of a program using OpenArray is comparable to that of well-designed programs using Fortran. Based on OpenArray, we implement a numerical ocean model (GOMO) with high productivity, enhanced readability, and excellent scalable performance. Moreover, GOMO shows high scalability on both CPU systems and the Sunway platform. Although more realistic tests are needed, OpenArray may signal the beginning of a new frontier in future ocean modeling through ingesting basic operators and cutting-edge computing techniques.

Appendix A: Continuous governing equations

The equations governing the baroclinic (internal) mode in GOMO are the 3-dimensional hydrostatic primitive equations.

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial UD}{\partial x}}+{\displaystyle \frac{\partial VD}{\partial y}}+{\displaystyle \frac{\partial W}{\partial \mathit{\sigma }}}=\mathrm{0},\text{(A2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial UD}{\partial t}}+{\displaystyle \frac{\partial U^{\mathrm{2}}D}{\partial x}}+{\displaystyle \frac{\partial UVD}{\partial y}}+{\displaystyle \frac{\partial UW}{\partial \mathit{\sigma }}}\\ & -fVD+gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}={\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left({\displaystyle \frac{K_M}{D}}{\displaystyle \frac{\partial U}{\partial \mathit{\sigma }}}\right)\\ & +{\displaystyle \frac{g{D}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial x}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\mathit{\rho }\mathrm{d}{\mathit{\sigma }}^{\prime }-{\displaystyle \frac{gD}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial D}{\partial x}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}^{\prime }{\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\sigma }}^{\prime }}}\mathrm{d}{\mathit{\sigma }}^{\prime }+F_u,\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial VD}{\partial t}}+{\displaystyle \frac{\partial UVD}{\partial x}}+{\displaystyle \frac{\partial V^{\mathrm{2}}D}{\partial y}}+{\displaystyle \frac{\partial VW}{\partial \mathit{\sigma }}}+fUD\\ & +gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial y}}={\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left({\displaystyle \frac{K_M}{D}}{\displaystyle \frac{\partial V}{\partial \mathit{\sigma }}}\right)+{\displaystyle \frac{g{D}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial y}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\mathit{\rho }\mathrm{d}{\mathit{\sigma }}^{\prime }\\ & -{\displaystyle \frac{gD}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial D}{\partial y}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}^{\prime }{\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\sigma }}^{\prime }}}\mathrm{d}{\mathit{\sigma }}^{\prime }+F_v,\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial TD}{\partial t}}+{\displaystyle \frac{\partial TUD}{\partial x}}+{\displaystyle \frac{\partial TVD}{\partial y}}+{\displaystyle \frac{\partial TW}{\partial \mathit{\sigma }}}={\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left(K_{\mathrm{H}}{\displaystyle \frac{\partial T}{\partial \mathit{\sigma }}}\right)\\ & +F_T+{\displaystyle \frac{\partial R}{\partial \mathit{\sigma }}},\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial SD}{\partial t}}+{\displaystyle \frac{\partial SUD}{\partial x}}+{\displaystyle \frac{\partial SVD}{\partial y}}+{\displaystyle \frac{\partial SW}{\partial \mathit{\sigma }}}=\\ & {\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left(K_{\mathrm{H}}{\displaystyle \frac{\partial S}{\partial \mathit{\sigma }}}\right)+F_S,\end{array}\text{(A6)}& {\displaystyle }\mathit{\rho }=\mathit{\rho }\left(TSp\right),\text{(A7)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial q^{\mathrm{2}}D}{\partial t}}+{\displaystyle \frac{\partial U{q}^{\mathrm{2}}D}{\partial x}}+{\displaystyle \frac{\partial V{q}^{\mathrm{2}}D}{\partial y}}+{\displaystyle \frac{\partial W{q}^{\mathrm{2}}}{\partial \mathit{\sigma }}}=\\ & {\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left({\displaystyle \frac{K_q}{D}}{\displaystyle \frac{\partial q^{\mathrm{2}}}{\partial \mathit{\sigma }}}\right)+{\displaystyle \frac{\mathrm{2}{K}_M}{D}}\left[{\left({\displaystyle \frac{\partial U}{\partial \mathit{\sigma }}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial V}{\partial \mathit{\sigma }}}\right)}^{\mathrm{2}}\right]\\ & +{\displaystyle \frac{\mathrm{2}g}{{\mathit{\rho }}_{\mathrm{0}}}}{K}_{\mathrm{H}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\sigma }}}-{\displaystyle \frac{\mathrm{2}D{q}^{\mathrm{3}}}{B_{\mathrm{1}}l}}+F_{q^{\mathrm{2}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(A9)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial q^{\mathrm{2}}lD}{\partial t}}+{\displaystyle \frac{\partial U{q}^{\mathrm{2}}lD}{\partial x}}+{\displaystyle \frac{\partial V{q}^{\mathrm{2}}lD}{\partial y}}+{\displaystyle \frac{\partial W{q}^{\mathrm{2}}l}{\partial \mathit{\sigma }}}=\\ & {\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\left({\displaystyle \frac{K_q}{D}}{\displaystyle \frac{\partial q^{\mathrm{2}}l}{\partial \mathit{\sigma }}}\right)+E_{\mathrm{1}}l\left\{{\displaystyle \frac{K_M}{D}}\left[{\left({\displaystyle \frac{\partial U}{\partial \mathit{\sigma }}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial V}{\partial \mathit{\sigma }}}\right)}^{\mathrm{2}}\right]\right.\\ & \left.+{\displaystyle \frac{g{E}_{\mathrm{3}}}{{\mathit{\rho }}_{\mathrm{0}}}}{K}_{\mathrm{H}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\sigma }}}\right\}\stackrel{\mathrm{̃}}{W}-{\displaystyle \frac{D{q}^{\mathrm{3}}}{B_{\mathrm{1}}}}+F_{q^{\mathrm{2}}l},\end{array}\end{array}$$

where F_u, F_v, $F_{q^{\mathrm{2}}}$, and $F_{q^{\mathrm{2}}l}$ are horizontal kinematic viscosity terms of u, v, q², and q²l, respectively. F_T and F_S are horizontal diffusion terms of T and S, respectively. $\stackrel{\mathrm{̃}}{W}$ is the wall proximity function.

$$\begin{array}{}\text{(A10)}& {\displaystyle }{F}_u={\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{2}{A}_{M}D{\displaystyle \frac{\partial U}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left[A_{M}D({\displaystyle \frac{\partial U}{\partial y}}+{\displaystyle \frac{\partial V}{\partial x}})\right],\text{(A11)}& {\displaystyle }{F}_v={\displaystyle \frac{\partial }{\partial y}}\left(\mathrm{2}{A}_{M}D{\displaystyle \frac{\partial V}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial x}}\left[A_{M}D({\displaystyle \frac{\partial U}{\partial y}}+{\displaystyle \frac{\partial V}{\partial x}})\right],\text{(A12)}& {\displaystyle }{F}_T={\displaystyle \frac{\partial }{\partial x}}\left(A_{\mathrm{H}}H{\displaystyle \frac{\partial T}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(A_{\mathrm{H}}H{\displaystyle \frac{\partial T}{\partial y}}\right),\text{(A13)}& {\displaystyle }{F}_S={\displaystyle \frac{\partial }{\partial x}}\left(A_{\mathrm{H}}H{\displaystyle \frac{\partial S}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(A_{\mathrm{H}}H{\displaystyle \frac{\partial S}{\partial y}}\right),\text{(A14)}& {\displaystyle }{F}_{q^{\mathrm{2}}}={\displaystyle \frac{\partial }{\partial x}}\left(A_{M}H{\displaystyle \frac{\partial q^{\mathrm{2}}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(A_{M}H{\displaystyle \frac{\partial q^{\mathrm{2}}}{\partial y}}\right),\text{(A15)}& {\displaystyle }{F}_{q^{\mathrm{2}}l}={\displaystyle \frac{\partial }{\partial x}}\left(A_{M}H{\displaystyle \frac{\partial q^{\mathrm{2}}l}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(A_{M}H{\displaystyle \frac{\partial q^{\mathrm{2}}l}{\partial y}}\right),\text{(A16)}& {\displaystyle }\stackrel{\mathrm{̃}}{W}=\mathrm{1}+{\displaystyle \frac{E_{\mathrm{2}}l}{\mathit{\kappa }}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\eta }-z}}+{\displaystyle \frac{\mathrm{1}}{H-z}}\right).\end{array}$$

The equations governing the barotropic (external) mode in GOMO are obtained by vertically integrating the baroclinic equations.

$$\begin{array}{}\text{(A17)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial U_{A}D}{\partial x}}+{\displaystyle \frac{\partial V_{A}D}{\partial y}}=\mathrm{0},\text{(A18)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial U_{A}D}{\partial t}}+{\displaystyle \frac{\partial (U_A{)}^{\mathrm{2}}D}{\partial x}}+{\displaystyle \frac{\partial U_A{V}_{A}D}{\partial y}}-f{V}_{A}D\\ & +gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}={\stackrel{\mathrm{̃}}{F}}_{u_a}-wu\left(\mathrm{0}\right)+wu(-\mathrm{1})\\ & -{\displaystyle \frac{gD}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\left[D{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}-{\displaystyle \frac{\partial D}{\partial x}}{\mathit{\sigma }}^{\prime }{\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\sigma }}}\right]\mathrm{d}{\mathit{\sigma }}^{\prime }\mathrm{d}\mathit{\sigma }+G_{u_a},\end{array}\text{(A19)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial V_{A}D}{\partial t}}+{\displaystyle \frac{\partial U_A{V}_{A}D}{\partial y}}+{\displaystyle \frac{\partial (V_A{)}^{\mathrm{2}}D}{\partial y}}+f{U}_{A}D\\ & +gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial y}}={\stackrel{\mathrm{̃}}{F}}_{v_a}-wv\left(\mathrm{0}\right)+wv(-\mathrm{1})\\ & -{\displaystyle \frac{gD}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\left[D{\displaystyle \frac{\partial \mathit{\rho }}{\partial y}}-{\displaystyle \frac{\partial D}{\partial y}}{\mathit{\sigma }}^{\prime }{\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\sigma }}}\right]\mathrm{d}{\mathit{\sigma }}^{\prime }\mathrm{d}\mathit{\sigma }+G_{v_a},\end{array}\end{array}$$

where ${\stackrel{\mathrm{̃}}{F}}_{u_a}$ and ${\stackrel{\mathrm{̃}}{F}}_{v_a}$ are the horizontal kinematic viscosity terms of U_A and V_A, respectively. $G_{u_a}$ and $G_{v_a}$ are the dispersion terms of U_A and V_A, respectively. The subscript “A” denotes vertical integration.

$$\begin{array}{}\text{(A20)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{̃}}{F}}_{u_a}={\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}H\left(A{A}_M\right){\displaystyle \frac{\partial U_A}{\partial x}}\right]\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[H\left(A{A}_M\right)\left({\displaystyle \frac{\partial U_A}{\partial y}}+{\displaystyle \frac{\partial V_A}{\partial x}}\right)\right],\end{array}\text{(A21)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{̃}}{F}}_{v_a}={\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}H\left(A{A}_M\right){\displaystyle \frac{\partial V_A}{\partial y}}\right]\\ & +{\displaystyle \frac{\partial }{\partial x}}\left[H\left(A{A}_M\right)\left({\displaystyle \frac{\partial U_A}{\partial y}}+{\displaystyle \frac{\partial V_A}{\partial x}}\right)\right],\end{array}\text{(A22)}& {\displaystyle }\begin{array}{rl}& G_{u_a}={\displaystyle \frac{{\partial }^{\mathrm{2}}(U_A{)}^{\mathrm{2}}D}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{U}_A{V}_{A}D}{\partial x\partial y}}-{\stackrel{\mathrm{̃}}{F}}_{u_a}\\ & -{\displaystyle \frac{{\partial }^{\mathrm{2}}(U^{\mathrm{2}}{)}_{A}D}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}(UV{)}_{A}D}{\partial y^{\mathrm{2}}}}+(F_u{)}_A,\end{array}\text{(A23)}& {\displaystyle }\begin{array}{rl}& G_{v_a}={\displaystyle \frac{{\partial }^{\mathrm{2}}{U}_A{V}_{A}D}{\partial x\partial y}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}(V_A{)}^{\mathrm{2}}D}{\partial y^{\mathrm{2}}}}-{\stackrel{\mathrm{̃}}{F}}_{v_a}\\ & -{\displaystyle \frac{{\partial }^{\mathrm{2}}(UV{)}_{A}D}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}(V^{\mathrm{2}}{)}_{A}D}{\partial y^{\mathrm{2}}}}+(F_v{)}_A,\end{array}\text{(A24)}& {\displaystyle }{U}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}Ud\mathit{\sigma },\text{(A25)}& {\displaystyle }{V}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}Vd\mathit{\sigma },\text{(A26)}& {\displaystyle }(U^{\mathrm{2}}{)}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}{U}^{\mathrm{2}}\mathrm{d}\mathit{\sigma },\text{(A27)}& {\displaystyle }(UV{)}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}UVd\mathit{\sigma },\text{(A28)}& {\displaystyle }(V^{\mathrm{2}}{)}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}{V}^{\mathrm{2}}\mathrm{d}\mathit{\sigma },\text{(A29)}& {\displaystyle }(F_u{)}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}{F}_u\mathrm{d}\mathit{\sigma },\text{(A30)}& {\displaystyle }(F_v{)}_A=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}{F}_v\mathrm{d}\mathit{\sigma },\text{(A31)}& {\displaystyle }A{A}_M=\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\left(A_M\right)\mathrm{d}\mathit{\sigma }.\end{array}$$

Appendix B: Discrete governing equations

The discrete governing equations of baroclinic (internal) mode expressed by operators are shown as below:

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\displaystyle \frac{{\mathit{\eta }}^{t+\mathrm{1}}-{\mathit{\eta }}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{x}U\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{y}V\right)+{\mathit{\delta }}_{\mathrm{f}}^z\left(W\right)=\mathrm{0},\text{(B2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{x}U{)}^{t+\mathrm{1}}-({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{x}U{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{b}}^x\left[({\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{x}U)}}_{\mathrm{f}}^x{\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^x\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{y}V\right)}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{U}}_{\mathrm{b}}^y\right]+{\mathit{\delta }}_{\mathrm{f}}^z\left({\stackrel{\mathrm{‾}}{W}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{U}}_{\mathrm{b}}^z\right)-{\stackrel{\mathrm{‾}}{\left(\stackrel{\mathrm{̃}}{f}{\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^{y}D\right)}}_{\mathrm{b}}^x\\ & -{\stackrel{\mathrm{‾}}{\left(f{\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^{y}D\right)}}_{\mathrm{b}}^x+g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x\left(\mathit{\eta }\right)={\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{K}}_{M_{\mathrm{b}}}^x}{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{)}^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z\left(U^{t+\mathrm{1}}\right)\right]\\ & +{\displaystyle \frac{g({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{)}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\left[{\mathit{\delta }}_{\mathrm{b}}^x\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^z\right)-{\displaystyle \frac{\mathit{\sigma }{\mathit{\delta }}_{\mathrm{b}}^x\left(D\right)}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x}}{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^x\right)\right]\mathrm{d}{\mathit{\sigma }}^{\prime }+F_u,\end{array}\text{(B3)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{y}V{)}^{t+\mathrm{1}}-({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{y}V{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left[({\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{x}U)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{V}}_{\mathrm{b}}^x\right]\\ & +{\mathit{\delta }}_{\mathrm{b}}^y\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^{y}V\right)}}_{\mathrm{f}}^y{\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^y\right]+{\mathit{\delta }}_{\mathrm{f}}^z\left({\stackrel{\mathrm{‾}}{W}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{V}}_{\mathrm{b}}^z\right)+{\stackrel{\mathrm{‾}}{\left(\stackrel{\mathrm{̃}}{f}{\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^{x}D\right)}}_{\mathrm{b}}^y\\ & +{\stackrel{\mathrm{‾}}{\left(f{\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^{x}D\right)}}_{\mathrm{b}}^y+g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y\left(\mathit{\eta }\right)={\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{K}}_{M_{\mathrm{b}}}^y}{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{)}^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z\left(V^{t+\mathrm{1}}\right)\right]\\ & +{\displaystyle \frac{g({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{)}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}\left[{\mathit{\delta }}_{\mathrm{b}}^y\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^z\right)-{\displaystyle \frac{\mathit{\sigma }{\mathit{\delta }}_{\mathrm{b}}^y\left(D\right)}{{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y}}{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^y\right)\right]\mathrm{d}{\mathit{\sigma }}^{\prime }+F_v,\end{array}\text{(B4)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{(TD{)}^{t+\mathrm{1}}-(TD{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{T}}_{\mathrm{b}}^{x}U{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{T}}_{\mathrm{b}}^{y}V{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\right)\\ & +{\mathit{\delta }}_{\mathrm{f}}^z\left({\stackrel{\mathrm{‾}}{T}}_{\mathrm{b}}^{z}W\right)={\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{K_{\mathrm{H}}}{D^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z\left(T^{t+\mathrm{1}}\right)\right]+F_T+{\mathit{\delta }}_{\mathrm{f}}^{z}R,\end{array}\text{(B5)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{(SD{)}^{t+\mathrm{1}}-(SD{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{S}}_{\mathrm{b}}^{x}U{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{S}}_{\mathrm{b}}^{y}V{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\right)+{\mathit{\delta }}_{\mathrm{f}}^z\left({\stackrel{\mathrm{‾}}{S}}_{\mathrm{b}}^{z}W\right)=\\ & {\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{K_{\mathrm{H}}}{D^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z\left(S^{t+\mathrm{1}}\right)\right]+F_S,\end{array}\text{(B6)}& {\displaystyle }\mathit{\rho }=\mathit{\rho }\left(TSp\right),\text{(B7)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{(q^{\mathrm{2}}D{)}^{t+\mathrm{1}}-(q^{\mathrm{2}}D{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{U}}_{\mathrm{b}}^z{\stackrel{\mathrm{‾}}{(q^{\mathrm{2}}{)}_{\mathrm{b}}}}^x{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}}}^x\right)\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{V_{\mathrm{b}}}}^z{\stackrel{\mathrm{‾}}{(q^{\mathrm{2}}{)}_{\mathrm{b}}}}^y{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}}}^y\right)+{\mathit{\delta }}_{\mathrm{f}}^z{\stackrel{\mathrm{‾}}{\left(W{q}^{\mathrm{2}}\right)}}_{\mathrm{b}}^z=\\ & {\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{(K_q{)}_{\mathrm{f}}}}^z}{D^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z(q^{\mathrm{2}}{)}^{t+\mathrm{1}}\right]+{\displaystyle \frac{\mathrm{2}{K}_M}{D}}\left\{{\left[{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^x\right)\right]}^{\mathrm{2}}\right.\\ & \left.+{\left[{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^y\right)\right]}^{\mathrm{2}}\right\}+{\displaystyle \frac{\mathrm{2}g}{{\mathit{\rho }}_{\mathrm{0}}}}{K}_{\mathrm{H}}{\mathit{\delta }}_{\mathrm{b}}^z\left(\mathit{\rho }\right)-{\displaystyle \frac{\mathrm{2}D{q}^{\mathrm{3}}}{B_{\mathrm{1}}l}}+F_{q^{\mathrm{2}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(B8)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{(q^{\mathrm{2}}lD{)}^{t+\mathrm{1}}-(q^{\mathrm{2}}lD{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}ti}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{U}}_{\mathrm{b}}^z{\stackrel{\mathrm{‾}}{\left(q^{\mathrm{2}}l\right)}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{V}}_{\mathrm{b}}^z{\stackrel{\mathrm{‾}}{q^{\mathrm{2}}l}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y\right)+{\mathit{\delta }}_{\mathrm{f}}^z{\stackrel{\mathrm{‾}}{\left(W{q}^{\mathrm{2}}l\right)}}_{\mathrm{b}}^z=\\ & {\mathit{\delta }}_{\mathrm{b}}^z\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{(K_q{)}_{\mathrm{f}}}}^z}{D^{t+\mathrm{1}}}}{\mathit{\delta }}_{\mathrm{f}}^z(q^{\mathrm{2}}l{)}^{t+\mathrm{1}}\right]+l{E}_{\mathrm{1}}{\displaystyle \frac{K_M}{D}}\left\{{\left[{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{U}}_{\mathrm{f}}^x\right)\right]}^{\mathrm{2}}\right.\\ & \left.+{\left[{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{V}}_{\mathrm{f}}^y\right)\right]}^{\mathrm{2}}\right\}\stackrel{\mathrm{̃}}{W}+{\displaystyle \frac{l{E}_{\mathrm{1}}{E}_{\mathrm{3}}g}{{\mathit{\rho }}_{\mathrm{0}}}}{K}_{\mathrm{H}}{\mathit{\delta }}_{\mathrm{b}}^z\left(\mathit{\rho }\right)\stackrel{\mathrm{̃}}{W}-{\displaystyle \frac{D{q}^{\mathrm{3}}}{B_{\mathrm{1}}}}\\ & +F_{q^{\mathrm{2}}l},\end{array}\end{array}$$

where F_u, F_v, $F_{q^{\mathrm{2}}}$, and $F_{q^{\mathrm{2}}l}$ are horizontal kinematic viscosity terms of u, v, q², and q²l, respectively. F_T and F_S are horizontal diffusion terms of T and S, respectively.

$$\begin{array}{}\text{(B9)}& {\displaystyle }\begin{array}{rl}& F_u={\mathit{\delta }}_{\mathrm{b}}^x\left[\mathrm{2}{A}_{M}D{\mathit{\delta }}_{\mathrm{f}}^x\left(U^{t-\mathrm{1}}\right)\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left\{{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{\left(A_M\right)}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y\left[{\mathit{\delta }}_{\mathrm{b}}^x\left(V^{t-\mathrm{1}}\right)+{\mathit{\delta }}_{\mathrm{b}}^y\left(U^{t-\mathrm{1}}\right)\right]\right\},\end{array}\text{(B10)}& {\displaystyle }\begin{array}{rl}& F_v={\mathit{\delta }}_{\mathrm{b}}^y\left[\mathrm{2}{A}_{M}D{\mathit{\delta }}_{\mathrm{f}}^y\left(V^{t-\mathrm{1}}\right)\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^x\left\{{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{\left(A_M\right)}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y\left[{\mathit{\delta }}_{\mathrm{b}}^x\left(V^{t-\mathrm{1}}\right)+{\mathit{\delta }}_{\mathrm{b}}^y\left(U^{t-\mathrm{1}}\right)\right]\right\},\end{array}\text{(B11)}& {\displaystyle }\begin{array}{rl}& F_T={\mathit{\delta }}_{\mathrm{f}}^x\left[{\stackrel{\mathrm{‾}}{\left(A_{\mathrm{H}}\right)}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x\left(T^{t-\mathrm{1}}\right)\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{\left(A_{\mathrm{H}}\right)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y\left(T^{t-\mathrm{1}}\right)\right],\end{array}\text{(B12)}& {\displaystyle }\begin{array}{rl}& F_S={\mathit{\delta }}_{\mathrm{f}}^x\left[\left({\stackrel{\mathrm{‾}}{\left(A_{\mathrm{H}}\right)}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x\right(S^{t-\mathrm{1}})\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{\left(A_{\mathrm{H}}\right)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y\left(S^{t-\mathrm{1}}\right)\right],\end{array}\text{(B13)}& {\displaystyle }\begin{array}{rl}& F_{q^{\mathrm{2}}}={\mathit{\delta }}_{\mathrm{f}}^x\left[{\stackrel{\mathrm{‾}}{({\stackrel{\mathrm{‾}}{\left(A_M\right)}}_{\mathrm{b}}^x{)}_{\mathrm{b}}}}^z{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x(q^{\mathrm{2}}{)}^{t-\mathrm{1}}\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{({\stackrel{\mathrm{‾}}{\left(A_M\right)}}_{\mathrm{b}}^y{)}_{\mathrm{b}}}}^z{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y(q^{\mathrm{2}}{)}^{t-\mathrm{1}}\right],\end{array}\text{(B14)}& {\displaystyle }\begin{array}{rl}& F_{q^{\mathrm{2}}l}={\mathit{\delta }}_{\mathrm{f}}^x\left[{\stackrel{\mathrm{‾}}{({\stackrel{\mathrm{‾}}{(A_M{)}_{\mathrm{b}}}}^x{)}_{\mathrm{b}}}}^z{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x(q^{\mathrm{2}}l{)}^{t-\mathrm{1}}\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{({\stackrel{\mathrm{‾}}{(A_M{)}_{\mathrm{b}}}}^y{)}_{\mathrm{b}}}}^z{\stackrel{\mathrm{‾}}{H}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y(q^{\mathrm{2}}l{)}^{t-\mathrm{1}}\right].\end{array}\end{array}$$

The discrete governing equations of barotropic (external) mode expressed by operators are shown as below:

$$\begin{array}{}\text{(B15)}& {\displaystyle }{\displaystyle \frac{{\mathit{\eta }}^{t+\mathrm{1}}-{\mathit{\eta }}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}te}}+{\mathit{\delta }}_{\mathrm{f}}^x\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{U}_A\right)+{\mathit{\delta }}_{\mathrm{f}}^y\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{V}_A\right)=\mathrm{0},\end{array}$$

$$\begin{array}{}\text{(B16)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{U}_A{)}^{t+\mathrm{1}}-({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{U}_A{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}te}}+{\mathit{\delta }}_{\mathrm{b}}^x\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{U}_A\right)}}_{\mathrm{f}}^x{\stackrel{\mathrm{‾}}{\left(U_A\right)}}_{\mathrm{f}}^x\right]\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{V}_A\right)}}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{\left(U_A\right)}}_{\mathrm{b}}^y\right]-{\stackrel{\mathrm{‾}}{\left[{\stackrel{\mathrm{̃}}{f}}_A{\stackrel{\mathrm{‾}}{\left(V_A\right)}}_{\mathrm{f}}^{y}D\right]}}_{\mathrm{b}}^x\\ & -{\stackrel{\mathrm{‾}}{\left[f{\stackrel{\mathrm{‾}}{\left(V_A\right)}}_{\mathrm{f}}^{y}D\right]}}_{\mathrm{b}}^x+g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^x\left(\mathit{\eta }\right)=\\ & {\mathit{\delta }}_{\mathrm{b}}^x\left\{\mathrm{2}\left(A{A}_M\right)D{\mathit{\delta }}_{\mathrm{f}}^x\left[(U_A{)}^{t-\mathrm{1}}\right]\right\}\\ & +{\mathit{\delta }}_{\mathrm{f}}^y\left\{{\stackrel{\mathrm{‾}}{\left[{\stackrel{\mathrm{‾}}{\left(A{A}_M\right)}}_{\mathrm{b}}^x\right]}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y{\left[{\mathit{\delta }}_{\mathrm{b}}^x\left(V_A\right)+{\mathit{\delta }}_{\mathrm{b}}^y\left(U_A\right)\right]}^{t-\mathrm{1}}\right\}\\ & +{\mathit{\varphi }}_x,\end{array}\text{(B17)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{V}_A{)}^{t+\mathrm{1}}-({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{V}_A{)}^{t-\mathrm{1}}}{\mathrm{2}\mathrm{d}te}}+{\mathit{\delta }}_{\mathrm{f}}^x\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{U}_A\right)}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left(V_A\right)}}_{\mathrm{b}}^x\right]\\ & +{\mathit{\delta }}_{\mathrm{b}}^y\left[{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{V}_A\right)}}_{\mathrm{f}}^y{\stackrel{\mathrm{‾}}{\left(V_A\right)}}_{\mathrm{f}}^y\right]+{\stackrel{\mathrm{‾}}{\left[{\stackrel{\mathrm{̃}}{f}}_A{\stackrel{\mathrm{‾}}{\left(U_A\right)}}_{\mathrm{f}}^{x}D\right]}}_{\mathrm{b}}^y\\ & +{\stackrel{\mathrm{‾}}{\left[f{\stackrel{\mathrm{‾}}{\left(U_A\right)}}_{\mathrm{f}}^{x}D\right]}}_{\mathrm{b}}^y+g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^y\left(\mathit{\eta }\right)=\\ & {\mathit{\delta }}_{\mathrm{b}}^y\left\{\mathrm{2}\left(A{A}_M\right)D{\mathit{\delta }}_{\mathrm{f}}^y\left[(V_A{)}^{t-\mathrm{1}}\right]\right\}+\\ & {\mathit{\delta }}_{\mathrm{f}}^x\left\{{\stackrel{\mathrm{‾}}{\left[{\stackrel{\mathrm{‾}}{\left(A{A}_M\right)}}_{\mathrm{b}}^x\right]}}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x\right)}}_{\mathrm{b}}^y{\left[{\mathit{\delta }}_{\mathrm{b}}^x\left(V_A\right)+{\mathit{\delta }}_{\mathrm{b}}^y\left(U_A\right)\right]}^{t-\mathrm{1}}\right\}\\ & +{\mathit{\varphi }}_y,\end{array}\end{array}$$

where

$$\begin{array}{}\text{(B18)}& {\displaystyle }\begin{array}{rl}& {\mathit{\varphi }}_x=-WU\left(\mathrm{0}\right)+WU(-\mathrm{1})\\ & -{\displaystyle \frac{g({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{)}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\left\{\left[\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\mathit{\delta }}_{\mathrm{b}}^x{\stackrel{\mathrm{‾}}{\left(\mathit{\rho }\right)}}_{\mathrm{b}}^z\mathrm{d}{\mathit{\sigma }}^{\prime }\right]\mathrm{d}\mathit{\sigma }\right\}\\ & +{\displaystyle \frac{g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^x{\mathit{\delta }}_{\mathrm{b}}^{x}D}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\left\{\left[\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{b}}^z{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^x\right)\right]\mathrm{d}\mathit{\sigma }\right\}+G_x,\end{array}\text{(B19)}& {\displaystyle }\begin{array}{rl}& {\mathit{\varphi }}_y=-WV\left(\mathrm{0}\right)+WV(-\mathrm{1})\\ & -{\displaystyle \frac{g({\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{)}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\left\{\left[\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\mathit{\delta }}_{\mathrm{b}}^y{\stackrel{\mathrm{‾}}{\left(\mathit{\rho }\right)}}_{\mathrm{b}}^z\mathrm{d}{\mathit{\sigma }}^{\prime }\right]\mathrm{d}\mathit{\sigma }\right\}\\ & +{\displaystyle \frac{g{\stackrel{\mathrm{‾}}{D}}_{\mathrm{b}}^y{\mathit{\delta }}_{\mathrm{b}}^{y}D}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}\left\{\left[\underset{\mathit{\sigma }}{\overset{\mathrm{0}}{\int }}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{b}}^z{\mathit{\delta }}_{\mathrm{b}}^z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{b}}^y\right)\right]\mathrm{d}\mathit{\sigma }\right\}+G_y.\end{array}\end{array}$$

Appendix C: Descriptions of symbols

The description of each symbol in the governing equations is listed below.

Table C1 Descriptions of symbols.