Development of a potential vorticity based dynamical core for general circulation models using the diabatic contour-advective semi-Lagrangian algorithm

The development of the dynamical core of a potential-vorticity-based atmospheric general circulation model is explored. There are some advantages of using potential vorticity (PV) as a prognostic variable in that the resulting model can give more accurate simulation of the evolution of PV, as arguably the most fundamental dynamical quantity. Further, there is possibility of explicit representation of unbalanced part of the flow during time integration, though in an approximate manner, by making proper choice of the prognostic variables used alongside PV. In this way, the model is equipped with some built-in form of the balance relation for PV inversion, which helps to maintain the underlying balance. A closed set of equations is constructed using the variables , a PV-like variable described below, and which are, respectively, the horizontal velocity divergence and an approximate form of horizontal acceleration divergence. For the primitive equations linearized around a rest state, it can be shown that there is a direct correspondence between the Rossby modes and the Q variable, on the one hand, and between the inertia-gravity modes and the and  variables, on the other hand. Linearizing the primitive equations in the generalized vertical coordinate around a resting basic state, the symmetric matrix  relating the column vector of the time tendency of modified pressure to the column vector of horizontal divergence is found. Here, the modified pressure is defined by  with and respectively, the perturbation geopotential, temperature and potential temperature, specific heat capacity at constant pressure and the basic state Exner function. The eigenvectors of are used to define the vertical modes and the projection of any given column vector from the physical space to vertical mode space and vice versa. This facilitates to generalize the Boussinesq PV-based multi-layer primitive-equation models to the corresponding non-Boussinesq set of equations. A PV-like quantity is defined by in which is the Coriolis parameter, the relative vertical vorticity, and the normalized perturbation pressure thickness. Here and are, respectively, the perturbation and the basic state pressure. The variable becomes the same as Rossby–Ertel PV whenever coincides with. Further, with the definition of modified pressure given above, the variable becomes equal to with the northward gradient of When coincides withbecomes equal to acceleration divergence. To use the variables and as the prognostic variables, one has to implement an inversion procedure to obtain the velocity field and the thermodynamic variables at each time step. Making use of the definition of and and projecting onto the vertical-mode space results in a modified Helmholtz equation for which is solved by spectral transform in longitude and fourth-order compact in latitude following the procedure introduced by Mohebalhojeh and Dritschel in 2007. Solving for the modified pressure can then be obtained either directly through the matrix relation  or through projection onto vertical-mode space. The task is then to find the thermodynamic variables using the information available for at each column of fluid.
The PV as a determining variable for vortical flows is given the highest priority in terms of accuracy. For this purpose, the Contour-Advective Semi-Lagrangian (CASL) algorithm, previously implemented for various settings and models including the many-layer Boussinesq primitive equation models on the sphere, provides the natural choice. An extension of CASL called DCASL has already been applied to the thermally-forced shallow water equations (SWEs) on the sphere. In generalized vertical coordinate, the evolution equation of  is similar to that of PV in the thermally-forced SWEs. Therefore, the available DCASL can be generalized for the non-Boussinesq equations with little effort.   
The generalized vertical coordinate is set as  with defined in such a way as to increase monotonically with geometrical height from zero at the surface to one at the top level. The functions and g are determined in such a way that (i) tends to  and  when pressure  tends to its value at the surface and the top of the model, respectively, and (ii) the condition  is satisfied to ensure monotonicity whenever and  where  and  are prescribed values of the lowest value of potential temperature and the vertical gradient of potential temperature with respect to sigma, respectively.
The time evolution of a two-layer baroclinically unstable midlatitude jet over a 30-day period is investigated as a test case to examine the performance of the algorithm developed. It should be mentioned that various experiments using different basic-state structures have been carried out. The experiment reported is however for the one with a uniform stratification obtained by setting a constant lapse rate of from to  This choice of the basic-state structure leads to a flow regime with order one Rossby and Froude numbers. Results show the formation and development of an intense baroclinic wave with zonal wave number 3. Further, embedded in the baroclinic wave there are inertia-gravity waves generated by vortical flow in a manner resembling what has previously observed for the Boussinesq primitive-equation model. The successful integration of model in extreme flows gives us confidence to further develop the algorithm to a dynamical code for atmospheric general circulation models.