Department of Mathematics
High-order compact finite difference schemes
for the spherical shallow water equations
Sarmad Ghader and Jan Nordstrm
LiTH-MAT-R–2013/09–SE
High-order compact finite difference schemes for the spherical
shallow water equations
Sarmad Ghader
a∗and Jan Nordstr¨
om
baInstitute of Geophysics, University of Tehran, Iran
bDivision of Computational Mathematics, Department of Mathematics, Link¨oping University, Sweden
SUMMARY
This work is devoted to the application of the super compact finite difference (SCFDM) and the combined compact finite difference (CCFDM) methods for spatial differencing of the spherical shallow water equations in terms of vorticity, divergence and height. Five high-order schemes including the fourth-order compact, the sixth-order and eighth-order SCFDM and the sixth-order and eighth-order CCFDM schemes are used for spatial differencing of the spherical shallow water equations. To advance the solution in time, a semi-implicit Runge-Kutta method is used. In addition, to control the non-linear instability and avoiding the polar problem a high-order spatial filter is proposed. An unstable barotropic mid-latitude zonal jet is employed as an initial condition. For the numerical solution of the elliptic equations in the problem, a direct hybrid method which consists of using a high-order compact scheme for spatial differencing in the latitude coordinate and a fast Fourier transform in longitude coordinate is utilized. The convergence rate for all methods is studied and verified. Qualitative and quantitative assessment of the results, such as measures of maximum vorticity gradient, power spectrum of total energy, relative change in potential enstrophy and potential palinstrophy, reveal that the sixth-order and eighth-order CCFDM and the sixth-order and eighth-order SCFDM methods lead to a remarkable improvement of the solution over the fourth-order compact method. It is also shown that the performance of the sixth-order and eighth-order CCFDM methods are superior to the sixth-order and eighth-order SCFDM methods. Copyright c⃝ 2013 John Wiley & Sons, Ltd. key words: Compact finite difference; High-order methods; Spherical shallow water equations; Numerical accuracy
1. Introduction
The shallow water equations are widely used in various oceanic and atmospheric problems. This model is applied to a fluid layer of constant density in which the horizontal scale of the flow is much greater than the layer depth (e.g., [1, 2]). Even though, the dynamics of two dimensional shallow water model is less general than the three dimensional general circulation models, it is often preferred because of its mathematical and computational simplicity. The
∗Correspondence to: Sarmad Ghader, Institute of Geophysics, University of Tehran, North Kargar Ave., Tehran,
spherical shallow water equations are usually used as the first step to examine new numerical algorithms to be used in more complex global climate and weather prediction models.
Compact finite-difference schemes are simple and powerful ways to reach the objectives of high accuracy and low computational cost (e.g., [3, 4, 5]). Compared with the traditional explicit finite-difference schemes of the same-order, compact schemes are significantly more accurate along with the benefit of using smaller stencil sizes. The compact finite difference methods can not be formulated on so called summation-by-parts (SBP) form [6, 7, 8, 9] which is a serious drawback for finite domain problems since the stability can not be proved. However, in this work we only consider periodic domains and this drawback is not important.
One of the main problems with the spherical harmonic method [10], which is the dominant global modelling method, is its rapid increase of computational cost with resolution. Therefore, examination of other high-order schemes, such as compact finite difference methods, as alternatives is important. Previous applications of compact schemes to idealized models of the atmosphere and oceans have shown that they are a promising alternative for the numerical simulation of geophysical fluid dynamics problems (e.g., [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]). Motivated by the remarkable performance of the sixth-order and eighth-order super compact finite difference methods compared with the fourth-order compact method in Ref. [18], this paper is devoted to the application of two families of high-order compact finite difference methods to the spatial differencing of the spherical shallow water equations. The super compact finite difference (SCFDM) and the combined compact finite difference (CCFDM) methods [26, 27, 21, 22, 23, 24, 25] are considered in this study.
The remainder of this paper is organized as follows. The spherical shallow water equations are presented in section 2. Details of the semi-implicit Runge-Kutta time advancing method and formulation of the high-order compact methods used for spatial differencing of the spherical shallow water equations are given in sections 3 and 4, respectively. The procedure used for the numerical solution of the elliptic equations is described in section 5. The spatial filter formulation is given in section 6. The order of convergence for the high-order compact methods used in this study is presented in section 7. Details of the test case used for the numerical solution, numerical results and discussions are given in section 8. Finally, section 9 presents the conclusions.
2. The shallow water equations
The momentum and mass continuity representation of the spherical shallow water equations in vector form can be written as (e.g., [28])
Dv
Dt + f ˆz× v + c
2∇˜h = 0 (1)
∂˜h
∂t +∇ · [(1 + ˜h)v] = 0. (2)
In (1) and (2), v = uˆi + vˆj is the horizontal velocity vector with u and v being the velocity components in longitudinal (λ) and latitudinal (ϕ) directions, respectively. ˆi and ˆj are the unit vectors in longitudinal and latitudinal directions, respectively. D()/Dt = ∂()/∂t + (v· ∇)() is the substantial derivative and f = 2Ω sin ϕ is the Coriolis parameter in which Ω is the rotation rate of the earth. ˜h = (h−H)/H is the dimensionless perturbation depth where h is the actual
depth where H is the domain area average depth. c = √gH is the gravity wave speed and g
denotes the acceleration due to gravity. The unit vector in vertical direction is denoted by ˆz.
Equations (1) and (2) can be rewritten in terms of relative vorticity, horizontal divergence and height as ∂ζ ∂t = Sζ (3) ∂δ ∂t +H˜h = Sδ (4) ∂˜h ∂t + δ = Sh. (5) In equations (3)-(5), Sζ =−∇ · [(ζ + f)v], Sδ= f (ζ− f˜h) − ∇2(|v| 2 2 )− ∇ · (−ζvˆi + ζuˆj) − u ˆβ, Sh=−∇ · (˜hv) and H = c2∇2− f2is the modified Helmholtz operator, ˆβ = (df /dϕ)/a = (2Ω/a) cos ϕ and a
is the radius of the earth. The relative vorticity and divergence are defined as
ζ = 1 a cos ϕ ( ∂v ∂λ − ∂u cos ϕ ∂ϕ ) , δ = 1 a cos ϕ ( ∂u ∂λ + ∂v cos ϕ ∂ϕ ) . (6)
In addition, to find the velocity components and to close the system of equations the Helmholtz decomposition v = ˆz×∇ψ+∇χ is employed where ψ and χ denote streamfunction and velocity
potential, respectively. To find the streamfunction and velocity potential, the following Poisson equations must be solved:
∇2ψ = ζ, ∇2χ = δ. (7)
3. The semi-implicit Runge-Kutta method
In the present work the semi-implicit Runge-Kutta (hereafter, SIRK) time advancing method developed in [29] for the flux convergence form of the shallow water equations in Cartesian coordinates, is extended to the vorticity, divergence and height representation of the spherical shallow water equations. For the vorticity-divergence form of the shallow water equations, the SIRK formulation is applied to equations (4) and (5) [30]. The Runge-Kutta formulation used in the present work is a third-order explicit three step method [31].
Equations (4) and (5) can be rewritten in the following vector form:
∂q
∂t = L(q) + N (q) (8)
where q = (δ, ˜h)T, L(q) = (−H˜h, −δ)T and N (q) = (S
δ, Sh)T. The superscript T denotes
transpose. The semi-implicit approach which employs the third-order Runge-Kutta method for time advancement of equation (8) leads to
qn+αk− qn αk∆t = L(q n+αk) + L(qn) 2 + N (q n+αk−1) (9)
where the superscript n denotes the time step. The three steps of the Runge-Kutta method are written in a compact form where subscript k (k = 1, 2, 3) denotes the three steps of the Runge-Kutta method and the values of αk are α0 = 0, α1= 1/3, α2= 1/2 and α3= 1. The following auxiliary variables are defined
¯ δn+αk= δ n+αk+ δn 2 , ¯ hn+αk= ˜ hn+αk+ ˜hn 2 . (10)
By substituting the auxiliary variables (10) into equation (9), the following equations for the auxiliary variables can be found
(H − 4 α2 k∆t2 )¯hn+αk = 2 αk∆t (δn− Shn+αk−1) + Sδn+αk−1− 4 α2 k∆t2 ˜ hn (11) ¯ δn+αk = 2 αk∆t (˜hn− ¯hn+αk) + Sn+αk−1 h . (12)
Equation (11) is a modified Helmholtz equation. The unknown variables δn+αk and ˜hn+αk are
found by substitution of the auxiliary variables into equation (10). The time discretization of the vorticity equation (3) by the third-order Runge-Kutta method leads to:
ζn+αk= ζn+ α
k∆tS n+αk−1
ζ (13)
where the definitions of n and αk are given above.
4. Spatial differencing
For spatial differencing of the first and second derivatives in equations (11), (12), (13) and the diagnostic equations (7), two families of compact finite difference schemes, i.e., the super compact finite difference method (SCFDM) and the combined compact finite difference method (CCFDM) are used.
In the SCFDM, approximations of the first (second) derivatives depends on all the nodal values and their odd (even) derivatives, but the CCFDM approximates the first and second derivatives simultaneously. See references [26, 27, 16] and [21, 22, 23, 24, 25] for more information on details of the derivation and alternative forms of the SCFDM and CCFDM, respectively.
4.1. The super compact method (SCFDM)
The SCFDM relations for approximation of the first derivative of an arbitrary function Φ in a uniform grid are
where Fj = Φ<1>j Φ<3>j .. . Φ<2Mj −1> , EFj = D0Φj 0 .. . 0 , Qf = 1 d 1 1! 1 3! 1 5! . . . 1 (2M− 1)! −d2 2D 2 1 2! 1 4! . . . 1 [2(M− 1)]! 0 −d22D2 1 2! . . . 1 [2(M− 2)]! . . . . . . . . . . . . . . . 0 0 . . . −d22D2 1 2! In (14), the difference operators are defined as
D0Φj= (Φj+1− Φj−1)/2d, D2Φj = (Φj+1− 2Φj+ Φj−1)/d2
in which d denotes the grid spacing and Φ<l>j = dl(∂lΦ
∂xl
)
j
. Qf is a M× M matrix, F and EFare M dimensional vectors. The expression Φ<2lj −1>/d2l−1 (l = 1, 2, . . . , M ) approximates
∂2l−1Φ
∂x2l−1 with an accuracy of order 2(M− l + 1) [26, 27].
In a similar way, the SCFDM relations for approximation of the second derivative on a uniform grid are
QsSj= ESj (15) in which Sj= Φ<2>j Φ<4>j .. . Φ<2M >j , E S j = 1 2D 2Φ j 0 .. . 0 , Qs= 1 d2 1 2! 1 4! 1 6! . . . 1 (2M )! −d2 2D 2 1 2! 1 4! . . . 1 [2(M− 1)]! 0 −d2 2D 2 1 2! . . . 1 [2(M− 2)]! . . . . . . . . . . . . . . . 0 0 . . . −d22D2 1 2! Here, Qs is a M× M matrix, S and ESare M dimensional vector. The expression Φ<2l>j /d2l
approximates ∂ 2lϕ
∂x2l with an accuracy of order 2(M− l + 1) [26, 27].
Using different values of M (M = 1, 2, 3,· · · ) in equations (14) and (15) leads to different orders of the SCFDM relations for approximation of the first and second derivatives. For example, M = 1 leads to the conventional second-order central finite difference formulation. The fourth-order Pade type compact finite difference are found by using M = 2. Higher values of M leads to the higher order SCFDM relations.
4.1.1. The SCFDM relations Using M = 2, 3, 4 in equation (14) leads to the fourth-order, sixth-order and eighth-order SCFDM relations for the approximation of the first derivative. The fourth-order relation reads
1 6Φ (1) j−1+ 2 3Φ (1) j + 1 6Φ (1) j+1− 1 2d(Φj+1− Φj−1) = 0, (16)
the sixth-order relations are { Φj+1− Φj−1− 2dΦ (1) j − d3 60Φ (3) j−1− 3d3 10Φ (3) j − d3 60Φ (3) j+1= 0 Φ(1)j+1− 2Φ(1)j + Φj(1)−1−d122Φ(3)j−1−5d62Φ(3)j −d122Φ(3)j+1= 0, (17) and the eighth-order relations are found as
Φj+1− Φj−1− 2dΦ (1) j − d3 3Φ (3) j − 2d5 5040Φ (5) j−1− 2d5 126Φ (5) j − 2d5 5040Φ (5) j+1= 0 Φ(1)j+1− 2Φ(1)j + Φ(1)j−1− d2Φ(3)j −360d4 Φ(5)j−1−7d904Φj(5)−360d4 Φ(5)j+1= 0 Φ(3)j+1− 2Φ(3)j + Φ(3)j−1−d2 12Φ (5) j−1− 5d2 6 Φ (5) j − d2 12Φ (5) j+1= 0. (18)
Similarly, using M = 2, 3, 4 in equation (15) leads to the fourth-order, sixth-order and eighth-order SCFDM relations for the approximation of the second derivative. The fourth-order relation is 1 12Φ (2) j−1+ 5 6Φ (2) j + 1 12Φ (2) j+1− 1 d2(Φj+1− 2Φj+ Φj−1) = 0, (19) the sixth-order relations are
{ Φj+1− 2Φj+ Φj−1− d2Φ (2) j − d4 360Φ (4) j−1− 7d4 90Φ (4) j − d4 360Φ (4) j+1= 0 Φ(2)j+1− 2Φ(2)j + Φ(2)j−1−d2 12Φ (4) j−1− 5d2 6 Φ (4) j − d2 12Φ (4) j+1= 0, (20)
and the eighth-order relations are found as Φj+1− 2Φj+1+ Φj−1− d2Φ (2) j − d4 12Φ (4) j − d6 20160Φ (6) j−1− 3d6 1120Φ (6) j − d6 20160Φ (6) j+1= 0 Φ(2)j+1− 2Φ(2)j + Φ(2)j−1− d2Φ(4) j − d4 360Φ (6) j−1− 7d4 90Φ (6) j − d4 360Φ (6) j+1= 0 Φ(4)j+1− 2Φ(4)j + Φ(4)j−1−d122Φj(6)−1−5d62Φj(6)−d122Φ(6)j+1= 0. (21)
In these equations, the superscript (k) denotes kth spatial derivative.
4.2. The combined compact method (CCFDM)
The CCFDM relations for the first and second derivatives on a uniform grid can be written as
QcTj= Gj (22) where Tj = Φ<1>j Φ<2>j .. . Φ<2Nj −1> Φ<2N >j , Gj = D0Φj D2Φj 0 .. . 0
Qc= a 1! 0 a 3! 0 . . . a (2N− 1)! 0 0 2b 2! 0 2b 4! . . . 0 2b (2N )! −D0 a 1! 0 a 3! . . . 0 a (2N − 1)! −D2 0 2b 2! 0 . . . 2b (2N )! 0 0 −D0 a 1! 0 . . . 0 a (2N − 3)! 0 −D2 0 2b 2! . . . 2b (2N− 2)! 0 0 0 −D0 a 1! . . . 0 a (2N − 5)! 0 0 −D2 0 . . . 2b (2N− 4)! 0 .. . ... ... ... ... 0 0 0 0 . . . . . . . . . and a = 1/d, b = 1/d2. For the CCFDM method, Q
c is a 2N× 2N matrix and T and G are 2N dimensional vectors. Equation (22) approximates both the first and second derivatives of Φ with an accuracy of order 2N .
Substitution of different values of N (N = 1, 2, 3,· · · ) in equation (22) leads to various orders of the CCFDM relations for approximation of the first and second derivatives simultaneously. Similar to the SCFDM, N = 1 and N = 2 lead to the second-order central finite difference and fourth-order Pade type compact finite difference relations. Higher values of N leads to the higher order CCFDM relations.
4.2.1. The CCFDM relations For the CCFDM method by using N = 3, 4 in equation (22), the sixth-order and eighth-order relations for the approximation of the first and second derivatives are obtained. The sixth-order relations are
{ 7 16(Φ (1) j+1+ Φ (1) j−1) + Φ (1) j − d 16(Φ (2) j+1− Φ (2) j−1)− 15 16d(Φj+1− Φj−1) = 0 9 8d(Φ (1) j+1− Φ (1) j−1)− 1 8(Φ (2) j+1+ Φ (2) j−1) + Φ (2) j − 3 d2(Φj+1− 2Φj+ Φj−1) = 0, (23)
and the eighth-order relations are found as 19 32(Φ (1) j+1+ Φ (1) j−1) + Φ (1) j − d 8(Φ (2) j+1− Φ (2) j−1) + d2 96(Φ (3) j+1+ Φ (3) j−1)− 35 32d(Φj+1− Φj−1) = 0 29 16d(Φ (1) j+1− Φ (1) j−1)− 5 16(Φ (2) j+1+ Φ (2) j−1) + Φ (2) j + d 48(Φ (3) j+1− Φ (3) j−1)− 4 d2(Φj+1− 2Φj+ Φj−1) = 0 −105 16d2(Φ (1) j+1+ Φ (1) j−1) + 15 8d(Φ (2) j+1− Φ (2) j−1)− 3 16(Φ (3) j+1+ Φ (3) j−1) + Φ (3) j + 105 16d3(Φj+1− Φj−1) = 0. (24)
4.3. The numerical grid
The shallow water equations are solved on the doubly periodic longitude-latitude grid system (0≤ λ ≤ 2π, −π/2 ≤ ϕ ≤ π/2) of reference [32]. The grid points in this system are arranged such that no grid points are located at the poles. This grid system has the advantage of avoiding singularity at the poles and enabling the periodic boundary conditions to be imposed
in both longitude and latitude directions. It should be noted that the sign of a quantity must be adjusted properly when using the periodic boundary condition in the latitude direction, depending on whether the quantity is a scalar or a vector. The grid points are defined as (λi, ϕj) = (i∆λ,−π/2 + (j − 1/2)∆ϕ) where ∆λ = 2π/Nλ and ∆ϕ = π/Nϕ. Nλ and Nϕ are
the number of grid points in latitude and longitude directions, respectively.
5. The elliptic equations
The SIRK formulation used to advance the shallow water equations in time leads to a set of modified Helmholtz equations (11) at each successive step of the Runge-Kutta method. In addition, to find the velocity components, the two Poisson equations given in (7) must be solved at each step.
The elliptic equations (11) and (7) are solved using a direct method. The procedure is similar to one developed in [33, 34]. However, the following changes are applied to make the procedure easier. The spatial discretization of the derivatives in longitude are performed by discrete Fourier series and the FFT (Fast Fourier Transform) technique. In the latitudinal direction we discretize using different orders of the SCFDM and CCFDM schemes.
5.1. The Helmholtz equation
The modified Helmholtz equation (11) can be written as
∇2¯h− β¯h = R (25)
where the Laplacian operator is
∇2= 1 a2cos2ϕ ∂2 ∂λ2+ 1 a2cos ϕ ( ∂ ∂ϕcos ϕ ∂ ∂ϕ ) .
In (25), β is a function of ϕ and R denotes the known right hand side of the elliptic equation. The longitudinal parts of ¯h and R in (25) are approximated by truncated Fourier series. By
substitution of the Fourier components ¯h = ˆhmeimλ and R = ˆRmeimλ with wave number m
into equation (25) we have ( −m2 a2cos2ϕ− β ) ˆ h + 1 a2 d2ˆh dϕ2 − tan ϕ a2 dˆh dϕ = ˆR. (26)
Now, the fourth-order compact, the sixth- and eighth-order SCFDM and CCFDM methods are used to discretize the first and second latitudinal derivatives in equation (26). We arrive at the following block tridiagonal system of equations
ASjUSj−1+ BSjUSj + CSjUSj+1= DSj (27)
for the different methods are UC4j = ˆ hj ˆ h(1)j ˆ h(2)j , USC6 j = ˆ hj ˆ h(1)j ˆ h(2)j ˆ h(3)j ˆ h(4)j , USC8j = ˆ hj ˆ h(1)j ˆ h(2)j ˆ h(3)j ˆ h(4)j ˆ h(5)j ˆ h(6)j , UCC6j = ˆ hj ˆ h(1)j ˆ h(2)j , UCC8 j = ˆ hj ˆ h(1)j ˆ h(2)j ˆ h(3)j . (28) The superscripts C4, SC6, SC8, CC6, CC8 denote the fourth-order compact, the sixth-order SCFDM, the eighth-order SCFDM, the sixth-order CCFDM and the eighth-order CCFDM methods, respectively. Equation (27) is a three-point implicit formulation which solution leads to a linear block tridiagonal system of equations. The block size for the fourth-order compact is 3× 3, for the sixth-order SCFDM is 5 × 5, for the eighth-order SCFDM is 7 × 7, for the sixth-order CCFDM is 3× 3 and for the eighth-order CCFDM is 4 × 4.
To close the system at south pole (j = 1) and north pole (j = Nϕ) boundary conditions are
needed. The boundary values are obtained by using the symmetry constraint of the Fourier coefficients (e.g., [33]). For the grid system used in the present work, the symmetry constraints imply that at j = 1 we have ˆh(k)0 = (−1)k(−1)mˆh(k)1 and for j = Nϕˆh
(k)
Nϕ+1= (−1)
k(−1)mˆh(k)
Nϕ.
The superscript (k) denotes as before the kth derivative (k = 0, 1, 3,· · · ) and m denotes the wave number. There is one sign change due to odd and even wave numbers and another sign change related to the spatial derivative.
By inserting these values into equation (27) at boundaries and solving the resulting block tridiagonal system of equations, the values of ˆh are obtained and then by using FFT, the values
of ¯hi,j are found. The details of matrix coefficients ASj, BSj, CSj and right hand side vector DSj
are given in appendix I.
5.2. The Poisson equation
The Poisson equations in (7) have the following form
∇2ψ = ζ. (29)
The procedure for the numerical solution of the Poisson equations is similar to the procedure for the modified Helmholtz equation above and can be reached by letting β = 0 in equation (26). For nonzero wave numbers (m ̸= 0) equation (29) can be solved in a similar manner as the modified Helmholtz equation.
The procedure for the zeroth Fourier mode (m = 0) is different. It can be seen that by using
β = 0 and m = 0 in (26), ˆh vanishes and the Poisson equation for wave number m = 0 is
singular, since the solution is unique within a constant. In fact, a Poisson equation on a sphere must satisfy the following compatibility condition (e.g.,[35, 33])
∫ 2π 0
∫ π/2 −π/2
ζ cos ϕdϕdλ = 0 (30)
By substitution of the Fourier components in equation (29), and considering m = 0 we have 1 a2 d2ψ¯ dϕ2 − tan ϕ a2 d ¯ψ dϕ = ¯ζ (31)
in which ¯ψ = ˆψ(m = 0) and ¯ζ = ˆζ(m = 0). After application of the Fourier method in longitude
direction the equivalent form of the compatibility condition (30) for the latitudinal direction
is ∫
π/2 −π/2
¯
ζ cos ϕdϕ = 0. (32)
To impose the compatibility condition, first the fourth-order compact, the sixth- and eighth-order SCFDM and the sixth- and eighth-eighth-order CCFDM methods are used to discretize the first and second latitudinal derivatives in equation (31) in the same way as for (26), and then the right hand side of equation (31) in each discrete version is replaced by ¯ζj− P where P is a
constant. P can be found by using the discrete version of the compatibility condition (32) as
P =
j=N∑ϕ
j=1
¯
ζjcos ϕj. (33)
Now, by finding the values of ¯ψj for the zeroth Fourier mode and having the values of ˆψj
for other wave numbers we employ the FFT to find the solution of the Poisson equation. The matrix coefficients and right hand side vector for different compact methods used to discretize equation (31) are also described in appendix I.
6. The spatial filter
The generation of fine-scale vortical structures by the nonlinear advection terms is a dominant feature of the shallow water equations. A problem facing any grid-based numerical method developed for the shallow water equations is misinterpretation of the fine-scale structures due to insufficient resolution and the potential for non-linear numerical instability. In this work, a simple spatial filter is developed to control the non-linear instability. In addition, on a regular latitude-longitude system the lines of constant longitude converge at the poles which leads to a high concentration of grid points near poles. The cluster of grid points forces a short time step to satisfy the CFL stability condition and this is often referred to as the polar problem [36]. The spatial filter can act as a smoother to avoid the polar problem.
The general form of a one-dimensional spatial filter (for example in latitude, ϕ) to filter an arbitrary field, Υ, can be written as
˜ Υj = Υj+ σp(∆ϕ2)p ( ∂2Υ ∂ϕ2 )p j , p = 1, 2, 3, . . . (34)
where ˜Υj denotes the filtered values. This formulation is used to construct a low-pass filter
which removes wave-lengths shorter than 2∆ϕ. Using higher values of p moves the filtering effect towards the shorter waves. The order of accuracy of the filter is determined by the finite difference method used to discretize the even derivative in equation (34).
The order of even derivatives (p) and the discretization method in equation (34) can be selected by the requirements that (i) the dissipation is kept minimal for stable solutions and (ii) the desired order of accuracy is not degraded by the action of the filter.
The unknown coefficient σp is determined by setting the location of cutoff in the transfer
function of the filter (e.g., [37]). The transfer function, Tf, of the filter is found using Fourier
transformation of equation (34) and is defined as ˜Υ(eiκϕ) = Tf(κ)eiκϕ, where i =
√ −1 and κ denotes wave number. For a low-pass filter which removes the wave-lengths shorter than
2∆ϕ it is necessary to have Tf(κmax) = 0 where κmax = π/∆ϕ. Therefore, the spatial filter coefficient can be found as
σp= ( −1 T2(κmax) )p (35) where T2is the transfer function of the second derivative. Table I presents values of σp for the
various schemes used in this study.
Table I. values of the spatial filter coefficient, σp, for different methods.
Method σp Fourth-order compact (−16 )p Sixth-order SCFDM (7.05882353−1 )p Sixth-order CCFDM (−19.6)p Eighth-order SCFDM (7.67741935−1 )p Eighth-order CCFDM (9.84615385−1 )p
The same spatial differencing method that is applied to discretize the governing equations is used to estimate the even derivative in the spatial filter (34). In addition, the spatial filter is sequentially applied in latitude and longitude directions at each time step of the time integration. We have used p = 3 for the filter power in both latitude and longitude directions. However, for the longitude direction we have used a slightly more diffusive filter with p = 2 near the poles where 0.7π2 ≤ ϕ ≤ π2 and−π2 ≤ ϕ ≤ −0.7π2 , to better remove the restriction of the polar problem on the time step.
7. The order of convergence
The spatial order of convergence of high-order compact schemes used in this work is assessed by considering a Helmholtz equation with a known analytical solution. In fact, at a fixed time the main part of the total algorithm consists of the solution of the elliptic equations and it is therefore a meaningful test.
The following elliptic equation is solved
∇2Ψ− aΨ = R
Ψ (36)
where RΨ= 2ν sin ϕ− 30K sin ϕ cos4ϕ cos 4λ + a3ν sin ϕ− a3K cos4ϕ sin ϕ cos 4λ and a is the radius of the Earth. The constants ν and K are defined as ν = K = 7.848× 10−6s−1.
The Rossby-Haurwitz wave [28], is the analytical solution of equation (36) given by Ψ =−a2ν sin ϕ + a2K cos4ϕ sin ϕ cos 4λ. (37) To investigate the order of convergence, the numerical solution of equation (36) generated by the 4th-order compact, 6th-order SCFDM, 6th-order CCFDM, 8th-order SCFDM and 8th-order CCFDM methods is compared with the analytical solution (37). The following normalized global error is used for the error measurement [28]:
l2(Ψ) = { I(|Ψn(λ, ϕ)− ΨE(λ, ϕ)| 2)}1/2 {I (|ΨE(λ, ϕ)| 2)}1/2 (38)
where Ψnand ΨE are the numerical and exact solutions, respectively. The function I is defined
as I(Ψ) = 1 4π Nλ ∑ k=1 Nϕ ∑ l=1 Ψk,lcos ϕ∆ϕ∆λ (39)
in which Nλ and Nϕ are the number of grid points in longitude and latitude directions and
Ψk,l is the argument at grid point (λk, ϕl).
For the numerical solution of equation (36), the compact finite difference methods are only applied for discretization in the latitudinal direction (see section 5). Therefore, to find the order of convergence we use ∂Ψ/∂ϕ instead of Ψ. Figure 1 presents the normalized global error l2 for ∂Ψ/∂ϕ at different uniform horizontal resolutions in longitude and latitude directions calculated by the 4th-order compact, 6th-order SCFDM, 6th-order CCFDM, 8th-order SCFDM and 8th-8th-order CCFDM schemes. Figure 1 clearly shows that the convergence rates of the different solutions are in agreement with theoretical order of convergence. In addition, it is seen that the 6th-order CCFDM method produces less error than the 6th-order SCFDM method and the 8th-order CCFDM generates less error than the 8th-order SCFDM method, respectively.
8. Numerical results
8.1. Initial condition
The test case proposed in reference [38] which is an unstable barotropic mid-latitude zonal jet is used as the starting point for the numerical solution and evaluation of the spherical shallow water equations. The initial vorticity filed of the unstable barotropic zonal jet is
ζ = 1
a cos ϕ
(
u sin ϕ− cos ϕ∂u ∂ϕ
)
. (40)
For this test case v = 0 and δ = 0 and
u = 0 for ϕ≤ ϕ0 umax en ( 1 (ϕ− ϕ0)(ϕ− ϕ1) ) for ϕ0< ϕ < ϕ1 0 for ϕ≥ ϕ1 (41)
Nλ ( = 2 Nφ ) L2 ( E rr o r) 100 200 300 400 500 10-16 10-14 10-12 10-10 10-8 10-6 10-4 10-2 100 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM P=4 P=6 P=8
Figure 1. Normalized l2 error of ∂Ψ/∂ϕ field versus number of grid points obtained using 4th-order compact method, 6th-order SCFDM, 6th-order CCFDM, 8th-order SCFDM and 8th-order CCFDM. Lines with slops 4, 6 and 8 are also shown. A uniform grid resolution is used in both λ and ϕ directions.
is used in which the constants are
umax= 80 ms−1, ϕ0= π 7, ϕ1= π 2 − ϕ0, en= exp[ −4 (ϕ1− ϕ0)2 ]. The following relation is used to find the initial height field
h = 1 g { gh0− ∫ ϕ au(ϕ′) [ f +tan(ϕ ′) a u(ϕ ′)]dϕ′ } (42) where the value of constant h0needs to be found such that the domain area average of depth is H = 10000 m. Here, we have used the trapezoid rule to evaluate the integral in equation (42)
numerically and finding h to machine precision.
Next, the following perturbation is added to the basic height in equation (42) to trigger the barotropic instability
h′(λ, ϕ) = hpcos(ϕ)e−(λ/γ)
2
e−[(ϕ2−ϕ)/ ¯β]2, for − π < λ < π (43)
where λ is longitude and the constants are hp= 120 m, γ = 1/3, ϕ2= π/4 and ¯β = 1/15.
8.2. Grid resolution and time step
The grid resolutions used in the numerical experiments reported here are Nλ × Nϕ =
128×64, 256×128, 512×256. To find the time steps a fixed Courant number of less than unity based on gravity-wave speed is used. Using small time steps reduces the frequency distortion
of the semi-implicit time-stepping scheme on gravity waves. The time steps ∆t = 270, 60, 15 seconds have been used for the successive grids.
8.3. Results and discussions
We begin by providing some qualitative results for the vorticity field. Figure 2 presents the time evolution of the vorticity field for Nλ× Nϕ= 512× 256 resolution at t = 0, 2, 4, 6 days (1 day
= 86400 s) calculated by the eighth-order CCFDM. A qualitative comparison of these results with those presented in [38, 17] indicates the validity of computations. In addition, figures 3 and 4 provide a qualitative comparison of the vorticity field calculated by the compact schemes used in this study at time t = 6 for Nλ×Nϕ= 256×128 and Nλ×Nϕ= 512×256 resolutions.
It can also be seen that the sixth-order and eighth-order CCFDM methods produce better results than the sixth-order and eighth-order SCFDM methods. Furthermore, it is evident that the fourth-order compact method can not reproduce the structure of flow in particular the high gradients of vorticity field at Nλ× Nϕ = 256× 128 resolution. It can be seen that
the solutions generated by the sixth-order and eighth-order SCFDM and CCFDM methods are very similar and a quantitative assessment of the solutions is needed to better inspect and understand the properties of different methods.
Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (t = 2) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (t = 0) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (t = 6) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (t = 4)
Figure 2. The time evolution of the vorticity field at 2-day time intervals for unstable barotropic jet generated by the eighth-order CCFDM. The grid resolution is Nλ× Nϕ = 512× 256 and contour
interval is 1× 10−5 s−1. Solid and dashed lines are for positive and negative contours, respectively. The zero contour is not shown and only the Northern Hemisphere is shown in figure. Figure 5 presents time evolution of magnitude of the maximum vorticity gradient for different
Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (a) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (c) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (b) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (d) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (e)
Figure 3. The vorticity field at time t = 6 days, for (a) 4th-order compact method, (b) 6th-order SCFDM, (c) 6th-order CCFDM, (d) 8th-order SCFDM and (e) 8th-order CCFDM. The grid resolution is Nλ×Nϕ= 256×128 and the contour interval is 1×10−5s−1. Solid and dashed lines are for positive
and negative contours, respectively.
methods. Clearly, the sixth-order and eighth-order CCFDM exhibit more levels of the vorticity gradients than the other methods. As sharp gradients of vorticity field are a feature of the unstable barotropic mid-latitude zonal jet test case used in this study [38, 17], it can be seen that sixth-order and eighth-order CCFDM represent the vorticity field in a better way than the other methods.
Figure 6 shows the power spectrum of total energy, which is the spectrum of 12h(u2+ v2) + 1
2gh
2, at time t = 6 days for different methods and different resolutions. It is seen that the spectrum of total energy for different methods can be differentiated only at the lower end of the scales. The theory in [1] suggests that at inertial range, the spectrum of the energy decays as
Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (a) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (c) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (b) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (d) Longitude L a ti tu d e -150 -100 -50 0 50 100 150 0 20 40 60 80 (e)
Figure 4. The vorticity field at time t = 6 days, for (a) 4th-order compact method, (b) 6th-order SCFDM, (c) 6th-order CCFDM, (d) 8th-order SCFDM and (e) 8th-order CCFDM. The grid resolution is Nλ×Nϕ= 512×256 and the contour interval is 1×10−5s−1. Solid and dashed lines are for positive
and negative contours, respectively.
K−3, K being the total wave number. Figure 6 indicates that at the inertial range all methods used in this study exhibit this feature of the solution.
By inspecting the time evolution of potential enstrophy C2 and potential palinstrophy P defined by C2= 1 2 ⟨ (1 + ˜h)Q2 ⟩ , P = 1 2 ⟨ (1 + ˜h)|∇Q| ⟩ (44) we can compare the global spatial distribution of potential vorticity (Q = (ζ +f )/(1+˜h)) in the
solutions. In equation (44) the domain area average is denoted by⟨ ⟩. The time variations of the percentage changes in C2′ and P , i.e. (C2′(t)− C2′(0)) /C2′(0) and (P (t)− P (0)) /P (0) multiplied
Time (day) M a x im u m v o rt ic it y g ra d ie n t 0 1 2 3 4 5 6 0.001 0.002 0.003 0.004 0.005 0.006 0.007 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (a) Time (day) M a x im u m v o rt ic it y g ra d ie n t 0 1 2 3 4 5 6 0.001 0.002 0.003 0.004 0.005 0.006 0.007 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (b) Time (day) M a x im u m v o rt ic it y g ra d ie n t 0 1 2 3 4 5 6 0.002 0.003 0.004 0.005 0.006 0.007 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (c)
Figure 5. Time evolution of maximum vorticity gradient for different methods at grid resolutions (a) Nλ× Nϕ= 128× 64, (b) Nλ× Nϕ= 256× 128 and (c) Nλ× Nϕ= 512× 256
by 100, where C2′ = C2− ⟨
f2⟩/2 are shown in figures 7 and 8 for different methods and different grid resolutions. We have subtracted the time-independent contribution of⟨f2⟩/2 from C
2 to have a better quantification of the effects of damping/discretization error in destroying the global conservation of potential enstrophy.
Figures 7 and 8 show that the sixth-order and eighth-order CCFDM methods have better global conservation properties of the potential enstrophy than the other methods. For the sixth-order and eighth-order CCFDM methods at Nλ× Nϕ= 512× 256 resolution, a less than
1% reduction in C2′ during the first 3 days is followed by a much stronger loss of potential enstrophy to near 5% at time t = 6. For the sixth-order and eighth-order SCFDM methods at
Nλ× Nϕ= 512× 256 resolution, it is observed that a less than 2% reduction in C2′ during the first 3 days is followed by a stronger loss of the potential enstrophy to near 7% at t = 6. For
Log (wavenumber) E n e rg y p o w e r s p e c tr u m 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 30 35 40 45 50 Theory 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (a) Log (wavenumber) E n e rg y p o w e r s p e c tr u m 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 30 35 40 45 50 Theory 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (b) Log (wavenumber) E n e rg y p o w e r s p e c tr u m 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 30 35 40 45 50 Theory 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (c)
Figure 6. Power spectrum of total energy at time t = 6 day for (a) Nλ× Nϕ = 128× 64, (b)
Nλ× Nϕ= 256× 128 and (c) Nλ× Nϕ= 512× 256 resolutions. In each panel the straight line plots
a K−3spectrum, for comparison.
the fourth-order compact method at Nλ×Nϕ= 512×256 resolution, a less than 3% reduction
in C2′ during the first 3 days is followed by a much stronger loss of the potential enstrophy to near 9% in at time t = 6. It can be seen that a slightly better conservation of potential enstrophy is obtained when the order of spatial differencing of each scheme is increased.
For potential palinstrophy, as the grid resolution increases and the solution converges, a rapid growth of disturbances in the unstable jet which reaches a maximum at time t = 6 days can be observed. This is in agreement with the results presented in figure 5. Furthermore, it is observed that a slightly higher values of potential palinstrophy are obtained when the order of spatial differencing of each scheme is increased.
V1 P o te n ti a l e n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (a) V1 P o te n ti a l e n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (b) V1 P o te n ti a l e n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (c)
Figure 7. Time evolution of the percentage of relative change in potential enstrophy for different methods at (a) Nλ× Nϕ = 128× 64, (b) Nλ× Nϕ = 256× 128 and (c) Nλ× Nϕ = 512× 256
resolutions. defined by E = ⟨ 1 2h(u 2+ v2) +1 2gh 2 ⟩ (45) where h = H(1 + ˜h) and ⟨⟩ denotes the domain area average. Figure 9 shows a 20 days time
evolution of the percentage changes in E′defined by (E′(t)− E′(0)) /E′(0) multiplied by 100, where E′ = E−⟨c2H⟩/2. The time-independent part⟨c2H⟩/2 has been subtracted from E to have a better representation of changes in total energy. It is seen that the sixth-order and eighth-order CCFDM methods at time t = 20 days shows approximately 7 percent reduction in total energy, the eighth-order SCFDM method exhibits approximately 9 percent reduction in total energy, the sixth-order SCFDM method shows approximately 10 percent reduction in
V1 P a li n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -80 -60 -40 -20 0 20 40 60 80 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (a) V1 P a li n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -80 -60 -40 -20 0 20 40 60 80 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (b) V1 P a li n s tr o p h y ( re la ti v e ) 0 1 2 3 4 5 6 -80 -60 -40 -20 0 20 40 60 80 4th-order compact 6th-order CCFDM 6th-order SCFDM 8th-order CCFDM 8th-order SCFDM (c)
Figure 8. Time evolution of the percentage of relative change in potential palinstrophy for different methods at (a) Nλ×Nϕ= 128×64, (b) Nλ×Nϕ= 256×128 and (c) Nλ×Nϕ= 512×256 resolutions.
total energy and the fourth-order compact method shows approximately 13 percent reduction in total energy. We conclude that the sixth-order and eighth-order CCFDM methods better conserve total energy compared to the other methods.
To summarize, the qualitative and quantitative assessment of the results indicate noticeable improvement of the results by using the sixth-order and the eighth-order SCFDM and CCFDM schemes over the fourth-order compact method in terms of global conservation and solution quality.
Time (day) C h a n g e i n t o ta l e n e rg y ( re la ti v e ) 0 5 10 15 20 25 -15 -10 -5 0 5 4th-order Compact 6th-order SCFDM 6th-order CCFDM 8th-order SCFDM 8th-order CCFDM
Figure 9. Time evolution of the percentage of relative change in total energy for different methods at Nλ× Nϕ= 256× 128 resolution.
9. Conclusions
The focus of this work was on the application and qualitative and quantitative comparison of two families of high-order compact finite difference methods, namely, the super compact finite difference and the combined compact finite difference methods, to spatial differencing of the spherical shallow water equations in terms of vorticity, divergence and height on a regular latitude-longitude grid. They were compared with the standard fourth-order compact scheme. In addition, a semi-implicit Runge-Kutta method was developed for time advancing of the vorticity, divergence and height representation of the spherical shallow water equations.
Furthermore, to control the non-linear instability and avoiding the polar problem a spatial filter was developed. The assessment of the five high-order compact schemes was carried out for the test case developed in reference [38]. For numerical solution of the elliptic equations appeared in formulation of the problem the direct procedure proposed in references [33, 34] with minor changes was extended to the high-order compact finite difference schemes.
A study of the convergence rate showed that all spatial schemes convergence in agreement with the formal order of accuracy of the method. For qualitative assessment of the solutions the measurements of maximum vorticity gradient, power spectrum of total energy, relative change in potential enstrophy and potential palinstrophy were used.
Qualitative and quantitative measurements for the test case employed in this work, shows that using the the sixth-order and eighth-order CCFDM and the sixth-order and eighth-order SCFDM for spatial differencing of the spherical shallow water equations lead to a noticeable improvement of the solution over the fourth-order compact method. However, the performance of sixth-order and order CCFDM methods is better than the sixth-order and
eighth-order SCFDM methods.
Furthermore, the computational cost of the spatial differencing for finite difference methods used in this study is linearly proportional to the number of grid points, therefore, the methods in terms of computational cost from the lowest to highest are: the fourth-order compact, the sixth-order CCFDM, the order CCFDM, the sixth-order SCFDM and the eighth-order SCFDM. Therefore, in terms of qualitative and quantitative measurements including computational cost we conclude that the sixth-order and eighth-order CCFDM methods are superior to the sixth-order and eighth-order SCFDM methods for spatial differencing of the spherical shallow water equations.
ACKNOWLEDGEMENTS
Authors would like to thank university of Tehran, Link¨oping university and National Supercomputer Center in Sweden (NSC) for supporting this research work.
APPENDIX
I. Details of block tridiagonal system for numerical solution of elliptic equations The matrix coefficients ASj, BSj, CSj and right hand side vector DSj of the block tridiagonal system
obtained from different compact method used in this study are presented here. I.1. Helmholtz equation
The matrix coefficients for the fourth-order compact method for 2 < j < Nϕ−1 are
BC4j = α0j − tan ϕ2 j 1 3 0 2 d2 0 5 6 , AC4 j = 01 0 0 2d 1 6 0 −1 d2 0 1 12 , CC4 j = −10 0 0 2d 1 6 0 −1 d2 0 1 12 where d is grid space and
αj= −m 2 cos2ϕ j − a 2 βj (46)
and coefficients for j = 1 and j = Nϕare
BC41 = α1 − tan ϕ1 1 (−1)m 2d 2 3+ (−1)m+1 6 0 −(−1)m d2 +d22 0 (−1)m 12 + 5 6 , CC4 1 = −10 0 0 2d 1 6 0 −1 d2 0 121 BC4Nϕ= αNϕ − tan ϕNϕ 1 −(−1)m 2d 2 3 + (−1)m+1 6 0 −(−1)m d2 + 2 d2 0 −(−1) m 12 + 5 6 , AC4 Nϕ = 01 0 0 2d 1 6 0 −1 d2 0 1 12 and the right hand side vector Dj for 1 < j < Nϕis
DC4j = a 2Rˆ j 0 0
The matrix coefficients for the sixth-order SCFDM method for 2 < j < Nϕ−1are BSC6j = αj − tan ϕj 1 0 0 0 −2d 0 −103d3 0 0 −2 0 −56d2 0 −2 0 −d2 0 −907d4 0 0 −2 0 −56d2 , A SC6 j = 0 0 0 0 0 −1 0 0 −d3 60 0 0 1 0 −d122 0 1 0 0 0 −d4 360 0 0 1 0 −d2 12 CSC6j = 0 0 0 0 0 1 0 0 −d3 60 0 0 1 0 −d2 12 0 1 0 0 0 −360d4 0 0 1 0 −d2 12 and coefficients for j = 1 and j = Nϕare
BSC61 = α1 − tan ϕ1 1 0 0 −(−1)m −2d 0 −3d103−(−1)m+160 d3 0 0 −2 + (−1)m+1 0 −5d62−(−1)m+112 d2 0 −2 + (−1)m 0 −d2 0 −7d4 90 − (−1)md4 360 0 0 −2 + (−1)m 0 −5d62 −(−1)12md2 CSC61 = 0 0 0 0 0 1 0 0 −d3 60 0 0 1 0 −d122 0 1 0 0 0 −d4 360 0 0 1 0 −d2 12 BSC6Nϕ = αNϕ − tan ϕNϕ 1 0 0 (−1)m −2d 0 −3d103−(−1)m+160 d3 0 0 −2 + (−1)m+1 0 −5d62−(−1)m+112 d2 0 −2 + (−1)m 0 −d2 0 −7d4 90 − (−1)md4 360 0 0 −2 + (−1)m 0 −5d62 −(−1)12md2 ASC6Nϕ = 0 0 0 0 0 −1 0 0 −d3 60 0 0 1 0 −d2 12 0 1 0 0 0 −d4 360 0 0 1 0 −d122 and the right hand side vector Dj for 1 < j < Nϕis
DSC6j = a2Rˆj 0 0 0 0
The matrix coefficients for the eighth-order SCFDM method for 2 < j < Nϕ−1are BSC8j = αj − tan ϕj 1 0 0 0 0 0 −2d 0 −d33 0 −2d1265 0 0 d 0 d3 2 0 7d5 180 0 0 0 0 d3 0 5d125 0 −2 0 −d2 0 −d4 12 0 −6d 6 2240 0 0 d2 0 d24 0 7d1806 0 0 0 0 d4 0 5d126 ASC8j = 0 0 0 0 0 0 0 −1 0 0 0 0 −2d50405 0 0 −d2 0 0 0 720d5 0 0 0 0 −d23 0 d5 24 0 1 0 0 0 0 0 40320−2d6 0 0 −d22 0 0 0 720d6 0 0 0 0 −d24 0 d6 24 CSC8j = 0 0 0 0 0 0 0 1 0 0 0 0 −2d50405 0 0 −d2 0 0 0 720d5 0 0 0 0 −d23 0 d5 24 0 1 0 0 0 0 0 40320−2d6 0 0 −d22 0 0 0 720d6 0 0 0 0 −d24 0 d6 24 and coefficients for j = 1 and j = Nϕare
BSC81 = α1 − tan ϕ1 1 0 0 0 0 −(−1)m −2d 0 −d3 3 0 −2d5126 − (−1)m+12d5 5040 0 0 d−(−1)m+1d2 0 d32 0 7d5180 +(−1)m+1d5720 0 0 0 0 d3−(−1)m+1d3 2 0 5d512 + (−1)m+1 24 0 −2 + (−1)m 0 −d2 0 −d4 12 0 −6d62240 − (−1)m2d6 40320 0 0 d2−(−1)md2 2 0 d42 0 7d6180+ (−1)md6 720 0 0 0 0 d4−(−1)md42 0 5d612 +(−1)m24 CSC81 = 0 0 0 0 0 0 0 1 0 0 0 0 −2d50405 0 0 −d2 0 0 0 d5 720 0 0 0 0 −d23 0 d245 0 1 0 0 0 0 0 40320−2d6 0 0 −d22 0 0 0 d6 720 0 0 0 0 −d24 0 d246 BSC8Nϕ = αNϕ − tan ϕNϕ 1 0 0 0 0 (−1)m −2d 0 −d33 0 −2d5126 −(−1)m+12d55040 0 0 d−(−1)m+1d 2 0 d32 0 7d5180 + (−1)m+1d5 720 0 0 0 0 d3−(−1)m+1d3 2 0 5d5 12 + (−1)m+1 24 0 −2 + (−1)m 0 −d2 0 −d4 12 0 −6d62240 − (−1)m2d6 40320 0 0 d2−(−1)md2 2 0 d4 2 0 7d6 180+ (−1)md6 720 0 0 0 0 d4−(−1)md4 2 0 5d612 + (−1)m 24
ASC8Nϕ = 0 0 0 0 0 0 0 −1 0 0 0 0 −2d50405 0 0 −d2 0 0 0 720d5 0 0 0 0 −d23 0 d5 24 0 1 0 0 0 0 0 40320−2d6 0 0 −d22 0 0 0 720d6 0 0 0 0 −d24 0 d6 24 and the right hand side vector Dj for 1 < j < Nϕis
DSC8j = a2Rˆj 0 0 0 0 0 0
The matrix coefficients for the sixth-order CCFDM method for 2 < j < Nϕ−1 are
BCC6j = α0j − tan ϕ1 j 10 6 d2 0 1 , ACC6 j = 150 0 0 16d 7 16 d 16 −3 d2 −98d −18 , CCC6 j = −150 0 0 16d 7 16 −d 16 −3 d2 9 8d −18 and coefficients for j = 1 and j = Nϕare
BCC61 = α1 − tan ϕ1 1 (−1)m15 16d 1 + (−1)m+17 16 (−1)md 16 6 d2 − (−1)m3 d2 −(−1) m+19 8d 1− (−1)m 8 , CCC6 1 = −150 0 0 16d 7 16 −d16 −3 d2 9 8d −1 8 BCC6Nϕ = αNϕ − tan ϕNϕ 1 −(−1)m15 16d 1 + (−1)m+17 16 −(−1)md 16 6 d2 − (−1)m3 d2 (−1)m+19 8d 1− (−1)m 8 , ACC6 Nϕ = 150 0 0 16d 7 16 d 16 −3 d2 −98d −18 and the right hand side vector for 1 < j < Nϕis
DCC6j = a 2Rˆ j 0 0
The matrix coefficients for the eighth-order CCFDM method for 2 < j < Nϕ−1are
BCC8j = αj − tan ϕj 1 0 0 1 0 0 8 d2 0 1 0 0 0 0 1 , ACC8 j = 0 0 0 0 35 32d 19 32 d 8 d2 96 −4 d2 −2916d −516 −d48 −105 16d3 −10516d2 −158d −3 16 CCC8j = 0 0 0 0 −35 32d 19 32 −d 8 d2 96 −4 d2 29 16d −516 d 48 105 16d3 −10516d2 15 8d −316
and coefficients for j = 1 and j = Nϕare BCC81 = α1 − tan ϕ1 1 0 (−1)m35 32d 1 + (−1)m+119 32 (−1)md 8 (−1)m+1d2 96 8 d2 − (−1)m4 d2 −(−1) m+129 16d 1− (−1)m5 16 −(−1)m+1d 48 −(−1)m105 16d3 −(−1) m+1105 16d2 −(−1) m15 8d 1− (−1)m+13 16 , C CC8 1 = 0 0 0 0 −35 32d 19 32 −d 8 d2 96 −4 d2 29 16d −516 d 48 105 16d3 −10516d2 15 8d −316 BCC8Nϕ = αNϕ − tan ϕNϕ 1 0 −(−1)m35 32d 1 + (−1)m+119 32 −(−1)md 8 (−1)m+1d2 96 8 d2 − (−1)m4 d2 (−1)m+129 16d 1− (−1)m5 16 (−1)m+1d 48 (−1)m105 16d3 −(−1) m+1105 16d2 (−1)m15 8d 1− (−1)m+13 16 , A CC8 Nϕ = 0 0 0 0 35 32d 19 32 d 8 d2 96 −4 d2 −2916d −516 −d48 −105 16d3 −10516d2 −158d −316 and the right hand side vector for 1 < j < Nϕis
DCC8j = a2Rˆj 0 0 0
I.2. Poisson equation
The matrix coefficients for the Poisson equation for nonzero wave numbers (m ̸= 0) are found by setting βj = 0 in (46). For wave number zero (m = 0) equation (31) is used and we need to modify
the first row of BSj and DSj.
REFERENCES 1. Pedlosky J. Geophysical fluid dynamics. Springer-Verlag, 1987.
2. Vallis G.K. Atmospheric and oceanic fluid dynamics: Fundamentals and large-scale circulation. Cambridge University Press, 2006.
3. Kreiss H. O., Oliger J. Comparison of accurate methods for the integration of hyperbolic equations. Tellus 1972; 24:199–215.
4. Hirsch S. R. Higher order accurate difference solutions of fluid mechanics problems by a compact differencing technique. J. Comput. Phys. 1975; 19:90–109.
5. Lele S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992; 103:16–42. 6. Kreiss H.O., Scherer G. Finite element and finite difference methods for hyperbolic partial differential equations, in: C. De Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equation, Academic Press, New York, 1974.
7. Carpenter M.H., Nordstr¨om J., Gottlieb D. A stable and conservative interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 1999; 148:341-365.
8. Mattsson K., Nordstr¨om J. Summation by parts operators for finite difference approximations of second derivatives, J. Comput. Phys. 2004; 199:503-540.
9. Strand B. Summation by parts for finite difference approximation for d/dx, J. Comput. Phys. 1994; 110:47-67.
10. Durran, D.R. Numerical methods for fluid dynamics with applications to geophysics. Second Edition, Springer, 2010.
11. Navon I.M., Riphagen H.A. An implicit compact fourth-order algorithm for solving the shallow water equations in conservative-law form. Mon. Wea. Rev. 1979;107:1107–1127.
12. Chang H.H., Shirer N., Compact spatial differencing techniques in numerical modeling. Mon. Wea. Reu. 1985; 113:409–423.
13. Spotz W.F., Taylor M.A., Swarztrauber P.N. Fast shallow-water solvers in latitude-longitude coordinates. Journal of Computational Physics, 1998; 145:432–444.
14. Blayo E. Compact finite difference schemes for ocean models: 1.ocean waves. Journal of Computational Physics 2000; 164:241–257.
